關(guān)于基因組組裝的流程可以參考我的簡(jiǎn)書:http://www.reibang.com/p/e02ecd537c83

前面介紹了SOAPdenovo2的使用，但是由于結(jié)果整體不理想掀抹，嘗試使用另外一個(gè)軟件SPAdes逛艰。

該軟件目前的版本支持paired-end reads, mate-pairs及unpaired reads探赫。SPAdes一般用于小基因組組裝，不適合大型基因組組裝(例:哺乳動(dòng)物)。官方文檔可知第一步和第二步耗內(nèi)存(cpu)，第三步耗磁盤空間(實(shí)際是緩存繁涂，屬于高I/O操作)。

file

SPAdes包含幾個(gè)單獨(dú)的模塊：

(1) BayesHammer:Illumina reads的read糾錯(cuò)二驰。在單細(xì)胞和標(biāo)準(zhǔn)數(shù)據(jù)集上效果也很好扔罪。

(2) IonHammer:IonTorrent(半導(dǎo)體測(cè)序，通過(guò)半導(dǎo)體芯片直接把化學(xué)信號(hào)轉(zhuǎn)為數(shù)字信號(hào))數(shù)據(jù)的read糾錯(cuò)桶雀。

(3) SPAdes: 迭代short reads基因組組裝模塊步势。K的值是根據(jù)read長(zhǎng)度和數(shù)據(jù)集類型自動(dòng)選擇的，也可以自己選擇背犯。

(4) MismatchCorrector:一種工具，可改善所產(chǎn)生的contig和scaffold的錯(cuò)配和較短的插入缺失盅抚。默認(rèn)是關(guān)的漠魏，建議打開(對(duì)于小型基因組而言)。

運(yùn)行SPAdes時(shí)建議運(yùn)行BayesHammer或 IonHammer以獲得高質(zhì)量的組裝(默認(rèn)打開)妄均，如果你使用別的read糾錯(cuò)工具柱锹，該模塊可以關(guān)閉。

1.組裝前準(zhǔn)備

(1)文件夾組織形式:

file

clean_data:存放測(cè)序得到的reads丰包。
fastqc_results:存放兩次fastqc結(jié)果禁熏。(看reads質(zhì)量)
kmergenie_results:存放kmergenie結(jié)果。(基因組調(diào)查)
QC_results:質(zhì)控后的結(jié)果邑彪。
quast_results:quast結(jié)果瞧毙。(組裝結(jié)果評(píng)價(jià))
reference:物種參考基因組。(用于quast，如果沒有參考基因組可以不用)
spades_results:spades組裝結(jié)果宙彪。
tools:組裝所需工具矩动。

(2)所有軟件寫進(jìn)環(huán)境變量。

vim ~/.bashrc
#注:把所需工具寫入到環(huán)境變量中释漆。(fastqc,trimmomatic,kmergenie,SPAdes,quast)
source ~/.bashrc

file

2.read質(zhì)量控制

所用軟件為fastqc和trimmomatic(軟件安裝略)悲没。

