text: user2884: 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
  
  
  user776: 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<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!
  
original_text: 2884: <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>
   776: <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>
  