Question: (Closed) Subscript error from TSS proximity function
0
12 months ago by
rbronste60
rbronste60 wrote:

Hi all,

Have a set of peaks and TSSs, trying to use the following function to retrieve the distance from peaks to nearest TSS and getting this error:

Wondering you anyone has an idea about why? Thanks!

# Get distance from peaks to the nearest TSS
tmp <- getNearestDist(all.peaks,tss,return.nearest=T)

Error in .Call2("map_positions", run_lens, pos, method, PACKAGE = "S4Vectors") :

subscript contains NAs

​function(gr1, gr2, return.nearest=F, r1.pos="mid", r2.pos="mid", strand=2)
{
r1 <- ranges(gr1)
r2 <- ranges(gr2)
map.id <- rep(NA, length(r1))
dist <- c()
for(chr in sort(unique(seqnames(gr1)))){
select1 <- seqnames(gr1)==chr
gr1.select <- gr1[select1,]
select2 <- seqnames(gr2)==chr
gr2.select <- gr2[select2,]
l <- nearest(ranges(gr1.select), ranges(gr2.select))
map.id[which(select1)] <- which(select2)[l]

start1=start(gr1.select)
end1 = end(gr1.select)
end5.1 <- start1
end3.1 <- end1
neg1 <- strand(gr1.select) == "-"
end5.1[which(neg1)] <- end1[which(neg1)]
end3.1[which(neg1)] <- start1[which(neg1)]
mid1 <- as.integer( (start1+end1)/2)

start2=start(gr2.select)
end2 = end(gr2.select)
end5.2 <- start2
end3.2 <- end2
neg2 <- strand(gr2.select) == "-"
end5.2[which(neg2)] <- end2[which(neg2)]
end3.2[which(neg2)] <- start2[which(neg2)]
mid2 <- as.integer( (start2+end2)/2)

pos1 <- switch(r1.pos,
mid=mid1,
start=start1,
end = end1,
end5 = end5.1,
end3 = end3.1)
pos2 <- switch(r2.pos,
mid=mid2,
start=start2,
end = end2,
end5 = end5.2,
end3 = end3.2)
tmp.dist <- as.integer(pos1 - pos2[l])
if(strand==2){
tmp.dist[which(neg2[l])] <- - tmp.dist[which(neg2[l])]
}
if(strand==1){
tmp.dist[which(neg1)] <- - tmp.dist[which(neg1)]
}
dist[which(select1)] <- tmp.dist
}
if(!return.nearest){
dist
}
else{
gr <- gr2[map.id,]
values(gr)$dist <- dist gr$map.id = map.id
gr
}
}​
diffbind s4vectors function tss • 268 views
modified 12 months ago by Michael Love25k • written 12 months ago by rbronste60

Hello rbronste!

We believe that this post does not fit the main topic of this site.

This support site is for Bioconductor packages, not for help debugging your own code.

For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.

If you disagree please tell us why in a reply below, we'll be happy to talk about it.

Cheers!

Is

PACKAGE = "S4Vectors"

not a Bioconductor package?

Yes. Is getNearestDist a function in S4Vectors?

You could try nearestTSS() in the edgeR package, although it works on genomic loci rather than genomic intervals.