How To Produce Simulated 'Synthetic' Sequences

Question

[Stefano Berri]

Hi.

I would like to find a program that produces simulated (Illumina GAII) reads starting from a fasta file of the genome and a series of parameters so that I can then test a method I am developing. The features (in descending order) are:

• Uses a fasta file as input (or something easily produced from a fasta file)
• Outputs random reads (or simulate as much as possible any known bias of GAII)
• Is written in (Bio)Perl, Python, ISO C/C++ or fully open source platforms.
• Produces errors (like GAII) at known, ideally tunable, rate
• Allows me to specify depth of coverage and length of reads
• Produces qseq or similar as output
• Allows to produce paired end reads

Often methods papers have some analysis of 'synthetic' or simulated data, but they usually don't bother to publish the program to produce such data. I guess I could write a "quick and dirty" program to do it, but I'd rather not reinvent the wheel if there is something available

Answer 1

[Brad Chapman]

In addition to MetaSim, samtools has wgsim located in the misc directory:

http://samtools.sourceforge.net/

http://bioinformatics.bc.edu/chuanglab/wiki/index.php/How_to_use_SAMtools

It's open source C, and has tunable parameters for error rates, read pair distribution and number of reads generated.

[Stefano Berri]

Thanks. That is just what I was looking for, and I already had it on my computer!

Answer 2

[Michael]

Try MetaSim. That adresses most of your features (except it's written in Java, which is maybe a plus).

Answer 3

[Jonathan Manning]

The author of BFAST has a read simulation program in the dnaa toolkit. A few enhancements over wgsim from samtools.

Answer 4

[User 59]

You might want to look at shuffleseq in EMBOSS.

This will allow you to specify an input fasta file and maintain nucleotide composition at least. You can then shred this into simulated 'reads' which should be relatively straightforward.

You might also be able to use msbar to introduce mutations, but whether they're platform appropriate or not is another question.

EDIT:

Why can't you just use an actual data set? Looking at the Abyss paper, they used experimental data :

sequence data for the genome of an African male individual (HapMap DNA identifier NA18507) (International HapMap Consortium 2003, 2007) from the NCBI short read archive (accession no. SRA000271). The sequence was generated by Illumina, Inc. using their Genome Analyzer platform (Bentley et al. 2008).

as well as synthetic data generated like this:

The first synthetic data set represented all possible error-free 36-mer paired sequences, using a fixed fragment size of 200 bp. We generated simulated reads by sliding a 200 bp window, with a step size of 1 bp, along each chromosome of the reference genome and reporting the first 36 bp and the reverse complement of the last 36 bp. This process produced a data set of perfectly tiled 72-fold read coverage of the reference genome.

[Stefano Berri]

I am afraid this is not what I am looking for. I need that from a series of fasta sequences (the human chromosomes) produces little bits (36-72 Bp long) as if it was "sequencing" the genome

[User 59]

In that case you don't want it to output 'random' reads :)

[Stefano Berri]

@Daniel. I can't use real data because I cannot change the underlying genome. I want to make, for instance, a duplication of a chromosome arm, then simulate "random reads" from that genome (with associated noise) and then use the program I am developing to see if I can pick up the duplication and with which resolution. getting every 200 bp is no good, because then it is trivial (and not realistic)

[steve]

I found this tool which sounds similar to what you are looking for: bamsurgeon. Havent tried it yet though