text: user174: I have a FASTA file with 100+ sequences like this:
  >Sequence1
  GTGCCTATTGCTACTAAAA ...
  >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?
  
  user31: I am not aware of any software that can do this directly, but I would split the fasta file into one sequence per file, loop over them in BASH and invoke ART the sequence simulator (or another) on each sequence.
  For more information about ART, please see their paper here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/
  
  user59: If you're ok with some randomness you can generate reads from your sequence file using a Poisson random variable.  You'll need to do some math to figure out what value of lambda to use in order for the expected coverage at each base pair in your read to match what you set in your text file.
  For example you have a sequence S of length 1,000, a read length of 50, and an insert size of 100.  For each base b in S generate a Poisson random variable p.  You will then generate p reads from base b to b+50.  Then, generate the paired read starting at b+50+100.  
  Again, you would have to play with it to figure out what lambda to use but this would give you basically what you want, as long as you're ok with not having exactly the coverage you're targeting for each read.
  
  user425: Simulating NGS reads while controlling sequence coverage is now easy with RNFtools (from version 0.3.1). See the tutorial, especially section Sequence extraction.
  Environment preparation
  First, install BioConda and add the required channels. Then either install RNFtools in the default Conda environment
  conda install rnftools
  
  or create and activate a separate Conda environment (preferable)
  conda create -n rnftools rnftools
  source activate rnftools
  
  Simulation
  Assume that you have a reference file ref.fa and a tab-separated coverage file coverage.tsv (e.g., those from your example). Then the following RNFtools Snakefile will do the job you want:
  import rnftools
  import csv
  
  rnftools.mishmash.sample("simulation_with_coverage_control", reads_in_tuple=1)
  
  fa = "ref.fa"
  tsv = "coverage.tsv"
  
  with open(tsv) as f:
      table = csv.reader(f, delimiter='\t')
      for seqname, cov in table:
  
          rnftools.mishmash.DwgSim(
              fasta=fa,
              sequences=[seqname],
              coverage=float(cov),
              read_length_1=10, # quick test with supershort reads
              read_length_2=0,
          )
  
  include: rnftools.include()
  rule: input: rnftools.input()
  
  When you save this file (Snakefile) and run snakemake, RNFtools will simulate reads using DWGsim with the coverages defined your text file, and save all the simulated reads in simulation_with_coverage_control.fq.
  You can play with all the parameters. In particular, you can use a different simulator (e.g., Art-Illumina using rnftools.mishmash.ArtIllumina). See the RNFtools documentation for more information.
  
  user235: The polyester bioconductor package can do this. It says it simulates RNA-seq reads, but I don't know if that's really any different from other NGS reads. 
  It can use a range of error and bias models, or learn them from a dataset.
  
  user57: This python script takes a fasta file and tsv file with counts and prints the sequences in the fasta files that many times as it is specified in the tsv file (assuming the format in the question). So if bar.tsv and foo.fasta will be your files:
  from Bio import SeqIO
  
  repeat = {}
  for line in open("bar.tsv"):
      seq_id, coverage = line.split()
      repeat[seq_id] = int(coverage)
  
  for seq_record in SeqIO.parse(foo.fasta, "fasta"):
      for i in range(repeat.get(seq_record.name, 0)):
          print(">",seq_record.name,"_",i,sep='')
          print(seq_record.seq)
  
  
  user96: The wgsim package by Heng Li (of BWA and samtools fame) is my go-to tool for simulating Illumina reads. It doesn't provide any convenient way to simulate differential coverage across different sequences, but it shouldn't be to hard to run wgsim multiple times, generating the desired level of coverage for each sequence of interest.
  I would implement a Python script to slurp up your test file, and call wgsim (using the subprocess module) for each sequence. This will probably require you to have each sequence in a separate file. :-(
  
  user179: 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
  >genome_A
  ATGC...
  >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 ;)
  
