question about summrizeOverlaps
1
0
Entering edit mode
Ou, Jianhong ★ 1.3k
@ou-jianhong-4539
Last seen 9 hours ago
United States

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    

 

summarizeoverlaps • 995 views
ADD COMMENT
2
Entering edit mode
@herve-pages-1542
Last seen 4 hours ago
Seattle, WA, United States

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.

ADD COMMENT
0
Entering edit mode

Good to know. I misunderstand the meaning of GRangesList in this function. Thank you.

ADD REPLY

Login before adding your answer.

Traffic: 636 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6