Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-06-02T01:35:14.480Z
Answers.LastActivityDate: 2018-07-02T04:56:52.460Z
Answers.LastEditDate: 2018-07-02T04:56:52.460Z
Answers.UserDisplayName: Alex Reynolds
Body: <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>
  
Tags: <python><bed><bedtools>
Title: Generating random BED intervals given constraints
ThreadText: # 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:
  
  1. minimum `start`
  2. maximum `end`
  3. fixed length
  4. maximum number of masked bases (similar to `-maxN` option in `faSplit`)
  5. set of intervals to avoid overlap with
  6. 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](https://en.wikipedia.org/wiki/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<threshold)' \
      | cut -f2- \
      | bedops -n 1 - regionsToAvoid.bed \
      | head -${SAMPLES} \
      > 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!
  
  
  ---
  
  