Hi,

I have question about summarizeOverlaps when set reads as a GRangesList object. In the Usage section, there is no method for signature 'GRanges,GRangesList'. However, in the arguments section, it is said that GRangesList could be the input of reads. Then I tried codes like this,

library(GenomicAlignments)
tg <- GRanges("chr1", IRanges(c(1, 101, 201), width=100))
reads <- GRangesList(A=GRanges("chr1", IRanges(1:100, width=10)))
tc <- summarizeOverlaps(features=tg, reads=reads)
assays(tc)[[1]]
     reads
[1,]     0
[2,]     0
[3,]     0

tc2 <- summarizeOverlaps(features=tg, reads=reads$A)
assays(tc2)[[1]]

     reads
[1,]    91
[2,]     0
[3,]     0

I hope summarizeOverlaps could accept GRangesList as inputs for reads. Is there anything I missed? Thank you.

Here is my sessionInfo

R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin16.1.0 (64-bit)
Running under: macOS Sierra 10.12.3

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] GenomicAlignments_1.10.0   Rsamtools_1.26.1           Biostrings_2.42.1          XVector_0.14.0             SummarizedExperiment_1.4.0
 [6] Biobase_2.34.0             GenomicRanges_1.26.2       GenomeInfoDb_1.10.3        IRanges_2.8.1              S4Vectors_0.12.1          
[11] BiocGenerics_0.20.0       

loaded via a namespace (and not attached):
[1] lattice_0.20-34    bitops_1.0-6       grid_3.3.2         zlibbioc_1.20.0    Matrix_1.2-8       BiocParallel_1.8.1 tools_3.3.2       
[8] RCurl_1.95-4.8    

Hi,

Your own examples show that it does accept reads represented as a GRangesList object. So maybe you're wondering why you get a different result when you replace reads with reads$A?

The GRangesList representation is interpreted as "one read per list element". For example, if an aligned read has a skipped region in it (i.e. an N in the cigar), then after coercing the GAlignments object to GRangesList the read is represented with a list element that is a GRanges object with 2 ranges:

> reads2 <- GAlignments("chr1", 10L, "25M30N25M", strand("+"))
> reads2
GAlignments object with 1 alignment and 0 metadata columns:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chr1      +   25M30N25M        50        10        89        80
          njunc
      <integer>
  [1]         1
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

> as(reads2, "GRangesList")
GRangesList object of length 1:
[[1]] 
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1  [10, 34]      +
  [2]     chr1  [65, 89]      +

-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths

Note that the list element representing a read will always be a disjoint GRanges object (i.e. a GRanges object where the ranges don't overlap), so your artificial read object is not realistic.  Anyway, the "one read per list element" interpretation means that your read object is seen as representing a single read that spans the region chr1:1-109 (note that you can do reduce(reads) to compute this region). So the result you get is what you would get if you were doing:

> assay(summarizeOverlaps(features=tg, reads=GRanges("chr1:1-109")))
     reads
[1,]     0
[2,]     0
[3,]     0

See the mode argument for why you get zero counts.

Hope this helps,

H.