<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>

CC BY-SA 3.0
2017-06-20T20:12:38.753Z
How do these algorithms generate an accurate VCF after randomly introducing SNPs and InDels?
ShanZhengYang
CC BY-SA 3.0
2017-06-20T20:21:16.283Z
The SNPs / InDels are generated directly by these programs so the programs know the coordinates, etc.
Karel Břinda
CC BY-SA 3.0
2017-06-20T20:29:36.160Z
I don't understand how they keep track of these though. If you randomly insert an InDel, it changes the entire reference sequence.
ShanZhengYang
CC BY-SA 3.0
2017-06-20T21:47:58.030Z
@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.
user172818
CC BY-SA 3.0
2017-06-20T20:04:39.883Z
2017-06-20T20:04:39.883Z
Karel Břinda
<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>

CC BY-SA 3.0
2017-06-21T04:21:49.597Z
"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?
ShanZhengYang
CC BY-SA 3.0
2017-06-21T08:06:04.697Z
That's going to give you the quickest results, because there are existing tools available for creating modified sequences from VCF files.
gringer
CC BY-SA 3.0
2017-06-21T10:12:09.940Z
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.
ShanZhengYang
CC BY-SA 3.0
2017-06-21T17:08:21.593Z
@ShanZhengYang Nope, I've only ever done SNVs, indels, and inversions. Nothing more complex.
Daniel Standage
CC BY-SA 3.0
2017-06-21T02:49:42.480Z
2017-06-22T20:27:43.903Z
2017-06-22T20:27:43.903Z
Daniel Standage
<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>

CC BY-SA 3.0
2017-06-20T19:27:47.220Z
What are you using to create the changes? Why not just keep track of their location with it?
Devon Ryan
CC BY-SA 3.0
2017-06-20T19:30:34.163Z
@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.
ShanZhengYang
CC BY-SA 3.0
2017-06-20T19:35:37.763Z
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.
Colin D
CC BY-SA 3.0
2017-06-20T19:39:53.427Z
"If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off." This doesn't seem possible
ShanZhengYang
CC BY-SA 3.0
2017-06-20T19:40:27.387Z
@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.
Devon Ryan
CC BY-SA 3.0
2017-06-20T20:41:24.813Z
@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?
ShanZhengYang
CC BY-SA 3.0
2017-06-20T20:46:35.367Z
@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.
Devon Ryan
CC BY-SA 3.0
2017-06-20T20:50:07.030Z
@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?
ShanZhengYang
CC BY-SA 3.0
2017-06-20T22:06:49.757Z
@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?
ShanZhengYang
CC BY-SA 3.0
2017-06-21T06:32:37.573Z
@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.
Devon Ryan
CC BY-SA 3.0
2017-06-21T09:42:39.950Z
"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.
bli
CC BY-SA 3.0
2017-06-21T10:00:44.627Z
@bli You're correct---I should edit.
ShanZhengYang
CC BY-SA 3.0
2017-06-21T10:11:22.610Z
@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.
ShanZhengYang
CC BY-SA 3.0
2017-06-21T10:21:11.683Z
@ShanZhengYang: Yeah, don't go crazy with those.
Devon Ryan
CC BY-SA 3.0
2017-06-22T15:49:56.310Z
@DevonRyan What do you do with overlapping intervals for events? Here, the order will matter. e.g. insertion and SNP vs. SNP and insertion
ShanZhengYang
CC BY-SA 3.0
2017-06-22T16:03:17.107Z
Any one location will only ever have one (or zero)  modifications.
Devon Ryan
CC BY-SA 3.0
2017-06-23T00:19:50.900Z
@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.
ShanZhengYang
CC BY-SA 3.0
2017-06-23T06:38:35.780Z
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.
Devon Ryan
CC BY-SA 3.0
2017-07-05T10:11:12.593Z
Note: I see in the documentation that there is a `SeqIO.to_dict` function that seems to do your data reading in one step
bli
<python><bed><structural-variation><indel><simulated-data>
After artificially creating events in a FASTA file, how do I keep track of the old coordinates?
# 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**: @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.  
> **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**: @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.  
> **username_0**: @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?  
> **username_1**: @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.  
> **username_0**: @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?  
> **username_0**: @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?  
> **username_1**: @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.  
> **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**: @bli You're correct---I should edit.  
> **username_0**: @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.  
> **username_1**: @ShanZhengYang: Yeah, don't go crazy with those.  
> **username_0**: @DevonRyan 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**: @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.  
> **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? @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?  
> **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**: @ShanZhengYang 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**: @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.  

---

