Start with:
v <- Views(coverage(gr2, weight="score"), as(gr1, "RangesList"))
That object is a list by chromosome, where each element is a set of "views" on the coverage vector. Perhaps it does what you need.
If you really need a list that aligns with your original GRanges, we need to do something different and it gets more complicated. First, we compute the coverage from gr2
, as before. It is a list of vectors, one per chromosome. We then subset the vectors for the values inside the ranges of gr1
.
cvg <- coverage(gr2, weight="score")
cvg.gr1 <- cvg[as(gr1, "RangesList")]
If your gr1 is sorted by chromosome, just unlist the list into a single vector (Rle):
rle <- unlist(cvg.gr1)
If it is NOT sorted, do this instead:
rle <- unsplit(cvg.gr1, rep(seqnames(gr1), width(gr1)))
Now, we just have to form the list with one element per range:
rle.range <- relist(rle, ranges(gr1))
And finally reverse the elements that are on the negative strand:
rle.range[strand(gr1)=="-"] <- revElements(rle.range[strand(gr1)=="-"])
Sorry that it's an insane 5 lines of code to do something so simple. In theory it could be made as simple as cvg[gr1]
. Actually, it looks like that actually works, except for the strand-sensitive reversing, so we are here:
scores <- coverage(gr2, weight="score")[gr1]
scores[strand(gr1)=="-"] <- revElements(scores[strand(gr1)=="-"])
Hi,
It's a little bit confusing. The question in the title of your post looks very much like the question you already asked here how much positions overlapped for each pair of the pairs obtained by findoverlap. 8 days ago and for which you got an answer (that you accepted). But then the details you give make it sound a little bit different, maybe.
Could you please clarify with an example? Would help a lot if you could show 2 small GRanges objects
gr1
andgr2
and what how expect the score vector to be. Thanks!H.
output
I want vectors for each genomic range of gr1 of length equal to width of gr1. Required vector is such that it have value of score of gr2 at overlapping positions of gr1 else zero at non-overlaping positions. eg. out put
vector for first genomic range of gr1 : [0 .4 .4 .3 .3 .3] reverse this, as it is on negative strand of gr1
vector for second genomic range of gr1 : [0 0 .7 .7 0 0] leave this as, it is as it is on +ve strand of gr1
..
vector for fourth genomic range of gr1 : [0 .1 .1 .1 .1 .1] leave this as, it is as it is on +ve strand of gr1