Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-29T22:38:11.933Z
Answers.Text: Could you expand on this? This really isn't a very useful answer in this state. What is ART? Where can I find it? How do I use it? How does iot work? Does it introduce errors? Can we control the error rate? Also, does it really not take multifasta files? If not, looping over several thousand files in bash will be very, very slow. There may be better ways.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T10:49:35.580Z
Answers.Text: >Also, does it really not take multifasta files? If not, looping over several thousand files in bash will be very, very slow. There may be better ways
  
  Where did anyone state that it does not take multifasta? OP wanted a very specific coverage per fasta record.
Answers.UserDisplayName: Gabriel Renaud
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T10:54:55.620Z
Answers.Text: No, I didn't mean explain it to me. I meant please edit your answer so it provides this information. We try to have answers as comprehensive as possible and written so they can stand alone without external references on the Stack Exchange sites. Comments are hard to read, easy to miss and result in the horrible mess of biostars where you have no idea where the actual answer is. Please edit your answer to clarify.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-05-17T16:36:36.997Z
Answers.LastActivityDate: 2022-03-04T09:52:25.653Z
Answers.LastEditDate: 2022-03-04T09:52:25.653Z
Answers.UserDisplayName: Gabriel Renaud
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-29T22:41:43.730Z
Answers.Text: Could you expand this by adding an example command showing how to take an input fasta file and produce the fastq? This feels more like (helpful) directions to where an answer can be found than an answer.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-22T19:13:49.987Z
Answers.LastActivityDate: 2017-05-22T19:13:49.987Z
Answers.UserDisplayName: Daniel Standage
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-29T22:49:44.313Z
Answers.Text: Could you explain what this does and how it works? It looks like you will just be printing the original sequence $coverage times, and that can't be right. The OP needs to simulate reads, not just repeat the same sequence N times.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T15:27:33.627Z
Answers.Text: Ouch, then I misunderstood the question. I though that the sequences are the reads that should be simulated and 40x is laterally I would like to have this sequence 40 times. I edited the answer to at least avoid confusion.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T15:38:48.167Z
Answers.Text: I don't know if you misunderstood, or I did. But one of us did :) Thanks for the edit, at least now the OP can see what this does and decide for themselves.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T17:32:11.053Z
Answers.LastActivityDate: 2017-05-30T15:24:35.630Z
Answers.LastEditDate: 2017-05-30T15:24:35.630Z
Answers.UserDisplayName: Kamil S Jaron
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-29T22:40:18.967Z
Answers.Text: Could you expand on this a little? How would the OP use it on their files? What range of error and bias? How would we use them? What would be a reasonable default to use? Perhaps a minimal working example would be helpful.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T08:57:06.240Z
Answers.LastActivityDate: 2017-05-18T08:57:06.240Z
Answers.UserDisplayName: Ian Sudbery
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T07:06:02.543Z
Answers.LastActivityDate: 2017-05-30T07:06:02.543Z
Answers.UserDisplayName: H. Gourlé
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-23T08:18:58.293Z
Answers.Text: And also will be more realistic to simulate insert distance from Poisson-like distribution! The thing is this is a reinvention of a bycicle. There are a lot of simulation tools.
Answers.UserDisplayName: German Demidov
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-23T15:09:36.160Z
Answers.Text: Wouldn't pulling from a normal distribution be better for the insertion size?  And yeah, I know there are lots of tools but this would be one way of doing it yourself if you were interested.
Answers.UserDisplayName: Greg
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-23T15:35:55.120Z
Answers.Text: from my experience, normal is a bit less flexible for fitting - the density is skewed, eg http://bioinformatics.ucdavis.edu/docs/2016-june-workshop/_images/picard-insertion-size-histogram.png . So probably a Negative Binomial will be a good choice, but I never tried to fit it to the actual insert size densities...
Answers.UserDisplayName: German Demidov
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-17T17:55:24.600Z
Answers.LastActivityDate: 2017-05-17T17:55:24.600Z
Answers.UserDisplayName: Greg
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-21T14:46:48.633Z
Answers.LastActivityDate: 2017-06-21T15:32:54.330Z
Answers.LastEditDate: 2017-06-21T15:32:54.330Z
Answers.UserDisplayName: Karel Břinda
Answers.Body: <p>Go to <a href="https://www.niehs.nih.gov/research/resources/software/biostatistics/art/index.cfm" rel="nofollow noreferrer">ART</a> website. Download your package, binary or source:
  <a href="https://i.sstatic.net/137Um.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/137Um.png" alt="screenshot" /></a>
  Uncompress the .gz file, open the folder &quot;art_bin_MountRainier&quot;, then &quot;examples&quot;. Open the text file &quot;run_test_examples_illumina.sh&quot;. This file contains numerous example of samples codes. That will give you an idea on how you can use the software. Here is an example:</p>
  <pre><code># 4) simulation of paired-end reads of 100bp with the mean fragment size 500 and standard deviation 10
  #   using the built-in seperated quality profiles for A, C, G, and T 
  $art -ss HS20 -i ./testSeq.fa -o ./paired_end_sep -l 100 -f 10 -p -m 500 -s 10 -sp -sam
  #convert both aln files to a bed file
  ../aln2bed.pl paired_end_sep.bed paired_end_sep1.aln paired_end_sep2.aln
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-03-04T15:29:42.903Z
Answers.LastActivityDate: 2022-03-04T15:29:42.903Z
Answers.UserDisplayName: Supertech
Answers.Body: <p>Another next-gen read simulation tool is <a href="https://popmodels.cancercontrol.cancer.gov/gsr/packages/gemsim/" rel="nofollow noreferrer">gemsim</a>.  I haven't tested it, but I would be interested if anyone has had any experience with it.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-11T20:27:23.053Z
Answers.LastActivityDate: 2017-06-11T20:27:23.053Z
Answers.UserDisplayName: winni2k
Answers.Body: <p>You can split your FASTA file sequence-by-sequence using</p>
  
  <pre><code>split -a 6 -p '^&gt;' your_file.fa seq_
  </code></pre>
  
  <p>and then use any existing read simulator supporting coverage (ART, DWGsim, etc.). If you want to have all the reads mixed (not ordered by the original sequence), you can use RNFtools.</p>
  
  <p><strong>Edit 1:</strong></p>
  
  <p>As @terdon pointed out, the previous command works on OS X only. An analogical one liner for Linux (but with a slightly different naming scheme using numbers rather than letters) can be</p>
  
  <pre><code>csplit -f seq_ -n 6 your_file.fa '/^&gt;/' {*}
  </code></pre>
  
  <p>To make this command work also on OS X, one needs to install coreutils (e.g., using brew) and then use <code>gcsplit</code> instead of <code>csplit</code>.</p>
  
  <p><strong>Edit 2:</strong></p>
  
  <p>Once FASTA is split by sequences, the simulation becomes straightforward and many different approaches can be used. My favorite one is using GNU Parallel. Imagine that you have your coverages in a text file called <code>covs.txt</code> on separate lines and in the same order as the sequences in <code>your_file.fa</code>, e.g.,</p>
  
  <pre><code>40
  30
  ...
  </code></pre>
  
  <p>Then you can simulate reads from the original sequences using DWGsim by</p>
  
  <pre><code>ls -1 seq_* | paste covs.txt - \
      | parallel -v --colsep '\t' dwgsim -1 100 -2 100 -C {1} {2} sim_{2}
  </code></pre>
  
  <p>and merge the obtained FASTQ files using:</p>
  
  <pre><code>cat sim_seq_*.bwa.read1.fastq &gt; reads.1.fq
  cat sim_seq_*.bwa.read2.fastq &gt; reads.2.fq
  </code></pre>
  
  <p>One possible danger of this approach is that we have assumed that the number of seq_* files is the same as the number of lines in <code>covs.txt</code>, which might not be true (by mistake). We should check this prior to the simulation step, e.g., by:</p>
  
  <pre><code>[[ "$(ls -1 seq_* | wc -l)" == "$(cat covs.txt | wc -l)" ]] \
      || echo "Incorrent number of lines in covs.txt"
  </code></pre>
  
  <p>Another caveat is that the simulated reads are not in a random order (they are grouped by source sequences).</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-29T22:45:40.923Z
