text: user146: I'm currently adding a few SNPs randomly into a FASTA within python using BioPython. In the following example from BioPython, I add an SNV at location "5"
  http://biopython.org/DIST/docs/api/Bio.Seq.MutableSeq-class.html
  from Bio.Seq import MutableSeq
  from Bio.Alphabet import generic_dna
  my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
  my_seq[5] = "A"
  print(my_seq)
  ## outputs ACTCGACGTCG 
  
  In my case, I'm using the reference hg38 as the input FASTA to manipulate. 
  After randomly inserting these SNPs, I would like to find a region whereby there are no SNPs within (let's say) +/- 50Kb, or the next nearest thing. 
  Is it possible to align the FASTA with variants against the reference within python? Otherwise, is there a quick way to check the variant FASTA against the reference FASTA, finding a region with no variants? What would be an efficient way to execute this task? 
  
  user776: If you store the positions of your random SNPs in a sort-bed sorted BED file, you can filter with BEDOPS bedmap:
  $ bedmap --count --echo --range 50000 --delim '\t' SNPs.bed | awk '$1==1' | cut -f2- > SNPs.filtered.bed
  
  The file SNPs.filtered.bed will only contain those elements from SNPs.bed which do not overlap within 50kb.
  Or you can use BEDOPS bedops and closest-features to do the same thing, but the logic isn't nearly as pretty:
  $ bedops --everything --range 0:50000 test.bed | closest-features --dist test.bed - | awk -v FS='|' '{ if ($NR > 50000) { print $1; } }' > SNPs.filtered.bed
  
  The result is the same as with the bedmap pipeline. IMO, the bedmap approach is cleaner.
  Once you have SNPs.filtered.bed, you could convert it to FASTA via samtools faidx and bed2faidx.pl or similar. Then you have FASTA upon which you can apply your random perturbations.
  
  user77: As others have said, doing whole-genome (or even whole-chromosome) alignments is the wrong solution. Simply track where you're creating SNVs. If you wrote your SNV locations to a sorted BED file you could use something like the following:
  #!/usr/bin/env python
  minWidth = 100000  # Only report regions of >=100kb
  last = [None, None]
  for line in open("some bed file"):
      cols = line.strip().split()
      if cols[0] != last[0]:
          if int(cols[2]) >= minWidth:
              print("{}:1-{}".format(cols[0], cols[2]))
      else:
          if int(cols[2]) - last[1] >= minWidth
              print("{}:{}-{}".format(cols[0], last[1], cols[2]))
      last = [cols[0], int(cols[2])]
  
  This will miss the last part of each chromosome/contig, since it doesn't know how long they are, but you should be able to figure out how to handle that. The memory requirements of this are trivially small.
  I suspect there's also a bedtools or bedops command that will directly give you the distance to the next downstream region. One could alternatively parse that (it's all the above code is doing, but using bedtools or bedops would produce more validated results than a quickly thrown together script).
  
original_text: 146: <p>I'm currently adding a few SNPs randomly into a FASTA within python using BioPython. In the following example from BioPython, I add an SNV at location "5"</p>
  
  <p><a href="http://biopython.org/DIST/docs/api/Bio.Seq.MutableSeq-class.html" rel="noreferrer">http://biopython.org/DIST/docs/api/Bio.Seq.MutableSeq-class.html</a></p>
  
  <pre><code>from Bio.Seq import MutableSeq
  from Bio.Alphabet import generic_dna
  my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
  my_seq[5] = "A"
  print(my_seq)
  ## outputs ACTCGACGTCG 
  </code></pre>
  
  <p>In my case, I'm using the reference hg38 as the input FASTA to manipulate. </p>
  
  <p>After randomly inserting these SNPs, I would like to find a region whereby there are no SNPs within (let's say) +/- 50Kb, or the next nearest thing. </p>
  
  <p>Is it possible to align the FASTA with variants against the reference within python? Otherwise, is there a quick way to check the variant FASTA against the reference FASTA, finding a region with no variants? What would be an efficient way to execute this task? </p>
   776: <p>If you store the positions of your random SNPs in a <code>sort-bed</code> sorted BED file, you can filter with BEDOPS <code>bedmap</code>:</p>
  
  <pre><code>$ bedmap --count --echo --range 50000 --delim '\t' SNPs.bed | awk '$1==1' | cut -f2- &gt; SNPs.filtered.bed
  </code></pre>
  
  <p>The file <code>SNPs.filtered.bed</code> will only contain those elements from <code>SNPs.bed</code> which do not overlap within 50kb.</p>
  
  <p>Or you can use BEDOPS <code>bedops</code> and <code>closest-features</code> to do the same thing, but the logic isn't nearly as pretty:</p>
  
  <pre><code>$ bedops --everything --range 0:50000 test.bed | closest-features --dist test.bed - | awk -v FS='|' '{ if ($NR &gt; 50000) { print $1; } }' &gt; SNPs.filtered.bed
  </code></pre>
  
  <p>The result is the same as with the <code>bedmap</code> pipeline. IMO, the <code>bedmap</code> approach is cleaner.</p>
  
  <p>Once you have <code>SNPs.filtered.bed</code>, you could convert it to FASTA via <code>samtools faidx</code> and <code>bed2faidx.pl</code> or similar. Then you have FASTA upon which you can apply your random perturbations.</p>
   77: <p>As others have said, doing whole-genome (or even whole-chromosome) alignments is the wrong solution. Simply track where you're creating SNVs. If you wrote your SNV locations to a sorted BED file you could use something like the following:</p>
  
  
  
  <pre><code>#!/usr/bin/env python
  minWidth = 100000  # Only report regions of &gt;=100kb
  last = [None, None]
  for line in open("some bed file"):
      cols = line.strip().split()
      if cols[0] != last[0]:
          if int(cols[2]) &gt;= minWidth:
              print("{}:1-{}".format(cols[0], cols[2]))
      else:
          if int(cols[2]) - last[1] &gt;= minWidth
              print("{}:{}-{}".format(cols[0], last[1], cols[2]))
      last = [cols[0], int(cols[2])]
  </code></pre>
  
  <p>This will miss the last part of each chromosome/contig, since it doesn't know how long they are, but you should be able to figure out how to handle that. The memory requirements of this are trivially small.</p>
  
  <p>I suspect there's also a bedtools or bedops command that will directly give you the distance to the next downstream region. One could alternatively parse that (it's all the above code is doing, but using <code>bedtools</code> or <code>bedops</code> would produce more validated results than a quickly thrown together script).</p>
  