Search
Question: (Closed) Subscript error from TSS proximity function
0
5 days 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
}
}​
modified 4 days ago by Michael Love19k • written 5 days 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!

ADD REPLYlink written 4 days ago by James W. MacDonald48k

Is

PACKAGE = "S4Vectors"

not a Bioconductor package?

ADD REPLYlink modified 4 days ago • written 4 days ago by rbronste60

Yes. Is getNearestDist a function in S4Vectors?

ADD REPLYlink written 4 days ago by James W. MacDonald48k

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