Can I get read counts by doing this???

[wetSEQer]

I want to quantify the reads I have from several .bam files from Tophat.

I used two different annotations(one at a time):
1: CuffMerge generated GTF file;
2: UCSC genome annotation
I made a transcript data base by using Bioconductor package "GenomicFeature" function "makeTranscitpDbFromGFF" or "makeTranscitpDbFromUCSC"

Then I make a object that store the exon or transcripts information, by using Rsamtools package function "exonsBy" or "transcriptsBy"

Then I do the overlap between my .bam file and the previous object ("countOverlaps")

Anyone see any problems by doing what I did?

I took this as a reference:

http://biowhat.ucsd.edu/homer/basicTutorial/rnaseqR.html

Please let me know what could go wrong with this~

Thanks a lot!

[dpryan]

Be careful using countOverlaps:

library(IRanges)
subject <- IRanges(c(1, 4, 9), c(5, 7, 10))
query <- IRanges(c(2, 2, 10), c(2, 3, 12))
#If this were RNAseq data, the output should be 0 0 1
countOverlaps(query, subject)

You'll see that the counts would be wrong for RNAseq data. You might just use featureCounts or htseq-counts.