如何用軟件模擬NGS數(shù)據(jù)
為了評(píng)價(jià)一個(gè)工具的性能，通常我們都需要先模擬一批數(shù)據(jù)。這樣相當(dāng)于有了參考答案，才能檢查工具的實(shí)際表現(xiàn)情況。因此對(duì)于我們而言，面對(duì)一個(gè)新的功能，可以先用模擬的數(shù)據(jù)測(cè)試下不同工具的優(yōu)缺點(diǎn)。有如下幾個(gè)工具值得推薦一下：
- 'wgsim/dwgsim': 從全基因組中獲取測(cè)序reads
- 'msbar': EMBOSS其中一個(gè)工具，能夠從單個(gè)序列中模擬隨機(jī)突變
- 'biosed': EMBOSS的一個(gè)工具，可以按照我們給定突變位點(diǎn)模擬
- 'ReadSim': 專門用于模擬PacBio/Nanopore這類儀器產(chǎn)生的long read
- 'Art': 目前最復(fù)雜的模擬工具，能夠模擬測(cè)序儀測(cè)序引入的錯(cuò)誤位點(diǎn)
- 'Metasim': 用于模擬宏基因組得到的reads
- 'Polyester': 用于模擬RNA-seq
值得注意的是，這些工具模擬效果是有限，比如建庫(kù)操作中超聲破碎會(huì)出現(xiàn)的誤差就很難模擬。但是最好的用途就是看看不同生物學(xué)事件在數(shù)據(jù)的情況，比如說發(fā)生了“大規(guī)模倒置”的基因組得到的數(shù)據(jù)比對(duì)到參考基因組上會(huì)是什么情況。
使用dwgsim進(jìn)行模擬
wgism和dwgsim能夠根據(jù)參考基因組模擬出測(cè)序reads，主要是二倍體基因組的SNPs和插入缺失(INDEL)多態(tài)位點(diǎn)。wgism容易安裝，但是參考答案是以簡(jiǎn)單的文本格式保存，不容易可視化。dwgsim受wgism啟發(fā)，雖然安裝稍微麻煩了點(diǎn)，但是參考答案是以VCF格式保存，很方便可視化。
# 請(qǐng)先安裝好ncurse
# 安裝dwgsim
mkdir -p ~/scr
mkdir -p ~/.local/bin
cd ~/src
git clone --recursive https://github.com/nh13/DWGSIM.git
cd DWGSIM
make
ln -s ~/src/DWGSIM/dwgsim ~/.local/bin/dwgsim
ln -s ~/src/DWGSIM/dwgsim_eval ~/.local/bin/dwgsim_eval
簡(jiǎn)單地模擬一批數(shù)據(jù)
# efetch 需要用到conda安裝啟動(dòng)
# conda create -n entrez entrez-direct
# conda activate entrez
# 獲取參考基因組
efetch -db=nuccore -format=fasta -id=AF086833 > genome.fa
# 模擬數(shù)據(jù)
~/.local/bin/dwgsim genome.fa data
會(huì)得到如下數(shù)據(jù)
|-- data.bfast.fastq.gz # 用于bfast
|-- data.bwa.read1.fastq.gz # 用于BWA的R1
|-- data.bwa.read2.fastq.gz # 用于BWA的R2
|-- data.mutations.txt
|-- data.mutations.vcf # VCF形式擦
隨后將這批數(shù)據(jù)用BWA比對(duì)，以bcftools檢測(cè)變異和參考答案比較一下。
# conda install bwa samtools bcftools
bwa index genome.fa
bwa mem genome.fa data.bwa.read1.fastq.gz data.bwa.read2.fastq.gz | samtools sort -o data.bwa.bam
samtools mpileup -uf genome.fa data.bwa.bam | bcftools call -mv -o data.bwa.vcf
samtools index data.bwa.bam
利用使用IGV可視化，檢查分析結(jié)果和真集的一致性
IGV檢查
說明samtools+bcftools找變異這個(gè)組合其實(shí)還是靠譜的，至少在動(dòng)植物領(lǐng)域研究里應(yīng)該夠用。
