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/

  • $\begingroup$ What are you using to create the changes? Why not just keep track of their location with it? $\endgroup$
    – Devon Ryan
    Commented Jun 20, 2017 at 19:27
  • $\begingroup$ @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. $\endgroup$ Commented Jun 20, 2017 at 19:30
  • $\begingroup$ 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. $\endgroup$
    – Colin D
    Commented Jun 20, 2017 at 19:35
  • $\begingroup$ "If you can somehow do the benchmarking without having to manage changed coordinates, you might be better off." This doesn't seem possible $\endgroup$ Commented Jun 20, 2017 at 19:39
  • 1
    $\begingroup$ Note: I see in the documentation that there is a SeqIO.to_dict function that seems to do your data reading in one step $\endgroup$
    – bli
    Commented Jul 5, 2017 at 10:11

2 Answers 2


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.

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.