Picking Random Genomic Positions

Question

[PoGibas]

I do have a set of TF binding coordinates and want to see if there is any significant overlap with an open chromatin annotation.

Example of TF coord:
chr1 19280 19298
chr1 245920 245938
chr2 97290 97308
chr9 752910 752938
...

Example of open chrom. coord. (UCSC track):
chr2 33031543 33032779
chr3 2304169 2304825
chr5 330899 330940
...

I have checked the intersection with the Bedtools (open chrom. coord vs TF coord. -/+ 100bp) and now I want to check the intersection between random genomic coordinates and open chrom.

The idea is to:

1. Pick random genomic position (from the same chromosome as TF coordinate);
2. -/+9bp (binding site size);
3. -/+ 100bp;
4. Run this simulation for 1000 times (TF x 1000);
5. Bedtools;

Any ideas how can I do this simulation to pick random genomic positions from the same chromosome? I know a little bit of bash and Perl, but won't be able to write the script by myself.
Is it possible to measure the length of every chromosome;
Pick TF chromosome and from it's length get a random number which would represent a genomic position?

Can someone help me with the simulation and the pipeline.

Answer 1

[Aaronquinlan]

My sense is that you are trying to do a Monte Carlo simulation and as such, you actually want to use the shuffle command. Shuffling will choose a new location for each of your original intervals while preserving its size. The loop below will print out how many intersections were observed for each of 1000 shuffles. A P-value can be derived by counting the number of times that the number of shuffled intersections exceeds your observed intersections. If 0, then you can say that P is less than 0.001.

#!/bin/bash
for i in {1..1000}
do
        OBS=`bedtools shuffle -i openchrom.bed -g hg19.chrom.sizes | \
            bedtools intersect -a tf.bed -b - | wc -l`
        echo $OBS 
done

If you really want a new set of randomly generated intervals, there is a random function in bedtools

bedtools random -n 1000 -l 118 -g hg19.chrom.sizes

If you are trying to measure the evidence for statistically significant relationship, the pybedtools package has nice tools for this. See the docs here as well as the intersection_matrix.py tool.

Lastly, there is a great new package called GenometriCorr in R, as well as an associated paper describing some interesting new metrics for measuring relationships between sets of genomic intervals.

[Shicheng Guo]

Some time we need keep similar GC contents or some other metric, any idea to update the script?

[Sebastian Kurscheid]

Very interesting paper - thanks for the link!

[Istvan Albert]

I love the new bedtools works with commands like samtools - I did not realize it up to know, the manuals still list the old scripts - I have been meaning to suggest it but did not want to look like I don't appreciate the existing tools already. Great link to the paper as well!

[Aaronquinlan]

@Istvan - yeah sorry, the docs are out of date. in the process of hiring someone to take over all of that and we should have a new docs website up by Fall.

[Biomonika (Noolean)]

Where can I please download hg19.chrom.sizes?

[Aaronquinlan]

If you have mysql installed, you can do: mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e \ "select chrom, size from hg19.chromInfo" > hg19.chrom.sizes

Otherwise, you can use the UCSC Table Browser. http://genome.ucsc.edu/cgi-bin/hgTables?command=start Group = "All tables" Table=chromInfo

Answer 2

[Sebastian Kurscheid]

As you are already using bedtools, you could perhaps work with the function "shufflebed"? It is well documented on the developer's website: http://code.google.com/p/bedtools/wiki/Usage#shuffleBed

[Biomonika (Noolean)]

I have just tryed it as well and it works just perfect.

Answer 3

[Istvan Albert]

A simple one liner in perl would look like this:

perl -e '$x=int(rand(1000)); printf( "chr1\t%s\t%s", $x, $x + 10)'

You just need to add this into a loop and tweak it to your needs.