Hi all, I need to create simulated paired-end sequence data with fixed read-lengths on each end (e.g., 75mers on each end of a 500bp DNA fragment, a la Illumina). Does anyone know of a reliable simulator that can generate paired-end sequences to a requested depth, with a requested insert size/variance and error rate, for a requested genome in a FASTA file? The output would preferably be two FASTQ files, one for each end.
I can write my own, but do not want to re-invent this boring (though useful) wheel. Any clues?
•
written
by
Aaronquinlan
samtools wgsim does most of what you request:
Usage: wgsim [options] <in.ref.fa> <out.read1.fq> <out.read2.fq>
Options: -e FLOAT base error rate [0.020]
-d INT outer distance between the two ends [500]
-s INT standard deviation [50]
-N INT number of read pairs [1000000]
-1 INT length of the first read [70]
-2 INT length of the second read [70]
-r FLOAT rate of mutations [0.0010]
-R FLOAT fraction of indels [0.10]
-X FLOAT probability an indel is extended [0.30]
-c generate reads in color space (SOLiD reads)
-C show mismatch info in comment rather than read name
-h haplotype mode
Note: For SOLiD reads, the first read is F3 and the second is R3.
Perfect. I hadn't looked in the misc/ directory in awhile and I never saw documentation for this. Thanks Keith!
MetaSim may be a good option. It has platform specific error modeling and that makes it suited for generating realistic input data rather than "perfectly" random reads.
another solid choice, thank you.
You can also try dwgsim. This is a fork of the SAMtools wgsim and its creator is Nils Homer.
Usage: dwgsim [options] <in.ref.fa> <out.bwa.read1.fq> <out.bwa.read2.fq> <out.bfast.fq>
Options: -e FLOAT base error rate [0.020]
-E FILE base/color error rate file
-d INT outer distance between the two ends [500]
-s INT standard deviation [50]
-N INT number of read pairs [1000000]
-1 INT length of the first read [70]
-2 INT length of the second read [70]
-r FLOAT rate of mutations [0.0010]
-R FLOAT fraction of indels [0.10]
-X FLOAT probability an indel is extended [0.30]
-n INT maximum number of Ns allowed in a given read[0]
-c generate reads in color space (SOLiD reads)
-h haplotype mode
pIRS: Profile-based Illumina pair-end reads simulator.
Or ART
Or simNGS
There's more on this OmicsTools page.
Note the difference between Illumina's paired ends (just reading from each end of a clone), and circularized clones (mate pairs), which give longer inserts, but different directions - and probably more artifacts like chimerae.
(BTW, I've written a simulator for 454 data (flowsim), feel fee to contact me if you're interested in seeing this extended to paired end - or rather, mate paired - sequences.)
RandomReads, in the BBMap package, supports paired-ends. For example:
randomreads.sh ref=ref.fa out=reads.fq paired interleaved reads=100k length=150 mininsert=200 maxinsert=400 gaussian
I have started to have the feeling that everything is implemented in the BBMap package :-)
•
written
by
Istvan Albert
That's my ultimate goal... haven't quite reached it yet!
•
written
by
Brian Bushnell
Hi Brian! Is it possible to generate reads in specific intervals? WES-like read simulation?
No, unfortunately not. You'd have to use something like bedtools to pull out the exome fasta using the genome fasta and the bait coordinates, and then use RandomReads on the result. I don't currently have anything to parse bed, but that does seem like a good addition to RandomReads.
Thank you for the fast answer. I'll try the bedtools pre-step. Another issue.. I've realised that in PE mode, the names of the output reads in the two files are not paired, is there any option for this?
Yes - add the flag "illuminanames".
Is that possible to generate RNA-seq reads from BBmap?
Log in to answer this question.
See also the following thread discussing read simulation with quality scores: http://bit.ly/kNePbA