Answers.Body: <p>If you look for a program which would randomly introduce SNPs + short indels and then would save everything into a VCF file, <a href="https://github.com/nh13/DWGSIM" rel="nofollow noreferrer">DWGsim</a> or <a href="https://github.com/seqan/seqan/blob/master/apps/mason2/README.mason_variator" rel="nofollow noreferrer">Mason Variator</a> could be a good choice. Then you can create a corresponding Chain file using <code>bcftools consensus -c</code> and transform various formats between these two coordinate systems using <a href="http://crossmap.sourceforge.net/" rel="nofollow noreferrer">CrossMap</a>.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:12:38.753Z
Answers.Text: How do these algorithms generate an accurate VCF after randomly introducing SNPs and InDels?
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:21:16.283Z
Answers.Text: The SNPs / InDels are generated directly by these programs so the programs know the coordinates, etc.
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:29:36.160Z
Answers.Text: I don't understand how they keep track of these though. If you randomly insert an InDel, it changes the entire reference sequence.
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T21:47:58.030Z
Answers.Text: @ShanZhengYang Mason probably generates an initial read in the reference coordinate. If the initial read overlaps with an indel in VCF, it modifies the read sequence to include the indel. This way, the final read uses the reference coordinate and contains the indel. You don't need to think about another coordinate system at all.
Answers.UserDisplayName: user172818
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-20T20:04:39.883Z
Answers.LastActivityDate: 2017-06-20T20:04:39.883Z
Answers.UserDisplayName: Karel Břinda
Answers.Body: <p>I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences.</p>
  
  <p>The trick has always been to sort the variants by genomic coordinate, <strong>apply the variant with the largest coordinate first</strong>, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate.</p>
  
  <p>Indels and other variants that affect sequence structure only affect <em>subsequent</em> coordinates. So going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T04:21:49.597Z
Answers.Text: "going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence" This is ingenious, thanks! Remaining question: what would be the most efficient way to create the VCF? @DevonRyan above has recommended the approach of writing each variant into a VCF when you manipulate the sequence. Is this your approach, or would you suggest something different?
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T08:06:04.697Z
Answers.Text: That's going to give you the quickest results, because there are existing tools available for creating modified sequences from VCF files.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T10:12:09.940Z
Answers.Text: Separate question: Have you ever tried the above with structural variations? Let's say I wanted to simulate a small structural variation between chromosomes, like a chromosomal translocation....in this more complex example, it becomes more difficult to keep track of changes.
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T17:08:21.593Z
Answers.Text: @ShanZhengYang Nope, I've only ever done SNVs, indels, and inversions. Nothing more complex.
Answers.UserDisplayName: Daniel Standage
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T02:49:42.480Z
Answers.LastActivityDate: 2017-06-22T20:27:43.903Z
Answers.LastEditDate: 2017-06-22T20:27:43.903Z
Answers.UserDisplayName: Daniel Standage
Body: <p>I'm beginning with the reference genome in FASTA format, hg19. I am reading the sequence into a Python dictionary with <code>BioPython</code>:</p>
  
  <pre><code>genome_dictionary = {} 
  for seq_record.id in SeqIO.parse(input_fasta_file, "fasta"):
      genome_dictionary[seq_record.id] = seq_record
  </code></pre>
  
  <p>This creates a Python dictionary with keys as chromosome names and values as strings from the FASTA. </p>
  
  <p>My goal is to take the FASTA, and simulate artificial SNPs and InDels along the FASTA to benchmark various algorithms. If I was only creating SNPs in the FASTA, this would be a trivial problem. The reference FASTA has the original coordinates, and the SNPs would also be in a coordinate system which is the same. </p>
  
  <p>However, when I randomly insert an insertion (e.g. an insertion of 12 bp), I move all events into a new coordinate system such that </p>
  
  <pre><code>new coordinate system = old coordinate system + 15
  </code></pre>
  
  <p>This becomes very difficult to keep track off throughout the entire genome after several interations. I therefore cannot keep track of the new and old coordinate system without running an algorithm and creating a VCF, which may have errors. This defeats the purpose of the simulation.</p>
  
  <p>The output would be either a BED file or VCF which has keep track of the changes.</p>
  
  <p>What method/data structure is used to solve this problem? Would LiftOver work? Maybe CrossMap? </p>
  
  <p><a href="http://genome.sph.umich.edu/wiki/LiftOver" rel="nofollow noreferrer">http://genome.sph.umich.edu/wiki/LiftOver</a>
  <a href="http://crossmap.sourceforge.net/" rel="nofollow noreferrer">http://crossmap.sourceforge.net/</a></p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T19:27:47.220Z