cd .../genome_assembly/spades                                   #注:...是省略了前面的路徑,后面也一樣。
touch fq_trim.sh
vim fq_trim.sh
#以下均在fq_trim.sh這個(gè)shell腳本中男图。
trimmomatic=.../genome_assembly/tools/Trimmomatic-0.38/trimmomatic-0.38.jar     #所安裝的trimmomatic所在路徑
#1.first fastqc
for fq_file in clean_data/*                                            
do
fastqc -o ./fastqc_results/first $fq_file                            
done
echo "** first fastqc down! **"
cd ./clean_data
for id in $(seq 501 503)                                        #注:這里根據(jù)你自己data的文件名進(jìn)行設(shè)置示姿。       
do
    file1="PB-${id}_1.clean.fq.gz"                              #測(cè)序得到的兩個(gè)reads文件的文件名。                   
    file2="PB-${id}_2.clean.fq.gz"
    sample_id="PB-${id}"                                   
    if [ -f "$file1" ]                                          #判斷file1是否存在逊笆。
    then
    cd ..
    fq1=./clean_data/$file1
    fq2=./clean_data/$file2
    outdir=.../genome_assembly/spades/QC_results                #trim后存放結(jié)果的文件夾栈戳。
    sample=$sample_id
    #2.trimmomatic(QC)
    time java -jar ${trimmomatic} PE -threads 80 -phred33\
           $fq1 $fq2 \
           $outdir/${sample}.paired.1.fq.gz $outdir/${sample}.unpaired.1.fq.gz \
           $outdir/${sample}.paired.2.fq.gz $outdir/${sample}.unpaired.2.fq.gz \
           ILLUMINACLIP:.../genome_assembly/tools/Trimmomatic-0.38/adapters/TruSeq3-PE-2.fa:2:30:10 \
           SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50 && echo "** fq QC done **"
    cd ./clean_data
    fi
done
#3.second fastqc
cd ..
for fq_pair_file in $outdir/*
do
    fastqc -o ./fastqc_results/second $fq_pair_file
done
echo "** second fastqc down! **"
end=`date +%s`
dif=$[ end - start ]
echo 'the shell script run time is below:'
echo $dif

#運(yùn)行
nohup ./fq_trim.sh >fq_trim.log &                                #放服務(wù)器后臺(tái)運(yùn)行

關(guān)于fastqc結(jié)果解讀請(qǐng)參考簡(jiǎn)書:http://www.reibang.com/p/dacedb7f6e2f

3基因組調(diào)查

對(duì)于較為復(fù)雜的基因組，通常會(huì)在正式組裝前進(jìn)行kmer分析览露，以了解以下信息荧琼。

(1)基因組雜合度。

(2)重復(fù)序列比例差牛。

(3)預(yù)估基因組大小命锄。

(4)測(cè)序深度。

3.1kmergenie評(píng)估基因組

KmerGenie可以進(jìn)行k-mer分析及基因組大小評(píng)估偏化。KmerGenie的最大優(yōu)點(diǎn)在于可以實(shí)現(xiàn)在多個(gè)預(yù)設(shè)k-mer下的自動(dòng)分析脐恩，除了進(jìn)行常規(guī)的k-mer頻數(shù)統(tǒng)計(jì)之外，還能夠基于不同k-mer自動(dòng)計(jì)算基因組大小侦讨，并為基因組組裝評(píng)估一個(gè)最佳組裝k-mer數(shù)值作為備選驶冒。

touch fq_501.txt                                 #存儲(chǔ)fastq文件的路徑。
vim fq_501.txt                                   #編輯fq.txt韵卤，輸入PB_501經(jīng)過(guò)trim后fastq文件路徑(上一篇SOAPdenovo2的使用過(guò)程中我未進(jìn)行trim骗污，這里的話我前面第一步增加了trim)。
touch fq_503.txt                                 
vim fq_503.txt   
#寫kmergenie的腳本
touch kmergenie.sh
vim kmergenie.sh
cd /sevzone/home/ytzheng/genome_assembly/spades
kmergenie fq_501.list -o kmergenie_result/PB_501 -k 140 -l 71 -s 6 -t 12
kmergenie fq_503.list -o kmergenie_result/PB_503 -k 150 -l 100 -s 6 -t 12
#-o:輸出文件存放的目錄沈条。
#-k:最大的k值
#-l:最小的k值
#-s:從最小的k----最大的k需忿，每次增加的k。
#-t:線程數(shù)蜡歹。
#注:剛開始s可以設(shè)大點(diǎn)屋厘，可以根據(jù)判斷出的最佳k值，縮小s再進(jìn)行判斷月而。(如果你設(shè)置的k范圍內(nèi)沒有最佳值汗洒，請(qǐng)調(diào)整 k范圍)

在結(jié)果文件中會(huì)生成report。報(bào)告會(huì)以折線圖的形式給出每種長(zhǎng)度的kmer下父款，估計(jì)的基因組大小溢谤，同時(shí)給出了最佳的kmer(評(píng)估基因組總大小最高)瞻凤。

file

file

kmer頻數(shù)分布圖:

file

3.2jellyfish && GenomeScope評(píng)估基因組

這里請(qǐng)移步到:https://mp.weixin.qq.com/s/3tC4w9Z1LCpPueS6llcK4w

4基因組組裝

4.1SPAdes軟件下載

github：http://cab.spbu.ru/files/release3.14.0/SPAdes-3.14.0-Linux.tar.gz

#解壓
tar -xzf SPAdes-3.14.0-Linux.tar.gz
cd SPAdes-3.14.0-Linux/bin/
ls
#添加到環(huán)境變量(前面部分已經(jīng)添加)
#測(cè)試
spades.py --test

(1)可執(zhí)行文件

file

(2)test結(jié)果(說(shuō)明每個(gè)文件的代表什么)

file

4.2組裝

注:為了運(yùn)行read糾錯(cuò)，reads必須是FASTQ或BAM格式溯香。

4.2.1參數(shù)介紹

spades.py        #查看參數(shù)(怎么用具體可以看官網(wǎng)說(shuō)明鲫构，我僅介紹其中一些參數(shù)。)
#Basic options
-o:指定輸出文件目錄(必需)玫坛。
其他的參數(shù):根據(jù)你是什么數(shù)據(jù)進(jìn)行設(shè)置(如果你的數(shù)據(jù)是宏基因組结笨、質(zhì)粒、單細(xì)胞湿镀、RNA-seq炕吸、 IonTorrent等數(shù)據(jù)，可以指定相應(yīng)參數(shù))勉痴。
#Pipeline options
--only-error-correction:僅進(jìn)行read糾錯(cuò)赫模。
--only-assembler:僅進(jìn)行組裝。
--careful:減少錯(cuò)配和短插入失蒸矛。 讓程序運(yùn)行MismatchCorrector----不建議用于中型和大型基因組瀑罗。官方文檔中可以看到MismatchCorrector占用的disk space非常多。
注:默認(rèn)情況下執(zhí)行read糾錯(cuò)雏掠、組裝斩祭。
--continue:簡(jiǎn)單理解就是你如果跑其中一個(gè)k時(shí)程序停止，SPAdes會(huì)繼續(xù)運(yùn)行該k乡话。
#Input data(僅介紹paired-end)
--pe<#>-1 <file_name>:left reads, <#>代表第幾個(gè)文庫(kù)摧玫。
--pe<#>-2 <file_name>:right reads, <#>代表第幾個(gè)文庫(kù)。
note:只有一個(gè)文庫(kù)绑青，指定#等于1即可诬像。
#Advanced options
-t:線程數(shù)，根據(jù)服務(wù)器有多少線程設(shè)置闸婴。
-m:memory限制(Gb),達(dá)到該值程序終止坏挠，默認(rèn)是250Gb，當(dāng)然不指定它在運(yùn)行時(shí)候也會(huì)給出一個(gè)值邪乍，如果因?yàn)樵撝堤团懿幌氯ソ岛荩阋部梢宰约褐付ā?-k:使用的k-mer用逗號(hào)分隔(所有值必須是奇數(shù)，小于128并且按升序列出)溺欧。--sc：設(shè)置默認(rèn)值為21,33,55。
--cov-cutoff:read coverge的截止值柏肪，可以是float值姐刁、auto、off烦味。默認(rèn)是off聂使。
--phred-offset:堿基質(zhì)量體系壁拉，33或64，可以自己指定柏靶，程序也會(huì)自動(dòng)檢測(cè)弃理。

4.2.2正式組裝

#1.root用戶下設(shè)置打開文件的最大數(shù)量(對(duì)于高I/O,占用的是緩存,不占用CPU,為了提高速度,可通過(guò)線程數(shù)增加來(lái)提高，但這樣會(huì)增加文件打開數(shù)屎蜓，所以需要設(shè)高一點(diǎn))痘昌。
vi /etc/security/limits.conf
xx soft nofile 6000         #這里xx代表的是你用的linux用戶名。
xx hard nofile 6000

#2.非root用戶組裝腳本
touch spades.sh
vim spades.sh
spades.py --pe1-1 QC_results/PB-501.paired.1.fq.gz --pe1-2 QC_results/PB-501.paired.2.fq.gz -t 200 -k 97,107,117,127 -m 600 --careful --phred-offset 33 -o spades_results/PB_501_trim
#--pe1-1炬转、--pe1-2:pair-end測(cè)序得到的兩個(gè)文件名辆苔，這里我用的是經(jīng)過(guò)trim后的數(shù)據(jù)。
#-t:線程數(shù)(對(duì)于CPU密集型的程序扼劈，應(yīng)該使用少一點(diǎn)的線程驻啤，對(duì)于IO密集型任務(wù)，應(yīng)該使用多一點(diǎn)的線程)荐吵。
#-k:kmer選擇的值骑冗，一般不需要設(shè)置，默認(rèn)的話21,33,55,77先煎，但是對(duì)于我的data(水稻)組裝結(jié)果明顯不好贼涩，所以我決定以前面kmergenie預(yù)測(cè)出來(lái)的結(jié)果為參考設(shè)置k值。
#-m:內(nèi)存榨婆。---跑不下去加到跑的下去磁携，不能超。單位為Gb良风，注:free -m可以查看最大的內(nèi)存谊迄。
#--careful:運(yùn)行BayesHammer(read糾錯(cuò))、SPAdes烟央、MismatchCorrector统诺。
#--phred-offset:堿基質(zhì)量體系，33或64疑俭。
#-o:輸出文件目錄粮呢。
#note:spades還有一個(gè)優(yōu)勢(shì)，如果你有其他組裝工具產(chǎn)生的contigs的話可以钞艇，可以使用--trusted-contigs或 --untrusted-contigs選項(xiàng)啄寡。前者是當(dāng)獲得高質(zhì)量組裝中使用的，這些contigs將用于graph的構(gòu)建哩照、間隙填補(bǔ)挺物、處理重復(fù)序列。第二個(gè)選項(xiàng)用于相對(duì)不太可靠的contigs飘弧，這些contigs用于間隙填補(bǔ)识藤、處理重復(fù)序列砚著。

#3.一些問(wèn)題。注:由于我用的SPAdes跑的基因組與官方文檔的相比相對(duì)較大痴昧，在運(yùn)行過(guò)程中出現(xiàn)過(guò)以下問(wèn)題稽穆。

(1).前兩步運(yùn)行時(shí)間相對(duì)還好，占用內(nèi)存峰值約200G(相對(duì)而言是CPU密集型)赶撰，第三步占用的緩存過(guò)多舌镶，導(dǎo)致服務(wù)器的緩存幾乎達(dá)到緩存峰值，這里如果緩存占用太多的話可以用以下腳本釋放緩存(root用戶下)扣囊，但是這樣釋放緩存會(huì)導(dǎo)致程序跑的慢乎折，因此如果服務(wù)器內(nèi)存較大的話還是建議自己手動(dòng)釋放。
touch release_memory.sh
vim release_memory.sh
#以下內(nèi)容在shell腳本中
sleeptime=1800                      #這里可以設(shè)置多少秒清除一次緩存侵歇。
for i in {1..1000}
do
   sync
   echo "sync done"
   echo 3 > /proc/sys/vm/drop_caches
   echo "release done"
   sleep ${sleeptime}
done
#運(yùn)行
nohup ./release_memory.sh >release_memory.log &

(2).第三步為I/O密集型骂澄，我的數(shù)據(jù)跑了很久，為了提高速度惕虑，我把線程數(shù)調(diào)高坟冲，但是調(diào)高的同時(shí)導(dǎo)致我前兩步報(bào)錯(cuò)超過(guò)了文件打開的數(shù)量(默認(rèn)是1024，但還好沒蹦溃蔫，只是在最后結(jié)果給出warning)健提，這里為了防止程序crash掉，我是通過(guò)root用戶來(lái)增加打開文件的數(shù)量(見前面1)伟叛。

file

(3).內(nèi)存不足的問(wèn)題(out of memory)私痹。
free -m  #查看服務(wù)器的內(nèi)存情況，如果還要很多空閑的內(nèi)存统刮，把組裝腳本中的-m參數(shù)調(diào)高紊遵。

file

4.3結(jié)果文件

<output_dir>/corrected/:包含用BayesHammer糾錯(cuò)reads產(chǎn)生*.fastq.gz文件。
<output_dir>/scaffolds.fasta:包含產(chǎn)生的scaffolds侥蒙。
<output_dir>/contigs.fasta:包含產(chǎn)生的contigs暗膜。
<output_dir>/assembly_graph.gfa:包含SPAdes組裝圖和scaffolds的路徑。
<output_dir>/assembly_graph.fastg:包含SPAdes組裝圖鞭衩。(FASTG格式)
<output_dir>/contigs.paths:組裝圖中包含與contigs.fasta對(duì)應(yīng)的路徑学搜。
<output_dir>/scaffolds.paths:組裝圖中包含與scaffolds.fasta對(duì)應(yīng)的路徑。
#note:一般只要scaffolds和contigs就可以了论衍。

5.評(píng)價(jià)組裝結(jié)果

QUAST軟件下載略瑞佩。

5.1參考基因組下載

如果你組裝的物種有參考基因組，可以去下載相應(yīng)的參考基因組坯台。

5.2QUAST使用

touch quast.sh
vim quast.sh
#以下均在shell腳本中炬丸。
cd .../genome_assembly/spades
quast.py -o ./quast_results/PB_501_trim_quast -R ./reference/xxxx.fasta.gz -t 70 ./spades_results/PB_501_trim/scaffolds.fasta ./spades_results/PB_503_trim/scaffolds.fasta
#-o輸出目錄
#-R參考基因組序列
#-t線程數(shù)
#后面為要比較的contig/scaffold所在目錄。

5.3QUAST結(jié)果

具體結(jié)果解讀可以參考:http://blog.sciencenet.cn/blog-3406804-1163959.html

本文由博客一文多發(fā)平臺(tái) OpenWrite 發(fā)布捂人！