How to generate random short reads from a reference genome

[wilflugo]

Hi,

I am in the process of generating accuracy benchmarks against different assembly algorithms/tools (including one of my own). Is there an established procedure or common tool to generate these reads from a reference genome?

I know bedtools is able to generate reads based on random positions in the genome, but from what I read seems they will be exact matches (which are randomly placed) against the genome. http://seqanswers.com/forums/images/icons/icon8.gif

I am looking for a tool that is able to add 'noise' to these generated reads in the form of insertions, deletions or mutations on single or few nucleotides.

Such a tool or procedure exist?

Thanks in advance.

[Brian Bushnell]

I wrote a program for that purpose; it's part of BBTools. Basic usage:

randomreads.sh ref=genome.fasta out=reads.fq len=100 reads=10000

You can specify paired reads, an insert size distribution, read lengths (or length ranges), and so forth. But because I developed it to benchmark mapping algorithms, it is specifically designed to give excellent control over mutations. You can specify the number of snps, insertions, deletions, and Ns per read, either exactly or probabilistically; the lengths of these events is individually customizable, the quality values can alternately be set to allow errors to be generated on the basis of quality; there's a PacBio error model; and all of the reads are annotated with their genomic origin, so you will know the correct answer when mapping.

For usage information, run the shellscript with no arguments (or open it with a text editor).

I also have a couple of programs for grading sam files generated using these reads by parsing the read names (samtoroc.sh and gradesam.sh).

[HESmith]

I want to generate synthetic reads to benchmark some variant-calling pipelines. The strategy is to introduce variants into the reference genome, generate random reads from this reference, then align to the non-variant reference. BBTools randomreads.sh seems perfect for this application, esp. the ability to introduce errors into the reads. Two questions:

1) I'd like to generate ~40M 50bp SE reads. Is there an upper limit on the number of independent reads that can be produced?

2) I'd also like to create heterozygous variant reads (by combining the output from a second variant reference). I noticed that randomreads produces the identical output (chromosome/position) by default. Is there a way to change the random generator (e.g., by specifying a different seed) to obtain different output reads?

Thanks,
Harold

[GenoMax]

Harold: Take a look at ART as an additional option. http://www.niehs.nih.gov/research/re...tatistics/art/

[Brian Bushnell]

1) I'd like to generate ~40M 50bp SE reads. Is there an upper limit on the number of independent reads that can be produced?

There's no limit; just set "reads=40000000 length=50".

2) I'd also like to create heterozygous variant reads (by combining the output from a second variant reference). I noticed that randomreads produces the identical output (chromosome/position) by default. Is there a way to change the random generator (e.g., by specifying a different seed) to obtain different output reads?

Yes - "seed=-1" will use a random seed; any other value will use that specific number as the seed.

[HESmith]

Great! Thanks for the quick reply, Brian.

[wilflugo]

I was able to generate all the reads I want with more variabled that I could have asked for. Thanks!. This tool is great!.

However, how I can know the origins of a read?

For example, now I have:
>0
ATACTAAGCGATAGACAAATGACTTTTGCTTCCCTACCATTGA

How do I know which specific scaffold and specific position on the genome this read "mutation" was based from?

From the documentation it says "Read names indicate their genomic origin." However on all my runs the read names are just numbers between 0 and reads-1.

[HESmith]

It looks like your output is in FASTA format. To get chromosome/contig and position information in the identifier, output in FASTQ format. If necessary, you can convert to FASTA using the reformat command.

[Brian Bushnell]

Hmmm....

randomreads.sh reads=4 out=random.fq

Gave me:

