Valarie, By calling: olap <- findOverlaps(gff0, reads) foo <- countSubjectHits(olap) `countSubjectHits(x): Counts the number of hits for each subject, returning an integer vector.` foo is a vector with a length equal to the number of reads mapped? so: length(foo[foo == 0]) would tell me how many reads did not map to my annotated features. can you recommend a method that would allow me to map to the regions outside of my current features? I have been trying to devise a way to make a new set of features that are the spaces between my current features, but I have been unable to implement it successfully. Maybe a more genomics oriented approach and look at the coverage over each chromosome to try and identify regions that have a higher 'background' that could be then looked at in a genome browser like IGV to see the reads that did not map to a transcript? Thanks again, Sam On Wed, Aug 28, 2013 at 5:38 PM, Valerie Obenchain <vobencha@fhcrc.org>wrote: > After a second read I think I may have missed what you're really after. > Specifically the question about 'where did the reads go'? > > A simple check of whether or not the mapped reads hit any part of the > annotation would be to overlap them with the 'gff0' GRanges. > > olap <- findOverlaps(gff0, reads) > > The result is a 'Hits' object that tells you what query overlapped what > subject. See ?Hits for details. > > ## number of hits for each row of gff0 > countQueryHits(olap) > ## number of hits for each read > countSubjectHits(olap) > > > Valerie > > > On 08/28/2013 03:15 PM, Valerie Obenchain wrote: > >> Hi Sam, >> >> ## counting reads >> >> If you haven't done this already I'd confirm counts on a single file. A >> bam can be read in with readGAlignments(). If you have paired-end data >> you can use ?readGAlignmentPairs or ?readGAlignmentsList. >> >> bf <- BamFile("pathToOneBam") >> reads <- readGAlignments(bf) ## devel >> reads <- readGappedAlignments(bf) ## release >> >> Take a look at the makeTranscriptDbFromGFF() function in >> GenomicFeatures. It creates a TranscriptDb object from gff or gtf. The >> advantage of working with a TxDb is that you get the accessors such as >> cdsBy(..., "gene"), etc. for free (you may know this, just an fyi). >> >> Get an idea of how many reads hit each list element of the GRangesList. >> This method does double count so it's possible for the sum of 'olap' to >> be greater than the length of 'reads'. >> >> olap <- countOverlaps(cds_gene, reads) >> >> If the numbers in 'olap' are very small (and this is unexpected) you may >> have a problem with how you've created the 'cds_gene' annotation. If the >> numbers make sense then try summarizeOverlaps. >> >> se <- summarizeOverlaps(cds_gene, bf, mode="IntersectionNotEmpty") >> >> If the number of hits in assays(se)$counts are greatly reduced you may >> have many reads overlapping each other or a single feature in >> 'cds_gene'. This can be checked by looking at the numbers in 'olap'. If >> you are using the devel version (GenomicRanges >= 1.13.9) the >> 'inter.feature' argument allows reads to be counted for each feature >> they overlap. >> >> >> ## mapped vs unmapped reads >> >> readGAlignments() only imports mapped reads; unmapped are dropped. If >> you are interested in the number of mapped vs unmapped and some details >> about them you can use scanBam(). If your file is large you can use the >> 'yieldSize' argument to BamFile. See ?BamFile and ?ScanBamParam for >> param details. >> >> All reads in the file, mapped and unmapped. >> scn <- scanBam(bf) >> names(scn[[1]]) >> >> Unmapped only. >> flag1 <- scanBamFlag(isUnmappedQuery=**TRUE) >> param1 <- ScanBamParam(what=scanBamWhat(**), flag=flag1) >> unmapped <- (bf, param1) >> >> Mapped only. >> flag2 <- scanBamFlag(isUnmappedQuery=**FALSE) >> param1 <- ScanBamParam(what=scanBamWhat(**), flag=flag2) >> mapped <- (bf, param2) >> >> >> Valerie >> >> On 08/28/2013 01:50 PM, Sam McInturf wrote: >> >>> 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 >>> >>> >> ______________________________**_________________ >> Bioconductor mailing list >> Bioconductor@r-project.org >> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https: stat.e="" thz.ch="" mailman="" listinfo="" bioconductor=""> >> Search the archives: >> http://news.gmane.org/gmane.**science.biology.informatics.**conduct or<http: news.gmane.org="" gmane.science.biology.informatics.conductor=""> >> > > -- Sam McInturf [[alternative HTML version deleted]]