Thanks Varlerie for looking into these suggestions and improvements. I
really appreciate the time and effort you keep investing into my
unsolicited questions. - Obviously, read counting has become a major
activity for many of us and your software will definitley put to great
use. In our projects, flexiblity and accuracy of the read counter is
most important to us, which summarizeOverlaps and associates address
very well. Since research facilities like ours have to perform those
analyses on very large disk storage devices with several hundred TBs
of
space that are shared among many NGS users, disk I/O and network
speeds
on compute clusters have become a major bottleneck. Thus, we are also
very interested in solutions that collect counts for many feature
types
and counting modes in a single pass-through of these big files, as
opposed to accessing them over and over again when it is not really
necessary. Usually, I also like to delete the bam files early on in
most
of our RNA-Seq workflows to save storage space. However, our
collaborators often keep coming up with all kinds of creative and
exciting new ideas what features/modes to count on where the best
analysis strategy is often to generate by default counts for many
feature types right away (even if I don't intend and recommend or to
use
them) rather than only one prescibed solution that usually ends up not
being enough. I guess RNA-Seq is still very much in an engineering
stage
where most of us depend on software design that provides a very high
level of flexiblitly.
Again, thanks a lot.
Thomas
On Thu, May 23, 2013 at 09:06:10PM +0000, Valerie Obenchain wrote:
> Hi Thomas,
>
> On 05/20/2013 10:12 PM, Thomas Girke wrote:
> > I apologize up front to Valarie for posting additional questions
related
> > to her excellent read counter. This time my questions are related
to the
> > economics of accessing big files as few times as possible to
improve
> > workflow performance in large RNA-Seq projects and to subsequently
allow
> > deleting all input bam files since BIG DATA is killing us :).
> >
> > (1) What is currently the intended approach to obtain counts for
> > multiple feature types in a single pass-through of one or many bam
> > files? For instance, often we want to count the reads overlapping
with
> > many different features types (not just one), such as exonBygenes,
> > cdsBygenes, UTRs, introns and intergenics. Computing all those
feature
> > counts sequentially in a loop is easy, but can be time consuming
and
> > heavy on disk I/O (the most expensive resource in NGS analysis)
and
> > internal networks when running many counting processes in
parallel. As a
> > partial solution to this, I often organize all feature types in
one big
> > GRanges/GRangesList container, perform the counting and then split
the
> > resulting count tables accordingly. However, is this really the
best way
> > of doing this? A major limitation of this approach is that it
does not
> > support feature overlap aware counting modes.
>
> I would not recommend combining different annotations into a single
> GRanges/GRangesList. To summarizeOverlaps(), each row of a GRanges
or
> each list element of a GRangesList represents a feature. If a read
hits
> more than one feature we have a double hit situation that must be
> resolved. If you're interested in counting both exonsByGene and
> transcriptsByGene you probably want these to be independent counts.
By
> combining them into one object they can be confounding and create
> multiple hits situations where there really aren't any.
>
> Toy example, but let's say 'gr1' comes from exonsBy() and 'gr2' from
> 3UTRsByTranscript().
>
> gr1 <- GRanges(c("chr1", "chr1"), IRanges(c(1, 10), width=10))
> gr2 <- GRanges("chr1", IRanges(25, width=10))
> grl <- GRangesList(gr1, gr2)
> reads <- GAlignments("chr1", pos=18L, cigar="10M",
strand=strand("*"))
>
>
> The read hits gr1 and gr2. If counted separately,
>
> > assays(summarizeOverlaps(grl[1], reads, Union))$counts
> reads
> [1,] 1
> > assays(summarizeOverlaps(grl[2], reads, Union))$counts
> reads
> [1,] 1
>
> Now count as GRangesList,
>
> > assays(summarizeOverlaps(grl, reads, Union))$counts
> reads
> [1,] 0
> [2,] 0
>
>
> Using 'inter.feature=FALSE' will regain the hits with modes Union
and
> IntersectionStrict but not necessarily IntersectionNotEmpty because
> shared regions of the annotation are removed before counting.
>
> > assays(summarizeOverlaps(grl, reads, Union,
inter.feature=FALSE))$counts
> reads
> [1,] 1
> [2,] 1
>
> >
> > (2) Similarly, what if I want to generate counts for multiple
counting
> > modes in a single call to summarizeOverlaps, such as different
overlap
> > modes, or sense and antisense counts which we often want to report
when
> > working with strand specific RNA-Seq data? Are there plans to
support
> > this? BTW: currently, I generate antisense read counts by
inversing the
> > strand information of the features which works fine just not for
both
> > ways in one step.
>
> We don't have a built-in method to count Bams over multiple
annotations
> and multiple count modes. There are several approaches you could
take.
> After some discussion this was one we came up with.
>
> ## Count a single file.
> count1bam <-
> function(file, features, mode, ...,
> yieldSize=1000000L, reader=readGAlignments,
> param=ScanBamParam())
> {
> ## initialize
> bf <- open(BamFile(file, yieldSize=yieldSize))
> counts <- lapply(sapply(features, length), integer)
>
> ## iterate
> while(length(reads <- reader(bf, param=param)))
> for (i in seq_along(features))
> counts[[i]] <-
> counts[[i]] + mode(features[[i]], reads, ...)
> close(bf)
> counts
> }
>
> ## Count over files in parallel.
> countbams <-
> function(features, files, ...)
> {
> counts0 <- mclapply(files, count1bam, features=features,
...)
> counts <- lapply(names(features), function(i, counts) {
> sapply(counts0, "[[", i)
> }, counts0)
> }
>
> ## In practice.
> library(parallel); options(mc.cores=detectCores())
> library(GenomicRanges)
> library(Rsamtools)
> features <-
> list(coarse=GRanges("seq1", successiveIRanges(rep(100,
15))),
> fine=GRanges("seq2", successiveIRanges(rep(20, 45))))
> fl <- system.file("extdata", "ex1.bam", package="Rsamtools")
> files <- list(fl, fl) ## usually different files
>
> ## As 'mode' you can directly call Union, IntersectionStrict,
> ## and IntersectionNotEmpty.
> result <- countbams(features, files, mode=Union,
> ignore.strand=TRUE, inter.feature=FALSE)
>
> ## 'reader' specifies how the data are read in. For paired-end
> ## use readGAlignmentsList (qname sorted or
readGAlignmentPairs.
> result <- countbams(features, files, mode=Union,
> reader=readGAlignmentPairs)
>
> This code handles the case of multiple Bam with multiple
annotations. To
> handle combinations for multiple arguments ('mode', 'ignore.strand',
> 'inter.feature', etc.) you could add a Map() in the 'while' loop in
> count1bam. Instead of passing a single value you would pass lists
the
> same length as the list of annotations. For example,
>
> modeList <- c(Union, IntersectionStrict)
> ignoreStrandList <- c(TRUE, TRUE)
> interFeatureList <- c(TRUE, FALSE)
>
> >
> > (3) Is it currently possible to return the total numbers of
aligned
> > reads in bam files with summarizeOverlaps and store them in the
> > resulting counting object? There are many places from where the
numbers
> > of aligned reads can be obtained, but to me the most obvious step
to
> > generate and store them would be right here. They are useful
parameters
> > for alignment QC'ing and reporting proper RPKM/FPKM values to
biologists.
> >
>
> I've added the output of countBam() to the 'colData' slot of the
> SummarizedExperiment. The 'mapped' counts are computed with
countBam()
> and a param where 'isUnmappedQuery=FALSE'. This isn't a great
example
> but does demonstrate the new output.
>
> library(Rsamtools)
> library(pasillaBamSubset)
> library("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
> exbygene <- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
> se <- summarizeOverlaps(exbygene, BamFile(untreated1_chr4()))
>
> > colData(se)
> DataFrame with 1 row and 3 columns
> records nucleotides mapped
> <integer> <numeric> <integer>
> untreated1_chr4.bam 204355 15326625 204355
>
>
>
> Code updates are in GenomicRanges 1.13.15 and Rsamtools 1.13.16.
>
> Valerie
>
> > Thomas
> >
> > _______________________________________________
> > Bioconductor mailing list
> > Bioconductor at r-project.org
> >
https://stat.ethz.ch/mailman/listinfo/bioconductor
> > Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
> >
