After artificially creating events in a FASTA file, how do I keep track of the old coordinates?

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?

• What are you using to create the changes? Why not just keep track of their location with it? Jun 20 '17 at 19:27
• @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. Jun 20 '17 at 19:30
• 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. Jun 20 '17 at 19:35
• "If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off." This doesn't seem possible Jun 20 '17 at 19:39
• 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
Jul 5 '17 at 10:11

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.

• "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? Jun 21 '17 at 4:21
• That's going to give you the quickest results, because there are existing tools available for creating modified sequences from VCF files.
– gringer
Jun 21 '17 at 8:06
• 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. Jun 21 '17 at 10:12
• @ShanZhengYang Nope, I've only ever done SNVs, indels, and inversions. Nothing more complex. Jun 21 '17 at 17:08

If you look for a program which would randomly introduce SNPs + short indels and then would save everything into a VCF file, DWGsim or 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.

• How do these algorithms generate an accurate VCF after randomly introducing SNPs and InDels? Jun 20 '17 at 20:12
• The SNPs / InDels are generated directly by these programs so the programs know the coordinates, etc. Jun 20 '17 at 20:21
• I don't understand how they keep track of these though. If you randomly insert an InDel, it changes the entire reference sequence. Jun 20 '17 at 20:29
• @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. Jun 20 '17 at 21:47