Hello everyone, I am having trouble counting RNA seq reads with GenomicRanges and GenomicFeatures. Essentially the problem is I have many reads (10 million mapped) and I only can count 200K in features and 300K between features. So obviously I am doing something wrong. I am working with Arabidopsis thaliana, so I mapped to the TAIR10.17 release with tophat2. The gff file I use is TAIR10.17.gtf, I also tried using TAIR10.16.gtf. I also made a GRangesList with the 10.17 with no effect. This is what I have done while using 1 library #outside of R stuff tophat2 (no -G ) blah blah samtools sort accepted_hits.bam foo samtools index foo #In R #load bams bf <- "myBam.bam" bi <- "myBam.bai" bfl <- BamFileList(bf, bi) #set up TxDb gff <- makeTranscriptDbFromGFF(gff) tx <- transcriptsBy(gff) #summarize olap <- summarizeOverlaps(tx, bfl) #how many reads mapped? colSums(assays(olap)$counts) #only 222,835! are the reads inbetween features? abf <- BamFile("path to a single bam file") reads <- readGappedAlignments(bf) tx2 <- transcripts(gff) strand(tx2) <- "*" txs <- reduce(tx2) nontx <- gaps(txs) o <- summarizeOverlaps(nontx, reads, ignore.strand = TRUE) sum(o) #365401 So what happened? If the reads are not in the features, and not inbetween the features, where did they go? Anything helps, and thanks in advance Sam > sessionInfo() R version 3.0.2 (2013-09-25) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] AnnotationDbi_1.24.0 Biobase_2.22.0 BiocGenerics_0.8.0 [4] biomaRt_2.18.0 Biostrings_2.30.1 bitops_1.0-6 [7] BSgenome_1.30.0 DBI_0.2-7 GenomicFeatures_1.14.2 [10] GenomicRanges_1.14.4 IRanges_1.20.6 parallel_3.0.2 [13] RCurl_1.95-4.1 Rsamtools_1.14.2 RSQLite_0.11.4 [16] rtracklayer_1.22.0 stats4_3.0.2 tools_3.0.2 [19] XML_3.98-1.1 XVector_0.2.0 zlibbioc_1.8.0 -- Sam McInturf [[alternative HTML version deleted]]