Answers.Text: I am afraid you're not answering the question. First of all, you are using non-standard `sort` options. I assume you use a Mac? The `-p` is a BSD-only option and not available on other *nix systems AFAIK. That said, while this will probably work to split the file on a MacOS system, splitting it is trivial and there are many, many ways of doing so. The hard bit is simulating the reads and you're not actually explaining how the OP can do that.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T14:19:42.717Z
Answers.Text: Thanks for the edit. However, the main issue here is that your answer is not answering the question asked. You are answering "how can I split a multifasta file into many single-sequence files" but the question is about simulating reads. Splitting the file might be *part of* an answer, or it might not, but it certainly isn't an answer.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T14:31:03.860Z
Answers.Text: Oh, and for a portable way without using any of the many specialized tools for this, you can use `awk`. Something like `awk -vRS=">" -F'\n' 'NR>1{print ">"$0 >> $1".fa"}' file.fa`. Your `csplit` (and your `split`) approach is indeed simpler though.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-29T19:46:25.187Z
Answers.LastActivityDate: 2017-05-30T18:09:35.290Z
Answers.LastEditDate: 2017-05-30T18:09:35.290Z
Answers.UserDisplayName: Karel Břinda
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-17T17:48:59.563Z
Comments.Text: What is the read length you're using?  How long are the sequences?  Do you need to hit the coverage target exactly or with some probability?
Comments.UserDisplayName: Greg
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-30T15:34:04.493Z
Comments.Text: I would add a few more questions to greg's. Do I understand right that you would like to simulate sequencing of the template sequence from the file? So do the sequences represent genomes? Would you like to consider amplification bias? Sequencing errors? What sequencing platform you would like to simulate?
Comments.UserDisplayName: Kamil S Jaron
Tags: <fasta><ngs><simulated-data>
Title: How to simulate NGS reads, controlling sequence coverage?
ThreadText: # How to simulate NGS reads, controlling sequence coverage?
  
  (Asked by: username_0 on 2017-05-17)
  
  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?
  
  **Tags:** fasta, ngs, simulated-data
  
  > **Comments:**
  > **username_1**: What is the read length you're using?  How long are the sequences?  Do you need to hit the coverage target exactly or with some probability?  
  > **username_2**: I would add a few more questions to username_1's. Do I understand right that you would like to simulate sequencing of the template sequence from the file? So do the sequences represent genomes? Would you like to consider amplification bias? Sequencing errors? What sequencing platform you would like to simulate?  
  
  ---
  
  ## 10 answers
  
  ### Answer by username_3 on 2017-05-17. Score: 2 | Accepted answer
  
  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/>
  
  > **Comments:**
  > **username_4**: Could you expand on this? This really isn't a very useful answer in this state. What is ART? Where can I find it? How do I use it? How does iot work? Does it introduce errors? Can we control the error rate? Also, does it really not take multifasta files? If not, looping over several thousand files in bash will be very, very slow. There may be better ways.  
  > **username_3**: >Also, does it really not take multifasta files? If not, looping over several thousand files in bash will be very, very slow. There may be better ways  Where did anyone state that it does not take multifasta? OP wanted a very specific coverage per fasta record.  
  > **username_4**: No, I didn't mean explain it to me. I meant please edit your answer so it provides this information. We try to have answers as comprehensive as possible and written so they can stand alone without external references on the Stack Exchange sites. Comments are hard to read, easy to miss and result in the horrible mess of biostars where you have no idea where the actual answer is. Please edit your answer to clarify.  
  
  ---
  
  ### Answer by username_7 on 2017-05-30. Score: 7
  
  I am working on a Illumina sequencing simulator for metagenomics: [InSilicoSeq](https://github.com/HadrienG/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 ;)
  
  
  ---
  
  ### Answer by username_5 on 2017-05-22. Score: 4
  
  The [wgsim package](https://github.com/lh3/wgsim) 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. :-(
  
  > **Comments:**
  > **username_4**: Could you expand this by adding an example command showing how to take an input fasta file and produce the fastq? This feels more like (helpful) directions to where an answer can be found than an answer.  
  
  ---
  
  ### Answer by username_2 on 2017-05-18. Score: 4
  
  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)
  
  ```
  
  > **Comments:**
  > **username_4**: Could you explain what this does and how it works? It looks like you will just be printing the original sequence $coverage times, and that can't be right. The OP needs to simulate reads, not just repeat the same sequence N times.  
  > **username_2**: Ouch, then I misunderstood the question. I though that the sequences are the reads that should be simulated and 40x is laterally I would like to have this sequence 40 times. I edited the answer to at least avoid confusion.  
  > **username_4**: I don't know if you misunderstood, or I did. But one of us did :) Thanks for the edit, at least now the OP can see what this does and decide for themselves.  
  
  ---
  
  ### Answer by username_6 on 2017-05-18. Score: 4
  
  The [polyester](http://bioconductor.org/packages/release/bioc/html/polyester.html) 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.
  
  > **Comments:**
  > **username_4**: Could you expand on this a little? How would the OP use it on their files? What range of error and bias? How would we use them? What would be a reasonable default to use? Perhaps a minimal working example would be helpful.  
  
  ---
  
  ### Answer by username_1 on 2017-05-17. Score: 2
  
  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.
  
  > **Comments:**
  > **username_8**: And also will be more realistic to simulate insert distance from Poisson-like distribution! The thing is this is a reinvention of a bycicle. There are a lot of simulation tools.  
  > **username_1**: Wouldn't pulling from a normal distribution be better for the insertion size?  And yeah, I know there are lots of tools but this would be one way of doing it yourself if you were interested.  
  > **username_8**: from my experience, normal is a bit less flexible for fitting - the density is skewed, eg http://bioinformatics.ucdavis.edu/docs/2016-june-workshop/_images/picard-insertion-size-histogram.png . So probably a Negative Binomial will be a good choice, but I never tried to fit it to the actual insert size densities...  
  
  ---
  
  ### Answer by username_9 on 2017-06-21. Score: 2
  
  Simulating NGS reads while controlling sequence coverage is now easy with [RNFtools](http://username_9-brinda.github.io/rnftools/) (from version 0.3.1). See the [tutorial](http://rnftools.readthedocs.io/en/latest/tutorial.html), especially section [Sequence extraction](http://rnftools.readthedocs.io/en/latest/tutorial/02_simulation.html#sequence-extraction).
  
  **Environment preparation**
  
  First, install [BioConda and add the required channels](https://bioconda.github.io/). Then either install RNFtools in the default Conda environment
  
  ```bash
  conda install rnftools
  
  ```
  
  or create and activate a separate Conda environment (preferable)
  
  ```bash
  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:
  
  ```py
  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](http://rnftools.rtfd.org) for more information.
  
  
  ---
  
  ### Answer by username_10 on 2022-03-04. Score: 0
  
  Go to [ART](https://www.niehs.nih.gov/research/resources/software/biostatistics/art/index.cfm) website. Download your package, binary or source:
  [![screenshot](https://i.sstatic.net/137Um.png)](https://i.sstatic.net/137Um.png)
  Uncompress the .gz file, open the folder "art\_bin\_MountRainier", then "examples". Open the text file "run\_test\_examples\_illumina.sh". This file contains numerous example of samples codes. That will give you an idea on how you can use the software. Here is an example:
  
  ```
  # 4) simulation of paired-end reads of 100bp with the mean fragment size 500 and standard deviation 10
  #   using the built-in seperated quality profiles for A, C, G, and T 
  $art -ss HS20 -i ./testSeq.fa -o ./paired_end_sep -l 100 -f 10 -p -m 500 -s 10 -sp -sam
  #convert both aln files to a bed file
  ../aln2bed.pl paired_end_sep.bed paired_end_sep1.aln paired_end_sep2.aln
  
  ```
  
  
  ---
  
  ### Answer by username_11 on 2017-06-11. Score: 0
  
  Another next-gen read simulation tool is [gemsim](https://popmodels.cancercontrol.cancer.gov/gsr/packages/gemsim/). I haven't tested it, but I would be interested if anyone has had any experience with it.
  
  
  ---
  
  ### Answer by username_9 on 2017-05-29. Score: -1
  
  You can split your FASTA file sequence-by-sequence using
  
  ```
  split -a 6 -p '^>' your_file.fa seq_
  
  ```
  
  and then use any existing read simulator supporting coverage (ART, DWGsim, etc.). If you want to have all the reads mixed (not ordered by the original sequence), you can use RNFtools.
  
  **Edit 1:**
  
  As @username_4 pointed out, the previous command works on OS X only. An analogical one liner for Linux (but with a slightly different naming scheme using numbers rather than letters) can be
  
  ```
  csplit -f seq_ -n 6 your_file.fa '/^>/' {*}
  
  ```
  
  To make this command work also on OS X, one needs to install coreutils (e.g., using brew) and then use `gcsplit` instead of `csplit`.
  
  **Edit 2:**
  
  Once FASTA is split by sequences, the simulation becomes straightforward and many different approaches can be used. My favorite one is using GNU Parallel. Imagine that you have your coverages in a text file called `covs.txt` on separate lines and in the same order as the sequences in `your_file.fa`, e.g.,
  
  ```
  40
  30
  ...
  
  ```
  
  Then you can simulate reads from the original sequences using DWGsim by
  
  ```
  ls -1 seq_* | paste covs.txt - \
      | parallel -v --colsep '\t' dwgsim -1 100 -2 100 -C {1} {2} sim_{2}
  
  ```
  
  and merge the obtained FASTQ files using:
  
  ```
  cat sim_seq_*.bwa.read1.fastq > reads.1.fq
  cat sim_seq_*.bwa.read2.fastq > reads.2.fq
  
  ```
  
  One possible danger of this approach is that we have assumed that the number of seq\_\* files is the same as the number of lines in `covs.txt`, which might not be true (by mistake). We should check this prior to the simulation step, e.g., by:
  
  ```
  [[ "$(ls -1 seq_* | wc -l)" == "$(cat covs.txt | wc -l)" ]] \
      || echo "Incorrent number of lines in covs.txt"
  
  ```
  
  Another caveat is that the simulated reads are not in a random order (they are grouped by source sequences).
  
  > **Comments:**
  > **username_4**: I am afraid you're not answering the question. First of all, you are using non-standard `sort` options. I assume you use a Mac? The `-p` is a BSD-only option and not available on other *nix systems AFAIK. That said, while this will probably work to split the file on a MacOS system, splitting it is trivial and there are many, many ways of doing so. The hard bit is simulating the reads and you're not actually explaining how the OP can do that.  
  > **username_4**: Thanks for the edit. However, the main issue here is that your answer is not answering the question asked. You are answering "how can I split a multifasta file into many single-sequence files" but the question is about simulating reads. Splitting the file might be *part of* an answer, or it might not, but it certainly isn't an answer.  
  > **username_4**: Oh, and for a portable way without using any of the many specialized tools for this, you can use `awk`. Something like `awk -vRS=">" -F'\n' 'NR>1{print ">"$0 >> $1".fa"}' file.fa`. Your `csplit` (and your `split`) approach is indeed simpler though.  
  
  ---
  
  