Bioconductors,
    I am working on an Arabidopsis RNA sequencing experiment and have run
into some trouble.  I found that I was able to map ~94% of my ~35 million
reads / library to the Arabidopsis TAIR10 genome using tophat2 (100 bp
single end reads).  After I counted my mapped reads I found that only about
47% of my mapped reads were counted. I found out there is a strong
possibility that there may be genomic DNA contamination, and I would like
to verify this hypothesis.
I used GenomicRanges and rtracklayer to do my read counting:

library(GenomicRanges)
library(rtracklayer)
library(Rsamtools)
setwd("/myDir")
myBamFiles <- list.files(pattern = ".bam$")
myBamIndex <- list.files(pattern = ".bai$")

bfl <- BamFileList(myBamFiles, index = myBamIndex)
gff0 <- import("/myOtherDir/Arabidopsis_thaliana.TAIR10.17.gtf",
asRangedData = FALSE)
cds <- subset(gff0, type == "CDS")
cds_gene <- reduce(split(cds, cds$gene_id))
olap <- summarizeOverlaps(cds_gene, bfl, mode = "IntersectionNotEmpty")

My intention with subsetting for the CDS and then splitting on gene_id is
to create a GRangesList where each element is a unique AGI number and has
all of the features.
To extract the number of reads counted I use:

colSums(assays(olap)$counts)

The output of colSums(...) is how I have found that less than half my reads
mapped to transcripts.  The first question becomes "Is this how to count
the number of reads mapped?"  The second question becomes "how do I figure
out where the reads went?"

My initial inclination would be to make a GRangesList that contained
intervals between genes and then do the same thing as before and see if I
can find my other reads.  Something in the vein of:

maxLimits <- t(as.dataframe(lapply(cds_gene, function(x)
c(min(start(ranges(x))),
   max(end(  ranges(x)))
   )
)
))
which is just an elaborate lapply call to get the smallest start coordinate
and the largest end coordinate for each unique gene_id and then transform
it to a pretty data frame.  From there I could find a way to get the
'in-betweens' of these features, but I am not sure that is the most prudent
thing to do.  I looked at some of the ShortRead and GenomicRanges functions
for a guiding hand and found the disjointBins() function may be useful in
establishing features to map over, but am at a bit of a loss as I don't
have the expertise to figure this out straight away.
  Any suggestions or resources to find mapped but uncounted reads would be
helpful!


Thanks
-- 
Sam McInturf

	[[alternative HTML version deleted]]

