How to simulate NGS reads, controlling sequence coverage? I have a FASTA file with 100+ sequences like this:

&gt;Sequence1
GTGCCTATTGCTACTAAAA ...
&gt;Sequence2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA
......


I also have a text file like this:

Sequence1 40
Sequence2 30
......


I would like to simulate next-generation paired-end reads for all the sequences in my FASTA file. For Sequence1, I would like to simulate at 40x coverage. For Sequence2, I would like to simulate at 30x coverage. In other words, I want to control my sequence coverage for each sequence in my simulation.

Q: What is the simplest way to do that? Any software I should use? Bioconductor?
I am working on a Illumina sequencing simulator for metagenomics: InSilicoSeq

It is still in alpha release and very experimental, but given a multi-fasta and an abundance file, it will generate reads from your input genomes with different coverages.

From the documentation:

iss generate --genomes genomes.fasta --abundance abundance_file.txt \
    --model_file HiSeq2500 --output HiSeq_reads


Where:

# multi-fasta file
&gt;genome_A
ATGC...
&gt;genome_B
CCGT...
...

# abundance file (total abundance must be 1!)
genome_A    0.2
genome_B    0.4
...


I didn't design it to work with coverage but rather abundance of the genome in a metagenome, so you might have to do a tiny bit of math ;)