Search
Question: (Closed) Subscript error from TSS proximity function
0
gravatar for rbronste
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
  }
}​
ADD COMMENTlink 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.

ADD REPLYlink written 3 days ago by Gordon Smyth35k
Please log in to add an answer.
The thread is closed. No new answers may be added.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 444 users visited in the last hour