@0_chr1_0_3026641_3026740_3018641_Ecoli
GTTCGCCGGGGGCAGCAAACTCAATGCTACACCAACCCGTACCGATAAAAAGATTGCCATTTCCTTACAGGATCTGGAACTGGACTGGGTTGACTGGGAT
+
====@==@AAD?>D@@==DACBC?@BB@C==AB==A@D>AD==?CB==@=B?=A>D?=DB=?>>D@EB===??=@C=?C>@>@B>=?C@@>=====?@>=
@1_chr1_0_831114_831213_823114_Ecoli
TCCGTGGGCGCGGGTAAAAATCTCATCCAGTCCGGCCTGCACTTTTAACGGATGATTAGCTTTTTGCCGCCAGTCATTGAAATCACCCGCCACCAATACC
+
????@CCBECFFGGFEDFGEDFEDGEDFEFEFEFCFGFEEDECEFFADEEEEGCFDFEEEEBCCEFFFFDEDDFEECCFEDEGFFF@BFEDADAC????@
@2_chr1_1_15074406_15074505_4818091_Soneidensis
AGGAAAATGTCGCGTTATCGCCTCTACTTGCCTGTATTATCGAGGAGCAGCAACGGGTACATGGTCATGAGCAAAAAATTCAATTATTTTCTGACGATAC
+
@@BCA@FGEHFGGFGHHFHGEHGFFFHEGFGHGFGEEFEEFFFGFECDFFDHHEDFFDGFEGEEGDDGDEEEECGEDDECFGCFFEFDGFEFFDBAAB@@
@3_chr1_0_4666555_4666654_18580_Evietnamensis
ATTTTATAAATATAATTCCATTATTAACTATTTACTGCTTAAAAAATCATTTACAGCCAAAAATAAGTAAATTATTTACATTTTCAGACTGCTTGCAGCT
+
<<<>>@@BABCCDBDDBDBDDCABBBDB?BCCBBBCBBABC@@CBA@C@C=A=B?BCA<?AA@?B>@B?>CA@@BABBA>@@?BB?=CA<A=A<?><<><

The headers are a little annoying, but take this read, for example:

@2_chr1_1_15074406_15074505_4818091_Soneidensis

That's read 2; it came from the minus strand; and starts at position 4818091 in the sequence named Soneidensis. The stop location of the read is harder to figure out, but is 4818091+(15074505-15074406)=4818190.

I just tested it and it looks like with fasta output you don't get the custom names. In that case, generate fastq output, then you can change it to fasta format like this, keeping the names:

reformat.sh in=random.fq out=random.fa

I will fix that oversight soon.

Edit: Ninja'd by HESmith!

[HESmith]

Even the blind squirrel finds the occasional nut…

Does this mean I've earned my first milli-bushnell?

-Harold

[Brian Bushnell]

Yes, I award you +0.001

[wilflugo]

Thank you guys!. That was the problem. I was generating the output in FASTA.
thanks again.

[wilflugo]

Sorry to continue bothering you guys with this. However, I have been using randomreads.sh from BBmap a lot recently and want to see if what I need can be accomplished.

I am trying to get a based set of reads that are an exact match against the reference genome (no SNPs, no quality issues, just exact chunks of the genome). Since what I am doing is benchmarking different tools accuracy (including one of my own) I want to have a base test case where all tools can match all the reads. Then I would like to start adding SNPs to the reads generation to see how each tool behaves.

Seems I have not being able to generate a set of reads which could be exact matches. This is what I am trying to do:

$ ./randomreads.sh ref=<my_reference> maxsnps=0 maxinss=0 maxdels=0 maxsubs=0 adderrors=false out=reads.fastq reads=1000 minlength=18 maxlength=55 seed=-1

Seems that around 45% of the reads generated are exact matches (calculated by doing grep file matches of the generated read against the reference) but not all of them are. Are there any other parameters that could be added to for exact reads? I could create a simple program myself but I want to continue using random reads for adding SNPs once my base is obtained.

I am assuming that after I can get an exact read then I will just be adding SNPs (using maxsnps parameter). Correct?

[HESmith]

'grep' doesn't wrap lines, and most references have a defined number of characters per line. I suspect that your unmatched reads are split across two lines. Check and let us know.

[wilflugo]

Thanks, yes you are correct. However, I forgot to mention that for the purpose of this exercise I trimmed all scaffolds in the reference genome to have their sequence on a single line.

The reason I did that was to understand why the tools I am benchmarking produced so low accuracy when compare to the start position on the read. When I performed the grep exercise I found the the tools results matched perfectly on the reads that were exact, but bad on the non exact one. Since I want to get the same baseline I want to have a set of perfect matches.