Generate random genomic positions with R

Question

[Nicolas Rosewick]

Hi,

Is there a package/function to generate random genomic position in R ? If not it's not complicated to write something using sample() but I do not want to re-re-re-reinvent the wheel ..

Edit : I wrote that. Some benchmarking with other given solutions in this thread seems to be required !

generateRandomPos <- function(n,chr,chr.sizes,width=1){
    random_chr <- sample(x=chr,size=n,prob=chr.sizes,replace=T)
    random_pos <- sapply(random_chr,function(chrTmp){sample(chr.sizes[chr==chrTmp],1)}) 
    res <- GRanges(random_chr,IRanges(random_pos,random_pos+width))
    return(res)
}

Thanks

[GenoMax]

(Davetang is the author from gist.github.com/davetang/6548010):

	#how many regions to sample?
	number <- 50000
	#how many bp to add to start
	span <- 4
	#some necessary packages
	#install if necessary
	source("http://bioconductor.org/biocLite.R")
	biocLite("BSgenome")
	#load library
	library(BSgenome)
	#find available genomes
	available.genomes()
	#I'm interested in hg19
	biocLite("BSgenome.Hsapiens.UCSC.hg19")
	library(BSgenome.Hsapiens.UCSC.hg19)
	#chromosomes of interest
	my_chr <- c(1:22,'X','Y')
	my_chr <- gsub(pattern="^", replacement='chr', my_chr)
	#initialise list to store chromosome sizes
	my_chr_size <- list()
	for (i in my_chr){
	my_chr_size[[i]] <- length(BSgenome.Hsapiens.UCSC.hg19[[i]])
	}
	#checkout my_chr_size
	head(my_chr_size,2)
	#$chr1
	#[1] 249250621
	#
	#$chr2
	#[1] 243199373
	#initialise some vectors for storing random coordinates
	my_random_start <- vector()
	my_random_end <- vector()
	my_random_chr <- vector()
	my_random_strand <- vector()
	set.seed(12345)
	#loop through number of regions
	for(i in 1:number){
	my_random_chr[i] <- sample(x=my_chr,size=1)
	my_random_strand[i] <- sample(x=c('-','+'),size=1)
	my_max <- my_chr_size[[my_random_chr[i]]]-span
	my_random_start[i] <- runif(n=1, min=1, max=my_max)
	my_random_end[i] <- my_random_start[i] + span
	}
	my_random_loci <- data.frame(chr=my_random_chr,
	start=my_random_start,
	end=my_random_end,
	strand=my_random_strand)
	head(my_random_loci)
	# chr start end strand
	#1 chr18 59415403 59415407 +
	#2 chr22 8535632 8535636 -
	#3 chr8 106509865 106509869 +
	#4 chrY 9046958 9046962 -
	#5 chr18 30544079 30544083 -
	#6 chr12 53873398 53873402 -

view raw random_region.R hosted with ❤ by GitHub

[Nicolas Rosewick]

In fact it's not really random as my_random_chr[i] <- sample(x=my_chr,size=1) will pick a chromosome unregarding of its size. Anyway thank you I'll will rewrite the code for my usage and post it here ;)

[Giovanni M Dall'Olio]

Nice but I don't like the for loop. It would be easy and much more elegant to convert it to a function and apply.

[Nicolas Rosewick]

I wrote that : chr= chromosome names and chr.sizes = their associated size

generateRandomPos <- function(n,chr,chr.sizes,width=1){
    random_chr <- sample(x=chr,size=n,prob=chr.sizes,replace=T)
    random_pos <- sapply(random_chr,function(chrTmp){sample(chr.sizes[chr==chrTmp],1)}) 
    res <- GRanges(random_chr,IRanges(random_pos,random_pos+width))
    return(res)
}

[Macherki M E]

You need a random position .may be with sample(1:length sequence,1,true) will gives one sample

[ivivek_ngs]

what about bedtools [shuffle][1]. You can take a look at this link as well. Ideally they allow to pick up regions from the original master file. I guess you can also use this code to generate one.

Answer 1

[Giovanni M Dall'Olio]

I don't think there is a function for it already, but I suggest you to get chromosome lengths from the TxDb.Hsapiens.UCSC.hg19.knownGene object from Homo.sapiens.

> library(Homo.sapiens)
> seqlengths(TxDb.Hsapiens.UCSC.hg19.knownGene)[1:24]
     chr1      chr2      chr3      chr4      chr5      chr6      chr7      chr8
249250621 243199373 198022430 191154276 180915260 171115067 159138663 146364022
     chr9     chr10     chr11     chr12     chr13     chr14     chr15     chr16
141213431 135534747 135006516 133851895 115169878 107349540 102531392  90354753
    chr17     chr18     chr19     chr20     chr21     chr22      chrX      chrY
 81195210  78077248  59128983  63025520  48129895  51304566 155270560  59373566

> rr  = data.frame(chr=sample(seqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene)[1:24], 40, replace=F))
> rr$start = apply(rr, 1, function(x) { 
         round(runif(1, 0, seqlengths(TxDb.Hsapiens.UCSC.hg19.knownGene)[x][[1]]), 0) }
         )
> rr$end = rr$start + runif(1, 1, 1000)  # random intervals of 1-1000 bases
> rr$strand = sample(c("+", "-"))

> rr.gr = makeGRangesFromDataFrame(rr)
> rr.gr
GRanges object with 40 ranges and 0 metadata columns:
       seqnames                 ranges strand
          <Rle>              <IRanges>  <Rle>
   [1]     chr1 [225055397, 225056035]      +
   [2]    chr12 [125243701, 125244339]      -
   [3]    chr22 [ 42744305,  42744943]      +
   [4]     chr1 [ 79241981,  79242619]      -
   [5]     chr7 [151649796, 151650434]      +
   ...      ...                    ...    ...
  [36]    chr21   [33907192, 33907830]      -
  [37]    chr18   [73931535, 73932173]      +
  [38]     chrY   [14300243, 14300881]      -
  [39]    chr17   [26049438, 26050076]      +
  [40]    chr22   [26631730, 26632368]      -
  -------

[GabrielMontenegro]

Is there a way to make random genomic ranges, but avoiding overlap between the segments? I want to randomly place non-overlapping segments (with a given length[width]) across the genome. Thanks.