Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-19T10:25:25.493Z
Answers.LastActivityDate: 2017-09-19T20:31:34.907Z
Answers.LastEditDate: 2017-09-19T20:31:34.907Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-18T12:40:38.123Z
Answers.LastActivityDate: 2017-09-18T12:40:38.123Z
Answers.UserDisplayName: Devon Ryan
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-14T01:32:50.197Z
Comments.Text: If you are making the random mutations why not just keep track of the mutated locations to locate the regions of interest.  There is no need to do an alignment if you know where the mutations are since you induced the mutations
Comments.UserDisplayName: Bioathlete
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-14T07:48:59.640Z
Comments.Text: Yes, the obvious approach would be to simply keep track of what regions you have modified. Aligning the entire genome to itself just to find a few SNPs is possible, but much, *much* harder than just keeping track of your changes. Please [edit] your question and explain why you don't just want to track them yourself. Also, it would be helpful if you could show your actual code instead of the example given in the documentation.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-14T14:34:00.350Z
Comments.Text: @Bioathlete In that case, let's say I keep track of the changes I make. I could reserve a region whereby I make no changes in this region of course. Given a list of mutations though in a bed-list format within Python (let's say, using a pandas DataFrame, how would you calculate a region with no mutations?)
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-14T14:45:45.823Z
Comments.Text: @terdon Yes, one could reserve a region and not apply mutations here
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-14T17:31:04.850Z
Comments.Text: I would sort the entries by chr and position and then just compare the current position with the next position to see if it is above the threshold of interest (50kb in your example).
Comments.UserDisplayName: Bioathlete
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-14T18:35:58.010Z
Comments.Text: The khmer library supports aligning against a probabilistic assembly graph, and exposes this functionality to Python. I'm not sure this would be a great fit for your problem though.
Comments.UserDisplayName: Daniel Standage
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-15T05:03:58.333Z
Comments.Text: @Bioathlete Sorting in pandas like this and comparing against other positions can be a memory explosion, if I understand you correctly
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-15T07:45:13.830Z
Comments.Text: @ShanZhengYang please [edit] your question so it asks what you are actually trying to do. If you have a bed file of regions, it's a completely different proposition. Please show us your code, explain what exactly you need to do going forwards and we'll see if we can help.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-15T13:12:48.850Z
Comments.Text: There is no memory explosion.  You are storing a list of coordinates and comparing two of them in an ordered fashion.  If you are worried about memory then write the coordinates to a bed file and either use python to walk the file or use bedtools to write out the distance between the coordinates.  I agree with @terdon.  You need to reformat the question this is too broad and seems to be drifting in the comments
Comments.UserDisplayName: Bioathlete
Tags: <phylogenetics><python><reference-genome>
Title: Is it possible to do alignment within Python? Check variants against reference?
ThreadText: # Is it possible to do alignment within Python? Check variants against reference?
  
  (Asked by: username_0 on 2017-09-13)
  
  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?
  
  **Tags:** phylogenetics, python, reference-genome
  
  > **Comments:**
  > **username_1**: If you are making the random mutations why not just keep track of the mutated locations to locate the regions of interest.  There is no need to do an alignment if you know where the mutations are since you induced the mutations  
  > **username_2**: Yes, the obvious approach would be to simply keep track of what regions you have modified. Aligning the entire genome to itself just to find a few SNPs is possible, but much, *much* harder than just keeping track of your changes. Please [edit] your question and explain why you don't just want to track them yourself. Also, it would be helpful if you could show your actual code instead of the example given in the documentation.  
  > **username_0**: @username_1 In that case, let's say I keep track of the changes I make. I could reserve a region whereby I make no changes in this region of course. Given a list of mutations though in a bed-list format within Python (let's say, using a pandas DataFrame, how would you calculate a region with no mutations?)  
  > **username_0**: @username_2 Yes, one could reserve a region and not apply mutations here  
  > **username_1**: I would sort the entries by chr and position and then just compare the current position with the next position to see if it is above the threshold of interest (50kb in your example).  
  > **username_3**: The khmer library supports aligning against a probabilistic assembly graph, and exposes this functionality to Python. I'm not sure this would be a great fit for your problem though.  
  > **username_0**: @username_1 Sorting in pandas like this and comparing against other positions can be a memory explosion, if I understand you correctly  
  > **username_2**: @username_0 please [edit] your question so it asks what you are actually trying to do. If you have a bed file of regions, it's a completely different proposition. Please show us your code, explain what exactly you need to do going forwards and we'll see if we can help.  
  > **username_1**: There is no memory explosion.  You are storing a list of coordinates and comparing two of them in an ordered fashion.  If you are worried about memory then write the coordinates to a bed file and either use python to walk the file or use bedtools to write out the distance between the coordinates.  I agree with @username_2.  You need to reformat the question this is too broad and seems to be drifting in the comments  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_5 on 2017-09-18. Score: 3 | Accepted answer
  
  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).
  
  
  ---
  
  ### Answer by username_4 on 2017-09-19. Score: 2
  
  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.
  
  
  ---
  
  