I want to compute the nearest neighbour distances using a for loop. I can do this with a random subset by doing this:

df <-  read.table("tmp.txt", header=F)
colnames(df) <- c("gene", "seqnames", "start", "stop", "LOCID")
gr <- makeGRangesFromDataFrame(df, keep.extra.columns=TRUE)

# genes of interest in a list
gene.list <- read.table("careplastic_diffs.txt")

# filters the GenomicRanges object with genes of interest
sub.gr <- gr[(elementMetadata(gr)[, "gene"] %in% gene.list$V1)]

# compute distances
dists <- distanceToNearest(x = sub.gr)

mcols(dists)[,1] # contains a list of distances

I now want to essentially do this 5 times, using the same logic as above:

## create a rand dataframe, each column is a random list of 50 genes
rand.df <- data.frame(replicate(5, sample(gene.list$V1, 50)))

# loop through columns and gets distances and prints to a file
for(i in 1:ncol(rand.df)){
  rand.gr[i] <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists[i] <- distanceToNearest(x = rand.gr[i])
  output <- print(mcols(dists[i][,1]))

But I can't seem to make this work.

Error: subscript contains out-of-bounds indices In addition: Warning message: In recycleSingleBracketReplacementValue(value, x, nsbs) :
number of values supplied is not a sub-multiple of the number of values to be replaced

It's probably something really silly but any help would be appreciated ...


2 Answers 2


Is it because rand.gr has not been initialised, so when you ask to set rand.gr[i] it can't because there isn't one yet?

If you set i to a value (e.g. 1), and run the lines inside your loop, which of the lines produces the error?

If you don't need to retain all the random genomic ranges at the end, you could just leave out some of the [i]s, something like this:

for(i in 1:ncol(rand.df)){
  rand.gr <- gr[(elementMetadata(gr)[,"gene"] %in% rand.df[,i])]
  dists <- distanceToNearest(x = rand.gr)
  output <- print(mcols(dists[,1]))
  • $\begingroup$ Isn't output being overwritten with each iteration? $\endgroup$
    – haci
    Nov 8, 2019 at 12:38
  • $\begingroup$ Thanks, I do need a file with essentially columns of distance data (each column corresponding to a column in rand.df) $\endgroup$
    – R-MASHup
    Nov 8, 2019 at 13:17
  • $\begingroup$ this seemed to create another problem: dropbox.com/s/jl03i9nm87vwbxz/… $\endgroup$
    – R-MASHup
    Nov 8, 2019 at 13:23
  • $\begingroup$ haci is correct, ouput is overwritten on each time round the loop. $\endgroup$
    – Jay Moore
    Nov 12, 2019 at 10:38

You could do this much more easily and much faster on the command line, starting with an unsorted BED file of genes called genes.unsorted.bed:

$ sort-bed genes.unsorted.bed | shuf - | head -5 - | sort-bed - | closest-features --closest --dist --no-overlaps - > answer.bed

Replace the 5 in head -5 with however many random genes you want. Remove --no-overlaps if you want to allow consideration of overlapping genes (which returns a distance of zero).

Output columns are delimited with |. This can be replaced with another delimiter, but the output can be brought into R easily with read.table or similar.

Ref: https://bedops.readthedocs.io/en/latest/content/reference/set-operations/closest-features.html

  • $\begingroup$ i thought closest-features requires a reference and an input file? $\endgroup$
    – R-MASHup
    Nov 15, 2019 at 23:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.