我正在開髮用於宏基因組學的Illumina測序模擬器： InSilicoSeq

它仍處於alpha發行階段，並且仍處於實驗階段，但考慮到它具有多個Fasta和一個豐度文件，

從文檔起：

  iss generate --genomes基因組.fasta-豐度abundance_file.txt \- model_file HiSeq2500-輸出HiSeq_reads  

其中：

 ＃多fasta文件>genome_AATGC ... >genome_BCCGT ......＃冗余文件必須為1！）genome_A 0.2genome_B 0.4 ...  

我設計的功能不是覆蓋範圍大，而是元基因組中的基因組豐富，因此您可能需要做一個一點數學；）

聚酯生物導體包裝可以做到這一點。它說它模擬RNA-seq讀取，但是我不知道這是否與其他NGS讀取真正不同。

它可以使用一系列誤差和偏差模型，或者從數據集中學習它們。

此python腳本採用帶計數的fasta文件和tsv文件，並在fasta文件中打印序列，次數是tsv文件中指定的次數（假設問題的格式）。因此，如果 bar.tsv 和 foo.fasta 將是您的文件：

 來自Bio import SeqIOrepeat = {}，用於打開行（ “ bar.tsv”）：seq_id，coverage = line.split（）repeat [seq_id] = int（coverage）for SeqIO.parse（foo.fasta，“ fasta”）中的seq_record：適用於範圍內的我（repeat.get（ seq_record.name，0））：print（“ >”，seq_record.name，“ _”，i，sep =''）print（seq_record.seq） 

Heng Li（BWA和samtools出名）的 wgsim軟件包是我用來模擬Illumina讀取的工具。它沒有提供任何方便的方法來模擬不同序列之間的差異覆蓋，但是，應該不難多次運行wgsim，從而為每個感興趣的序列生成所需的覆蓋水平。

I會實現一個Python腳本來處理您的測試文件，並為每個序列調用wgsim（使用 subprocess 模塊）。這可能需要您將每個序列放在單獨的文件中。 ：-（

我不知道有任何軟件可以直接執行此操作，但是我會將fasta文件拆分為每個文件一個序列，在BASH中對其進行循環，然後在每個序列上調用ART（或另一個）序列模擬器。

如果可以接受某些隨機性，則可以使用Poisson隨機變量從序列文件中生成讀數。您需要做一些數學運算以確定要使用的lambda值，以使讀取的每個鹼基對的預期覆蓋率與您在文本文件中設置的值相匹配。

例如序列S的長度為1,000，讀取長度為50，插入大小為100。對於S中的每個鹼基b，都產生一個Poisson隨機變量p。然後，您將從基數b到b + 50生成p個讀數。然後，從b + 50 + 100開始生成配對的讀數。

同樣，您必須使用它來弄清楚要使用的lambda，但這將基本為您提供所需的內容，只要您可以確定要定位的範圍是否完全正確每次閱讀。

使用 RNFtools（版本0.3.1起），可以輕鬆地在控制序列覆蓋率的同時模擬NGS讀數。請參閱指南，特別是序列提取部分。

環境準備

首先，安裝 BioConda並添加所需的頻道。然後要么在默認的Conda環境中安裝RNFtools

  conda install rnftools  

，要么創建並激活一個單獨的Conda環境（首選）

  conda create -n rnftools rnftoolssource激活rnftools  

模擬

假設您有一個參考文件 ref.fa 和一個製表符分隔的覆蓋文件 coverage.tsv （例如，您示例中的參考文件） 。然後，以下RNFtools 將完成您想要的工作：

  import rnftoolsimport csvrnftools.mishmash.sample（“ simulation_with_coverage_control” ，reads_in_tuple = 1）fa =“ ref.fa” tsv =“ coverage.tsv”，其中open（tsv）為f：seqname的表= csv.reader（f，delimiter ='\ t'），表中的cov：rnftools .mishmash.DwgSim（fasta = fa，sequence = [seqname]，coverage = float（cov），read_length_1 = 10，＃使用超短讀取的快速測試read_length_2 = 0，）包括：rnftools.include（）規則：輸入：rnftools。 input（） 

保存此文件（）並運行 snakemake 時，RNFtools將使用DWGsim模擬讀取，並定義覆蓋範圍您的文本文件，並將所有模擬的讀數保存在 simulation_with_coverage_control.fq 中。

您可以使用所有參數。特別是，您可以使用其他模擬器（例如，使用 rnftools.mishmash.ArtIllumina 的Art-Illumina）。有關更多信息，請參見 RNFtools文檔。

另一種下一代讀取模擬工具是 gemsim。我還沒有測試過，但是如果有人有任何經驗，我會很感興趣。

您可以使用

  split -a 6 -p'^ >'your_file.fa seq _  

來按順序拆分FASTA文件。然後使用支持覆蓋率的任何現有讀取模擬器（ART，DWGsim等）。如果要混合所有讀取（而不是按原始順序排序），則可以使用RNFtools。

編輯1：

為@terdon指出，先前的命令僅適用於OSX。一個類似的Linux襯板（但使用數字而不是字母的命名方案略有不同）可以

  csplit -f seq_ -n 6 your_file.fa'/ ^ > /'{* }  

要使此命令也能在OS X上運行，需要安裝coreutils（例如，使用brew），然後使用 gcsplit 而不是 csplit 。

編輯2：

一旦FASTA被序列分割，模擬變得簡單明了，可以使用許多不同的方法。我最喜歡的是使用GNU Parallel。想像一下，您將覆蓋範圍放在名為 covs.txt 的文本文件中的不同行上，並且其順序與 your_file.fa 中的序列相同，例如

  4030 ...  

然後，您可以使用DWGsim通過以下方式模擬原始序列的讀取：

  ls -1 seq_ * |粘貼covs.txt-\ |並行-v --colsep'\ t'dwgsim -1 100 -2 100 -C {1} {2} sim_ {2}  

並使用以下命令合併獲得的FASTQ文件：

  cat sim_seq _ *。bwa.read1.fastq >讀取1.fqcat sim_seq _ *。bwa.read2.fastq >讀取2.fq  

一個可能這種方法的危險在於，我們假設seq_ *文件的數量與 covs.txt 中的行數相同，這可能不是正確的（錯誤地）。我們應該在仿真步驟之前進行檢查，例如：

  [[“” $（ls -1 seq_ * | wc -l）“ ==”“ $（cat covs.txt | wc -l）“]] \ ||迴聲“ covs.txt中的行數不正確”  

另一個需要注意的是，模擬讀段不是隨機排列的（按源序列分組）。