<p>You can perhaps collapse requirements for (1), (2), and (6), so long as bounds from (1) and (2) fall within the bounds of (6). </p>

<p>To help with this, you could use <code>fetchChromSizes</code> from the Kent UCSC utilities to quickly get chromosomes and chromosomal bounds for your genome of interest. </p>

<p>Given a chromosome and a <code>maximumEnd</code> that falls within its bounds, uniformly sample from <code>[0, maximumEnd-minimumStart-maxLength)</code>, adding <code>minimumStart</code> to determine start position and <code>maxLength</code> as length to get the stop position.</p>

<p>To obtain the desired number of samples in an efficient way, you could pipe these "candidate elements" into <code>bedmap --echo-map-size</code> operations on repeatmasked regions (4) and <code>bedops -n 1</code> on "regions-to-avoid" (5) operations to do <a href="https://en.wikipedia.org/wiki/Rejection_sampling" rel="nofollow noreferrer">rejection sampling</a>. </p>

<p>Perhaps something like:</p>

<pre><code>$ REPEATMASK_THRESHOLD=123
$ SAMPLES=1000
$ someScriptThatGeneratesCandidateIntervalsWithinBounds \
    | sort-bed - \
    | bedmap --echo-map-size --echo --delim '\t' - repeatmaskedRegions.bed \
    | awk -vthreshold=${REPEATMASK_THRESHOLD} '($1&lt;threshold)' \
    | cut -f2- \
    | bedops -n 1 - regionsToAvoid.bed \
    | head -${SAMPLES} \
    &gt; qualifyingCandidates.bed
</code></pre>

<p>This will be relatively efficient, because once the <code>head</code> process at the end writes the number of <code>SAMPLES</code> to <code>qualifyingCandidates.bed</code>, it will send a <code>SIGINT</code> signal up the pipeline that will terminate the upstream processes. So you're only doing sampling to the extent needed.</p>

<p>Hope this helps!</p>

CC BY-SA 4.0
2018-06-02T01:35:14.480Z
2018-07-02T04:56:52.460Z
2018-07-02T04:56:52.460Z
Alex Reynolds
<p>Problem is to generate a random BED interval given the following constraints:</p>

<ol>
<li>minimum <code>start</code></li>
<li>maximum <code>end</code></li>
<li>fixed length</li>
<li>maximum number of masked bases (similar to <code>-maxN</code> option in <code>faSplit</code>)</li>
<li>set of intervals to avoid overlap with</li>
<li>stay within chromosome size constraints</li>
</ol>

<p>I can think of computationally intensive ways to code it up from scratch, but am wondering if there is a more efficient approach before re-inventing the wheel.</p>

<pre><code>seq_records = {x.name: x for x in SeqIO.parse('path/to/genome.fa', 'fasta')}

def generate_random_interval(chrom, lower, upper, length, maxrep=1.0, avoid_intervals=None):
    while True:
        start = np.random.randint(lower, upper - length)
        end = start + length
        regenerate = False

        for interval in avoid_intervals:
            if (interval.start &lt; end) or (interval.end &gt; start):
                regenerate = True
        if 1. * seq_records[chrom].seq[lower, upper].count('N') / length &gt; maxrep:
            regenerate = True

        if not regenerate:
            break

    return chrom, start, end
</code></pre>

<python><bed><bedtools>
Generating random BED intervals given constraints
# Generating random BED intervals given constraints

(Asked by: username_0 on 2018-06-01)

Problem is to generate a random BED interval given the following constraints:

minimum start
maximum end
fixed length
maximum number of masked bases (similar to -maxN option in faSplit)
set of intervals to avoid overlap with
stay within chromosome size constraints

I can think of computationally intensive ways to code it up from scratch, but am wondering if there is a more efficient approach before re-inventing the wheel.
seq_records = {x.name: x for x in SeqIO.parse('path/to/genome.fa', 'fasta')}

def generate_random_interval(chrom, lower, upper, length, maxrep=1.0, avoid_intervals=None):
    while True:
        start = np.random.randint(lower, upper - length)
        end = start + length
        regenerate = False

        for interval in avoid_intervals:
            if (interval.start < end) or (interval.end > start):
                regenerate = True
        if 1. * seq_records[chrom].seq[lower, upper].count('N') / length > maxrep:
            regenerate = True

        if not regenerate:
            break

    return chrom, start, end

**Tags:** python, bed, bedtools

---

## 1 answers

### Answer by username_1 on 2018-06-02. Score: 2

You can perhaps collapse requirements for (1), (2), and (6), so long as bounds from (1) and (2) fall within the bounds of (6). 
To help with this, you could use fetchChromSizes from the Kent UCSC utilities to quickly get chromosomes and chromosomal bounds for your genome of interest. 
Given a chromosome and a maximumEnd that falls within its bounds, uniformly sample from [0, maximumEnd-minimumStart-maxLength), adding minimumStart to determine start position and maxLength as length to get the stop position.
To obtain the desired number of samples in an efficient way, you could pipe these "candidate elements" into bedmap --echo-map-size operations on repeatmasked regions (4) and bedops -n 1 on "regions-to-avoid" (5) operations to do rejection sampling. 
Perhaps something like:
$ REPEATMASK_THRESHOLD=123
$ SAMPLES=1000
$ someScriptThatGeneratesCandidateIntervalsWithinBounds \
    | sort-bed - \
    | bedmap --echo-map-size --echo --delim '\t' - repeatmaskedRegions.bed \
    | awk -vthreshold=${REPEATMASK_THRESHOLD} '($1 qualifyingCandidates.bed

This will be relatively efficient, because once the head process at the end writes the number of SAMPLES to qualifyingCandidates.bed, it will send a SIGINT signal up the pipeline that will terminate the upstream processes. So you're only doing sampling to the extent needed.
Hope this helps!


---