Comments.Text: What are you using to create the changes? Why not just keep track of their location with it?
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T19:30:34.163Z
Comments.Text: @DevonRyan These are just string manipulations, e.g. delete 10 characters is a 10 bp deletion. Keeping track isn't trivial----a delete of 10 bp means all bases must be moved over 10 from the original location.
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T19:35:37.763Z
Comments.Text: Very informally, I'd just suggest that many bioinformatics file formats just don't really handle changed coordinate systems very well. If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off.
Comments.UserDisplayName: Colin D
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T19:39:53.427Z
Comments.Text: "If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off." This doesn't seem possible
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T19:40:27.387Z
Comments.Text: @ShanZhengYang: Err, that is trivial to track, you only care where the changes are in the original file, not where the positions are in the result.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T20:41:24.813Z
Comments.Text: @DevonRyan Perhaps you could explain? I'm randomly inserting SNPs, insertions and deletions in the FASTA. After +3 million events, I need to know how the new coordinates relate to the old coordinates, e.g. a VCF. This sounds quite algorithmically difficult to track to me. How would you do this?
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T20:46:35.367Z
Comments.Text: @ShanZhengYang: As you write the the output, track (1) where you are and (2) whether you insert/delete/mutate the base. As you make a modification, write it to a VCF or other file. You don't need to do anything other than that.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T20:50:07.030Z
Comments.Text: @DevonRyan This still doesn't make sense. So, randomly draw all of the positions for the indels/SNPs. Then sort them. Then implement them one by one, keeping track of the changes made?
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-20T22:06:49.757Z
Comments.Text: @DevonRyan "As you make a modification, write it to a VCF or other file." I think these modifications must be in order. Otherwise, it doesn't work. Correct?
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T06:32:37.573Z
Comments.Text: @ShanZhengYang: You can define at each base what, if any, kind of variant to add there. Yes, you MUST do this in order. Whether you predefine where the variants will be (modify that from the back then) or do it as you write out the sequence (modify if from the front then) is entirely up to you and makes no real difference.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T09:42:39.950Z
Comments.Text: "with keys as chromosome names and values as strings": I think that the code you show, the values are not simple strings, but `SeqRecord` objects.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T10:00:44.627Z
Comments.Text: @bli You're correct---I should edit.
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T10:11:22.610Z
Comments.Text: @DevonRyan This is making sense now. Thank you for the help. A separate question is if I wanted to simulate a small structural variation between chromosomes, like a chromosomal translocation....in this more complex example, it becomes more difficult to keep track of changes.
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-21T10:21:11.683Z
Comments.Text: @ShanZhengYang: Yeah, don't go crazy with those.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-22T15:49:56.310Z
Comments.Text: @DevonRyan What do you do with overlapping intervals for events? Here, the order will matter. e.g. insertion and SNP vs. SNP and insertion
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-22T16:03:17.107Z
Comments.Text: Any one location will only ever have one (or zero)  modifications.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-23T00:19:50.900Z
Comments.Text: @DevonRyan My question may have been unclear. You sort the events. When you do this, it's possible several events will occur at the same (or mre likely, overlapping) genomic coordinates. For example, an insertion from [5,15] and a SNP at 10.
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-23T06:38:35.780Z
Comments.Text: That would be more conveniently described as an insertion at base 5 (10 bases) and an SNP at base 10. Don't think in terms of InDel end positions, just start and widths. Obviously a deletion can't overlap a SNP, so keep track of that.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-05T10:11:12.593Z
Comments.Text: Note: I see in the documentation that there is a `SeqIO.to_dict` function that seems to do your data reading in one step
Comments.UserDisplayName: bli
Tags: <python><bed><structural-variation><indel><simulated-data>
Title: After artificially creating events in a FASTA file, how do I keep track of the old coordinates?
ThreadText: # After artificially creating events in a FASTA file, how do I keep track of the old coordinates?
  
  (Asked by: username_0 on 2017-06-20)
  
  I'm beginning with the reference genome in FASTA format, hg19. I am reading the sequence into a Python dictionary with `BioPython`:
  
  ```
  genome_dictionary = {} 
  for seq_record.id in SeqIO.parse(input_fasta_file, "fasta"):
      genome_dictionary[seq_record.id] = seq_record
  
  ```
  
  This creates a Python dictionary with keys as chromosome names and values as strings from the FASTA.
  
  My goal is to take the FASTA, and simulate artificial SNPs and InDels along the FASTA to benchmark various algorithms. If I was only creating SNPs in the FASTA, this would be a trivial problem. The reference FASTA has the original coordinates, and the SNPs would also be in a coordinate system which is the same.
  
  However, when I randomly insert an insertion (e.g. an insertion of 12 bp), I move all events into a new coordinate system such that
  
  ```
  new coordinate system = old coordinate system + 15
  
  ```
  
  This becomes very difficult to keep track off throughout the entire genome after several interations. I therefore cannot keep track of the new and old coordinate system without running an algorithm and creating a VCF, which may have errors. This defeats the purpose of the simulation.
  
  The output would be either a BED file or VCF which has keep track of the changes.
  
  What method/data structure is used to solve this problem? Would LiftOver work? Maybe CrossMap?
  
  <http://genome.sph.umich.edu/wiki/LiftOver>
  <http://crossmap.sourceforge.net/>
  
  **Tags:** python, bed, structural-variation, indel, simulated-data
  
  > **Comments:**
  > **username_1**: What are you using to create the changes? Why not just keep track of their location with it?  
  > **username_0**: @username_1 These are just string manipulations, e.g. delete 10 characters is a 10 bp deletion. Keeping track isn't trivial----a delete of 10 bp means all bases must be moved over 10 from the original location.  
  > **username_2**: Very informally, I'd just suggest that many bioinformatics file formats just don't really handle changed coordinate systems very well. If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off.  
  > **username_0**: "If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off." This doesn't seem possible  
  > **username_1**: @username_0: Err, that is trivial to track, you only care where the changes are in the original file, not where the positions are in the result.  
  > **username_0**: @username_1 Perhaps you could explain? I'm randomly inserting SNPs, insertions and deletions in the FASTA. After +3 million events, I need to know how the new coordinates relate to the old coordinates, e.g. a VCF. This sounds quite algorithmically difficult to track to me. How would you do this?  
  > **username_1**: @username_0: As you write the the output, track (1) where you are and (2) whether you insert/delete/mutate the base. As you make a modification, write it to a VCF or other file. You don't need to do anything other than that.  
  > **username_0**: @username_1 This still doesn't make sense. So, randomly draw all of the positions for the indels/SNPs. Then sort them. Then implement them one by one, keeping track of the changes made?  
  > **username_0**: @username_1 "As you make a modification, write it to a VCF or other file." I think these modifications must be in order. Otherwise, it doesn't work. Correct?  
  > **username_1**: @username_0: You can define at each base what, if any, kind of variant to add there. Yes, you MUST do this in order. Whether you predefine where the variants will be (modify that from the back then) or do it as you write out the sequence (modify if from the front then) is entirely up to you and makes no real difference.  
  > **username_3**: "with keys as chromosome names and values as strings": I think that the code you show, the values are not simple strings, but `SeqRecord` objects.  
  > **username_0**: @username_3 You're correct---I should edit.  
  > **username_0**: @username_1 This is making sense now. Thank you for the help. A separate question is if I wanted to simulate a small structural variation between chromosomes, like a chromosomal translocation....in this more complex example, it becomes more difficult to keep track of changes.  
  > **username_1**: @username_0: Yeah, don't go crazy with those.  
  > **username_0**: @username_1 What do you do with overlapping intervals for events? Here, the order will matter. e.g. insertion and SNP vs. SNP and insertion  
  > **username_1**: Any one location will only ever have one (or zero)  modifications.  
  > **username_0**: @username_1 My question may have been unclear. You sort the events. When you do this, it's possible several events will occur at the same (or mre likely, overlapping) genomic coordinates. For example, an insertion from [5,15] and a SNP at 10.  
  > **username_1**: That would be more conveniently described as an insertion at base 5 (10 bases) and an SNP at base 10. Don't think in terms of InDel end positions, just start and widths. Obviously a deletion can't overlap a SNP, so keep track of that.  
  > **username_3**: Note: I see in the documentation that there is a `SeqIO.to_dict` function that seems to do your data reading in one step  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_6 on 2017-06-21. Score: 5 | Accepted answer
  
  I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences.
  
  The trick has always been to sort the variants by genomic coordinate, **apply the variant with the largest coordinate first**, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate.
  
  Indels and other variants that affect sequence structure only affect *subsequent* coordinates. So going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence.
  
  > **Comments:**
  > **username_0**: "going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence" This is ingenious, thanks! Remaining question: what would be the most efficient way to create the VCF? @username_1 above has recommended the approach of writing each variant into a VCF when you manipulate the sequence. Is this your approach, or would you suggest something different?  
  > **username_7**: That's going to give you the quickest results, because there are existing tools available for creating modified sequences from VCF files.  
  > **username_0**: Separate question: Have you ever tried the above with structural variations? Let's say I wanted to simulate a small structural variation between chromosomes, like a chromosomal translocation....in this more complex example, it becomes more difficult to keep track of changes.  
  > **username_6**: @username_0 Nope, I've only ever done SNVs, indels, and inversions. Nothing more complex.  
  
  ---
  
  ### Answer by username_4 on 2017-06-20. Score: 4
  
  If you look for a program which would randomly introduce SNPs + short indels and then would save everything into a VCF file, [DWGsim](https://github.com/nh13/DWGSIM) or [Mason Variator](https://github.com/seqan/seqan/blob/master/apps/mason2/README.mason_variator) could be a good choice. Then you can create a corresponding Chain file using `bcftools consensus -c` and transform various formats between these two coordinate systems using [CrossMap](http://crossmap.sourceforge.net/).
  
  > **Comments:**
  > **username_0**: How do these algorithms generate an accurate VCF after randomly introducing SNPs and InDels?  
  > **username_4**: The SNPs / InDels are generated directly by these programs so the programs know the coordinates, etc.  
  > **username_0**: I don't understand how they keep track of these though. If you randomly insert an InDel, it changes the entire reference sequence.  
  > **username_5**: @username_0 Mason probably generates an initial read in the reference coordinate. If the initial read overlaps with an indel in VCF, it modifies the read sequence to include the indel. This way, the final read uses the reference coordinate and contains the indel. You don't need to think about another coordinate system at all.  
  
  ---
  
  