original_text: 174: <p>I have a FASTA file with 100+ sequences like this:</p>
  
  <pre><code>&gt;Sequence1
  GTGCCTATTGCTACTAAAA ...
  &gt;Sequence2
  GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA
  ......
  </code></pre>
  
  <p>I also have a text file like this:</p>
  
  <pre><code>Sequence1 40
  Sequence2 30
  ......
  </code></pre>
  
  <p>I would like to simulate next-generation paired-end reads for all the sequences in my FASTA file. For <code>Sequence1</code>, I would like to simulate at 40x coverage. For <code>Sequence2</code>, I would like to simulate at 30x coverage. In other words, I want to control my sequence coverage for each sequence in my simulation.</p>
  
  <p><strong>Q:</strong> What is the simplest way to do that? Any software I should use? Bioconductor?</p>
   31: <p>I am not aware of any software that can do this directly, but I would split the fasta file into one sequence per file, loop over them in BASH and invoke ART the sequence simulator (or another) on each sequence.</p>
  <p>For more information about ART, please see their paper here: <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/" rel="nofollow noreferrer">https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/</a></p>
   59: <p>If you're ok with some randomness you can generate reads from your sequence file using a Poisson random variable.  You'll need to do some math to figure out what value of lambda to use in order for the expected coverage at each base pair in your read to match what you set in your text file.</p>
  
  <p>For example you have a sequence S of length 1,000, a read length of 50, and an insert size of 100.  For each base b in S generate a Poisson random variable p.  You will then generate p reads from base b to b+50.  Then, generate the paired read starting at b+50+100.  </p>
  
  <p>Again, you would have to play with it to figure out what lambda to use but this would give you basically what you want, as long as you're ok with not having exactly the coverage you're targeting for each read.</p>
   425: <p>Simulating NGS reads while controlling sequence coverage is now easy with <a href="http://karel-brinda.github.io/rnftools/" rel="nofollow noreferrer">RNFtools</a> (from version 0.3.1). See the <a href="http://rnftools.readthedocs.io/en/latest/tutorial.html" rel="nofollow noreferrer">tutorial</a>, especially section <a href="http://rnftools.readthedocs.io/en/latest/tutorial/02_simulation.html#sequence-extraction" rel="nofollow noreferrer">Sequence extraction</a>.</p>
  
  <p><strong>Environment preparation</strong></p>
  
  <p>First, install <a href="https://bioconda.github.io/" rel="nofollow noreferrer">BioConda and add the required channels</a>. Then either install RNFtools in the default Conda environment</p>
  
  <pre class="lang-bash prettyprint-override"><code>conda install rnftools
  </code></pre>
  
  <p>or create and activate a separate Conda environment (preferable)</p>
  
  <pre class="lang-bash prettyprint-override"><code>conda create -n rnftools rnftools
  source activate rnftools
  </code></pre>
  
  <p><strong>Simulation</strong></p>
  
  <p>Assume that you have a reference file <code>ref.fa</code> and a tab-separated coverage file <code>coverage.tsv</code> (e.g., those from your example). Then the following RNFtools <code>Snakefile</code> will do the job you want:</p>
  
  <pre class="lang-py prettyprint-override"><code>import rnftools
  import csv
  
  rnftools.mishmash.sample("simulation_with_coverage_control", reads_in_tuple=1)
  
  fa = "ref.fa"
  tsv = "coverage.tsv"
  
  with open(tsv) as f:
      table = csv.reader(f, delimiter='\t')
      for seqname, cov in table:
  
          rnftools.mishmash.DwgSim(
              fasta=fa,
              sequences=[seqname],
              coverage=float(cov),
              read_length_1=10, # quick test with supershort reads
              read_length_2=0,
          )
  
  include: rnftools.include()
  rule: input: rnftools.input()
  </code></pre>
  
  <p>When you save this file (<code>Snakefile</code>) and run <code>snakemake</code>, RNFtools will simulate reads using DWGsim with the coverages defined your text file, and save all the simulated reads in <code>simulation_with_coverage_control.fq</code>.</p>
  
  <p>You can play with all the parameters. In particular, you can use a different simulator (e.g., Art-Illumina using <code>rnftools.mishmash.ArtIllumina</code>). See the <a href="http://rnftools.rtfd.org" rel="nofollow noreferrer">RNFtools documentation</a> for more information.</p>
   235: <p>The <a href="http://bioconductor.org/packages/release/bioc/html/polyester.html" rel="nofollow noreferrer">polyester</a> bioconductor package can do this. It says it simulates RNA-seq reads, but I don't know if that's really any different from other NGS reads. </p>
  
  <p>It can use a range of error and bias models, or learn them from a dataset.</p>
   57: <p>This python script takes a fasta file and tsv file with counts and prints the sequences in the fasta files that many times as it is specified in the tsv file (assuming the format in the question). So if <code>bar.tsv</code> and <code>foo.fasta</code> will be your files:</p>
  
  
  
  <pre><code>from Bio import SeqIO
  
  repeat = {}
  for line in open("bar.tsv"):
      seq_id, coverage = line.split()
      repeat[seq_id] = int(coverage)
  
  for seq_record in SeqIO.parse(foo.fasta, "fasta"):
      for i in range(repeat.get(seq_record.name, 0)):
          print("&gt;",seq_record.name,"_",i,sep='')
          print(seq_record.seq)
  </code></pre>
   96: <p>The <a href="https://github.com/lh3/wgsim" rel="nofollow noreferrer">wgsim package</a> by Heng Li (of BWA and samtools fame) is my go-to tool for simulating Illumina reads. It doesn't provide any convenient way to simulate differential coverage across different sequences, but it shouldn't be to hard to run wgsim multiple times, generating the desired level of coverage for each sequence of interest.</p>
  
  <p>I would implement a Python script to slurp up your test file, and call wgsim (using the <code>subprocess</code> module) for each sequence. This will probably require you to have each sequence in a separate file. :-(</p>
   179: <p>I am working on a Illumina sequencing simulator for metagenomics: <a href="https://github.com/HadrienG/InSilicoSeq" rel="noreferrer">InSilicoSeq</a></p>
  
  <p>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.</p>
  
  <p>From the documentation:</p>
  
  <pre><code>iss generate --genomes genomes.fasta --abundance abundance_file.txt \
      --model_file HiSeq2500 --output HiSeq_reads
  </code></pre>
  
  <p>Where:</p>
  
  <pre><code># 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
  ...
  </code></pre>
  
  <p>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 ;)</p>
  