Generating random BED intervals given constraints

Question

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

Answer 1

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!