Question for making count table for edgeR

[lucer105]

I am trying to make a count table for following edgeR..... but after running the edgeR, my marker gene didn't change as it suppose to change. (It suppose to change, confirmed by qPCR, and cuffdiff, and load .bam file in IGV, and it is huge change). I checked the count table content I made and the reads in the table is weird. I got 0, 0 reads for control and 1, 7 reads for treatment.... Anyway, I think I made the table wrong, would anyone check what I did? I check it over and over again and keep changing settings but still didn't hit the core of the problem.

a=ctrl, b=ctrl, a&b are biological replicates; c=treatment, d=treatment, c&d are biological replicates



library(GenomicFeatures)
library(GenomicRanges)
library(Rsamtools)
library(edgeR)
library(org.Mm.eg.db)

txdb=makeTranscriptDbFromUCSC(genome="mm9",tablename="refGene")
ex_by_gene=exonsBy(txdb,"gene")

####read .bam file in object
readsa=readBamGappedAlignments("/home/vagrant/genedata/a_accepted_hits.bam")
readsb=readBamGappedAlignments("/home/vagrant/genedata/b_accepted_hits.bam")
readsc=readBamGappedAlignments("/home/vagrant/genedata/c_accepted_hits.bam")
readsd=readBamGappedAlignments("/home/vagrant/genedata/d_accepted_hits.bam")
#####count reads
countsa = countOverlaps(ex_by_gene,readsa)
countsb = countOverlaps(ex_by_gene,readsb)
countsc = countOverlaps(ex_by_gene,readsc)
countsd = countOverlaps(ex_by_gene,readsd)

#### making count table
countTableabcd = data.frame(conditiona=countsa,conditionb=countsb,conditionc=countsc,conditiond=countsd,stringsAsFactors=FALSE)
rownames(countTableabcd)=names(ex_by_gene)
write.table(countTableabcd,file="countTableabcd.txt",sep="\t")
#### eliminate empty rows
x <- rowSums(countTableabcd==0)!=ncol(countTableabcd)
newCountTableabcd <- countTableabcd[x,]
write.table(newCountTableabcd,file="newcountTableabcd.txt",sep="\t")

[dpryan]

Just use featureCounts or htseq-counts to count the reads. I've never been a fan of the R methods, particularly since they usually involve reading the whole file into memory.

[bruce01]

You can use featureCounts in R, if the OP wants to do it all in that wonderful enviRonment=)

[lucer105]

Sounds like it is easy to fix.

But I just wondering, why did my method failed? At least for counting that single gene....