基因組組裝
基因組組裝一般分為三個層次，contig, scaffold和chromosomes. contig表示從大規(guī)模測序得到的短讀(reads)中找到的一致性序列。組裝的第一步就是從短片段(pair-end)文庫中組裝出contig。進(jìn)一步基于不同長度的大片段(mate-pair)文庫，將原本孤立的contig按序前后連接，其中會調(diào)整contig方向以及contig可能會存在開口(gap,用N表示)，這一步會得到scaffolds,就相當(dāng)于supercontigs和meatacontigs。最后基于遺傳圖譜或光學(xué)圖譜將scaffold合并調(diào)整，形成染色體級別的組裝(chromosome).
目前基于二代測序的組裝存在挑戰(zhàn)：
- 全基因組測序得到的短讀遠(yuǎn)遠(yuǎn)小于原來的分子長度
- 高通量測序得到海量數(shù)據(jù)會增加組裝的計算復(fù)雜性，消耗更高的計算資源
- 測序錯誤會導(dǎo)致組裝錯誤，會明顯影響contig的長度
- 短讀難以區(qū)分基因組的重復(fù)序列
- 測序覆蓋度不均一，會影響統(tǒng)計檢驗和結(jié)果結(jié)果診斷
上述的問題可以嘗試從如下角度進(jìn)行解決
- 短讀長度：可以通過提供更多樣本，并且建庫時保證位置足夠隨機(jī)
- 數(shù)據(jù)集大小: 使用K-mers算法對數(shù)據(jù)進(jìn)行組裝。assembler不再搜尋overlap，而是搜索具有相同k-mers的reads。但是k-mer算法相比較overlap-based算法，靈敏度有所欠缺，容易丟失一些true overlaps。關(guān)鍵在于定義K。注: K-mer表示一條序列中長度為k的連續(xù)子序列,如ABC的2-mer為AB,BC
- 測序錯誤: 必須保證測序結(jié)果足夠正確, 如提高質(zhì)量控制的標(biāo)準(zhǔn)
- 基因組重復(fù)區(qū)： 測序深度要高，結(jié)果要正確。如果repeat短于read長度，只要保證有足夠多且特異的read。如果repeat長于read，就需要paired ends or “mate-pairs”
- 覆蓋度不均一： 提高深度，保證隨機(jī)
- 組裝結(jié)果比較：contig N50, scaffold N50, BUSCO
二代數(shù)據(jù)組裝的算法和工具
基因組組裝的組裝工具主要分為三類：基于貪婪算法的拼接方法，基于讀序之間的重疊序列(overlapped sequence)進(jìn)行拼接的OLC(Overlap-Layout-Consensus)拼接方法和基于德布魯因圖(de bruijn graph)的方法，這三種方法或多或少基于圖論。第一種是最早期的方法，目前已被淘汰，第二種適用于一代測序產(chǎn)生長片段序列，可以稱之為字符串圖(string graph),第三種是目前二代測序組裝基因組的工具的核心基礎(chǔ)，也就是要繼續(xù)介紹的de bruijn圖。
de bruijn圖由兩部分組成，節(jié)點(diǎn)(Nodes)和邊(Edges)，節(jié)點(diǎn)由k-mers組成，節(jié)點(diǎn)之間要想形成邊就需要是兩個k-mers存在K-1個完全匹配。比如說，ACTG, CTGC, TGCC在K=3時的k-mers為ACT,CTG,TGC,GCC，可以表示為ACT -> CTG -> TGC -> GC.
對于de brujin圖而言，冗余序列不會影響k-mers的數(shù)量，比如說ACTG,ACTG,CTGC,CTGC,CTGC,TGCC,TGCC在K=3時依舊表示為ACT -> CTG -> TGC -> GCC。
上面是理想情況，實(shí)際序列中的測序錯誤，序列之間的SNP以及基因組低復(fù)雜度(重復(fù)序列)就會出現(xiàn)如下de brujin圖
用圖的方式表示就是下面情況
組裝軟件的任務(wù)就是從k-mers形成的圖按照一定的算法組裝出可能的序列，根據(jù)"GAGE: A critical evaluation of genome assemblies and assembly algorithms"以及自己的經(jīng)驗，目前二代數(shù)據(jù)比較常用的工具有Velvet, ABySS, AllPaths/AllPaths-LG, Discovar, SOAPdenovo, Minia, spades,Genomic Assemblers這篇文章有比較好的總結(jié)，
- ALLPaths-LG是公認(rèn)比較優(yōu)秀的組裝工具，但消耗內(nèi)存大，并且要提供至少兩個不同大小文庫的數(shù)據(jù)
- SPAdes是小基因組(<100Mb)組裝時的首選
- SOAPdenovo是目前使用率最高的工具(華大組裝了大量的動植物基因組)，效率也挺好，就是錯誤率也高
- Minia是內(nèi)存資源最省的工具，組裝人類基因組contig居然只要5.7G的RAM，運(yùn)行23小時，簡直難以相信。
當(dāng)然工具之間的差別并沒有想象的那么大，也沒有想象中那么小，可能在物種A表現(xiàn)一般的工具可能在物種B里就非常好用，因此要多用幾個工具，選擇其中最好的結(jié)果。
數(shù)據(jù)準(zhǔn)備
這里使用來自于GAGE的金黃色葡萄球菌 Staphylococcus aureusa 數(shù)據(jù)進(jìn)行練習(xí)。一方面數(shù)據(jù)量小，服務(wù)器能承受并且跑得快，另一方面本身基因組就組裝的不錯，等于是考完試能夠自己對答案。
mkdir Staphylococcus_aureu && cd Staphylococcus_aureus
mkdir genome
curl ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/013/425/GCF_000013425.1_ASM1342v1/GCF_000013425.1_ASM1342v1_genomic.fna.gz > genome/Saureus.fna.gz
mkdir -p raw-data/{lib1,lib2}
curl http://gage.cbcb.umd.edu/data/Staphylococcus_aureus/Data.original/frag_1.fastq.gz > raw-data/lib1/frag_1.fastq.gz
curl http://gage.cbcb.umd.edu/data/Staphylococcus_aureus/Data.original/frag_2.fastq.gz > raw-data/lib2/frag_2.fastq.gz
curl http://gage.cbcb.umd.edu/data/Staphylococcus_aureus/Data.original/shortjump_1.fastq.gz > raw-data/lib2/shortjump_1.fastq.gz
curl http://gage.cbcb.umd.edu/data/Staphylococcus_aureus/Data.original/shortjump_2.fastq.gz > raw-data/lib2/shortjump_2.fastq.gz
基因組survey
在正式組裝之前，需要先根據(jù)50X左右的illumina測序結(jié)果對基因組進(jìn)行評估，了解基因組的大小，重復(fù)序列含量和復(fù)雜度。基于這些信息，確定后續(xù)策略以及是否真的需要對該物種進(jìn)行測序。
基因組survery的核心就是使用k-mers對整體進(jìn)行評估，k-mers時基因組里長度為k的子序列，當(dāng)k=17時，ATCG的組合數(shù)就有170億種，也就說理想條件下基因組大小只有超過17Gb才會出現(xiàn)2條一摸一樣的k-mers。比如說有一個長度為14的序列，給定k-mers的k為8，于是能產(chǎn)生7條長度為8的子序列，于是推測基因組大小為7bp，但是這似乎和實(shí)際的14bp偏離有點(diǎn)遠(yuǎn).
GATCCTACTGATGC (L=14), k-mers for k=8
n = (L-k) + 1 = 14 - 8 + = 7
GATCCTAC, ATCCTACT, TCCTACTG, CCTACTGA, CTACTGAT, TACTGATG, ACTGATGC
如果基因組大小為1MB, 那么當(dāng)k-mers的k=18時，會得到(1000000-18)+1=999983個不同的k-mers，與實(shí)際大小偏差僅僅只有0.0017%，也就說基因組越大，預(yù)測就越接近。這是對單條基因組的估計結(jié)果，實(shí)際上高通量測序會得到基因組30X到50X深度的測序結(jié)果，比如說10個拷貝(C)的“GATCCTACTGATGC”在k-mers=8時會有70條子序列，
n = [(L-K) + 1] * C = [(14-8)+1]*10 = 70
為了得到實(shí)際的基因組大小，既需要將70除以拷貝數(shù)10，那么就得到了和之前一樣的預(yù)測值7。當(dāng)然上述都是理想條件，實(shí)際上測序不均一，低復(fù)雜區(qū)域，重復(fù)序列等都會影響預(yù)測結(jié)果。舉個例子，"Genome sequencing reveals insights into physiology and longevity of the naked mole rat"的k=17, k_num=52,143,337,243，測序程度可以通過k-mers深度分布曲線來估計
圖中，深度為1的k-mers所占比例最高，表示絕大多數(shù)的k-mers僅僅出現(xiàn)了幾次，這可能是測序錯誤造成。后續(xù)在depth=20逐漸形成一個峰，說明測序測度大概是20x附近，實(shí)際上是19x有極大值。于是基因組的大小就是"52,143,337,243/19=2744386170", 差不多就是2.74Gb
k-mers一般選擇17即可，對于高度重復(fù)基因組或者基因組過大，可以選擇19甚至31也行。但不是越大越好，因為如果一條reads里有一個錯誤位點(diǎn)，越大的k-mers就會導(dǎo)致包含這個錯誤位點(diǎn)的k-mers個數(shù)增多。
根據(jù)上述的介紹，便可以使用jellyfish統(tǒng)計k-mer，然后用R作圖對基因組進(jìn)行評估。當(dāng)然這類工具其實(shí)已經(jīng)有人開發(fā)，比如說ALLPATHS-LG/FindErrors，它不但能夠修正低質(zhì)量的短讀，還能初步評估基因組，還有GCE(genome characteristics Estimation)，由華大基因開發(fā)出來的一款基因組評估工具等。為了避免重復(fù)造輪子，簡單就用這些工具即可。
使用GCE評估基因組: 先用kmer_freq_hash統(tǒng)計k-mer頻數(shù)
# Staphylococcus_aureus項目根目錄下
mkdir genome_survey && cd genome_survey
## 提供用于read的位置信息
ls raw-data/lib1/frag_*.fastq.gz > genome_survey/reads.list
## k-mer_freq_hash統(tǒng)計
~/opt/biosoft/gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash -k 15 -l genome_survey/reads.list -t 10 -o 0 -p genome_survey/sa &> genome_survey/kmer_freq.log
k-mer_freq_hash運(yùn)行結(jié)束后會有粗略估計基因組大小，粗略估計為4.22Mb。注意，Kmer_individual_num 數(shù)據(jù)用于gce的輸入?yún)?shù)。
隨后用gce程序基于前面的輸出結(jié)果進(jìn)行估計
~/opt/biosoft/gce-1.0.0/gce -f genome_survey/sa.freq.stat -c 16 -g 108366227 -m 1 -D 8 -b 0 > genome_survey/sa.table 2> genome_survey/sa.log
# -c為主峰對應(yīng)depth
# -g使用的就是Kmer_individual_num對應(yīng)值
# -m 選擇估算模型，真實(shí)數(shù)據(jù)選擇1，表示連續(xù)型
在這次的日志文件中有預(yù)測后的結(jié)果4.34Mb，但是根據(jù)NCBI的數(shù)據(jù)，這個物種的基因組大小是2.8M左右。因此使用k-mers通過數(shù)學(xué)方法預(yù)測存在一定的局限性，需要結(jié)合流式細(xì)胞儀和粗組裝的結(jié)果。
雖然也可以使用FindErrors對基因組進(jìn)行評估，但是我實(shí)際使用時出現(xiàn)了各種問題，這里不做介紹。其他的工具也是大同小異，不做額外推薦。
基因組正式組裝
當(dāng)你拿到測序數(shù)據(jù)后，就可以按照如下幾步處理數(shù)據(jù)。第一步是數(shù)據(jù)質(zhì)控控制，這一步對于組裝而言非常重要，處理前和處理后的組裝結(jié)果可能會天差地別；第二步，根據(jù)經(jīng)驗確定起始參數(shù)，如K-mer和覆蓋率；第三步，使用不同軟件進(jìn)行組裝；第四步，評估組裝結(jié)果，如contig N50, scaffold N50, 判斷是否需要修改參數(shù)重新組裝。
原始數(shù)據(jù)質(zhì)量控制
盡管目前的測序技術(shù)已經(jīng)非常成熟，公司提供的數(shù)據(jù)一般都可以直接用于普通的項目(特殊項目如miRNA-seq除外)。但由于 de novo 組裝對數(shù)據(jù)質(zhì)量比較敏感，因此需要通過質(zhì)控來降低偏差。原始數(shù)據(jù)質(zhì)量控制分為四個部分：
- 了解數(shù)據(jù)質(zhì)量: 了解質(zhì)量這一步可以暫時忽略，基本上基因組測序的結(jié)果都能通過FastQC的標(biāo)準(zhǔn)。
- 去接頭和低質(zhì)量reads過濾: 去接頭和低質(zhì)量reads過濾可供選擇的軟件非常之多，如NGSQCToolkit, Trimmomatic, cutadapter, 似乎都是國外開發(fā)的軟件，但其實(shí)國內(nèi)也有一款很優(yōu)秀的工具叫做fastp
- 去除PCR重復(fù): 去重一般都是在比對后根據(jù)位置信息進(jìn)行，沒有基因組的話只能根據(jù)PE的reads是否完全一樣進(jìn)行過濾。從理論上說，測序相當(dāng)于是從基因組上隨機(jī)抽樣，不太可能存在完全一摸一樣的兩條序列。不過貌似只有FastUniq能做這件事情，后來有一個人寫了sequniq。
- reads修正: 除了過濾或修剪低質(zhì)量的reads外，一般而言，還需要對reads中的錯誤堿基進(jìn)行修正。尤其當(dāng)測序的覆蓋度比較高時，錯誤的reads也就越來越多，會對 de novo 組裝造成不良的影響。工具有BLESS2, BFC, Musket等，其中BLESS2的效率最高，效果也不錯。
去接頭和低質(zhì)量reads過濾這一步推薦fastp，主要是因為它基于C/C++，運(yùn)行速度塊。
# 使用, 項目文件夾下
mkdir -p clean-data{lib1,lib2}
~/opt/biosoft/fastp/fastp -i raw-data/lib1/frag_1.fastq.gz -I raw-data/lib1/frag_2.fastq.gz -o clean-data/lib1/frag_1.fastq.gz -O clean-data/lib1/frag_2.fastq.gz
效果非常的驚人，直接干掉了90%的reads，從原來的1,294,104條變成77,375，一度讓我懷疑軟件是否出現(xiàn)了問題，直到我用同樣的代碼處理現(xiàn)在Illumina的測序結(jié)果以及看了FastQC的結(jié)果才打消了我的疑慮，沒錯，以前的數(shù)據(jù)質(zhì)量就是那么差。注，除非是去接頭，否則不建議通過刪除序列的方式提高質(zhì)量。
質(zhì)控另一個策略是對短讀中一些可能的錯誤堿基進(jìn)行糾正，測序錯誤會引入大量無意義的K-mers，從而增加運(yùn)算復(fù)雜度。此處使用BFC對測序質(zhì)量：
~/opt/biosoft/bfc/bfc -s 3m -t 16 raw-data/lib1/frag_1.fastq.gz | gzip -1 > clean-data/lib1/corrected_1.fq.gz
~/opt/biosoft/bfc/bfc -s 3m -t 16 raw-data/lib1/frag_2.fastq.gz | gzip -1 > clean-data/lib1/corrected_2.fq.gz
總之，質(zhì)控的目標(biāo)是在不引入的錯誤的情況下盡量提高整體質(zhì)量，這一步對后續(xù)的組裝影響很大，所以盡量做這一步，除非組裝軟件要求你別做，那你就不要手賤了。
使用不同工具和參數(shù)進(jìn)行組裝
二代組裝可供選擇的工具很多, 但是主流其實(shí)就那么幾個, 所以組裝的時候選擇3~5個工具運(yùn)行比較結(jié)果即可，比如說MaSuRCA
, IDBA-UD, SOAPdenovo2, Abyss, velvet和Spades。當(dāng)然一旦你選擇一個軟件準(zhǔn)備運(yùn)行的時候，你就會遇到參數(shù)選擇問題，比如怎么確定k-mers，組裝軟件最基礎(chǔ)也是最核心的參數(shù)。這里有幾條原則值得借鑒：
- k要大于log4(基因組大小)，如果數(shù)學(xué)不好，無腦選擇20以上
- 盡量減少測序錯誤形成的k-mers, 因為這是無意義的噪音, 也就是要求k不能過大
- 當(dāng)然k也不能太小，否則會導(dǎo)致重復(fù)壓縮,比如說ATATATA，在2kmers的情況下，就只有AT了
- 測序深度越高，K值也就可以選擇的越大
但是說了那么多，你依舊不知道應(yīng)該選擇什么樣的K，如果你的計算資源無限，那么窮舉法最簡單粗暴。如果窮舉法不行，那么建議先用k=21, 55,77 組裝一下contig, 對不同參數(shù)的contig N50有一個大致的了解，然后繼續(xù)調(diào)整。此外還有一個工具叫做KmerGenie可以預(yù)測一個初始值。總之，讓我們先運(yùn)行第一個工具--SPAdes，可通過bioconda安裝。
SPAdes全稱是圣彼得堡基因組組裝工具，包含了一系列組裝工具處理不同的項目，如高雜合度的dipSPAdes，宏基因組的metaSPAdes。官方文檔中以大腸桿菌為例運(yùn)行整個流程，花了將近1個小時。我們的數(shù)據(jù)集比較小，速度會更快
# 項目根文件夾下
mkdir assembly/spades
spades.py --pe1-1 raw-data/lib1/frag_1.fastq.gz --pe1-2 raw-data/lib1/frag_2.fastq.gz --mp1-1 raw-data/lib2/shortjump_1.fastq.gz --mp1-2 raw-data/lib2/shortjump_2.fastq.gz -o assembly/spades/
你會發(fā)現(xiàn)之前說的k-mers在這里根本沒出現(xiàn)，而且用的也是原始數(shù)據(jù)，這是因為spades.py有一個組件BayesHammer處理測序錯誤，并且它是多K類組裝工具(multi-k assembly), 也就是說它會自動選擇不同的K運(yùn)行，從而挑選比較合適的k值，當(dāng)然你還可以自己設(shè)置，比如說-k 21,55,77。最后結(jié)果為糾正后的短讀數(shù)據(jù)，組裝后的contig, 組裝后的scaffold, 不同格式的組裝graph。
同樣運(yùn)行多k-mers運(yùn)行后比較的工具還有IDBA，它也有一系列的工具。IDBA是基礎(chǔ)版，IDBA-UD適用于宏基因組和單細(xì)胞測序的數(shù)據(jù)組裝，IDBA-Hybrid則是基于相似的基因組提高組裝結(jié)果，IDBA-Tran是專門處理轉(zhuǎn)錄組數(shù)據(jù)。對于無參考基因組組裝，作者推薦使用IDBA-UD。
IDBA-UD工具要求將兩個配對的短讀文件合并成一個，我們的原始數(shù)據(jù)需要先用它提供的fq2fa先轉(zhuǎn)換格式
# 項目文件夾下
mkdir -p assembly/idba_ud
~/opt/biosoft/idba/bin/fq2fa --merge <(zcat clean-data/lib1/corrected_1.fq.gz) <(zcat clean-data/lib1/corrected_2.fq.gz) assembly/idba_ud/lib1.fa
~/opt/biosoft/idba/bin/fq2fa --merge <(zcat clean-data/lib2/corrected_1.fq.gz) <(zcat clean-data/lib2/corrected_2.fq.gz) assembly/idba_ud/lib2.fa
idba_ud和k-mers相關(guān)參數(shù)為--mink,--maxk,--step, 通過--read_level_x 傳入不同大小的文庫，也提供了短讀糾正的相關(guān)參數(shù)--no_correct,--pre_correction
~/opt/biosoft/idba/bin/idba_ud -r assembly/idba_ud/lib1.fa --read_level_2 assembly/idba_ud/lib2.fa -o assembly/idba_ud/ --mink 19 --step 10
運(yùn)行結(jié)束后在assembly/idba_ud下會生成一系列的文件，其中結(jié)果文件為contig.fa和scaffold.fa。
最后介紹一個要手動運(yùn)行不同k-mers的工具，如ABySS, 它有一個亮點(diǎn)，就是能夠可以使用多個計算節(jié)點(diǎn)。我們使用k=31進(jìn)行組裝
mkdir -p assembly/abyss
# 增加 /1,/2
sed 's/^@SRR.*/&\/1/' <(zcat raw-data/lib2/shortjump_1.fastq.gz) | gzip > raw-data/lib2/s1.fq.gz
sed 's/^@SRR.*/&\/2/' <(zcat raw-data/lib2/shortjump_2.fastq.gz) | gzip > raw-data/lib2/s2.fq.gz
~/opt/biosoft/abyss-2.0.2/bin/abyss-pe -C assembly/abyss k=31 n=5 name=asm lib='frag short' frag='../../raw-data/lib1/frag_1.fastq.gz ../../raw-data/lib1/frag_2.fastq.gz' short='../../raw-data/lib2/s1.fq.gz ../../raw-data/lib2/s2.fq.gz' aligner=bowtie
注意，首先ABYSS要求雙端測序的reads命名要以/1和/2結(jié)尾，其次第二個文庫才37bp, 所以比對軟件要選擇bowtie,否則你運(yùn)行一定會遇到
histogram xxx.hist is empty的報錯。當(dāng)然到最后，這個問題我都沒有解決掉，所以我放棄了。
雖然看起來abyss用起來很簡單，但其實(shí)背后的工作流程還是比較復(fù)雜，如下是它的流程示意圖
小結(jié)一下，這里用到了spades, idba,abyss三種工具對同一種物種進(jìn)行組裝，得到對應(yīng)的contig結(jié)果，重點(diǎn)在于k-mers的選擇。contig是組裝的第一步，也是非常重要的一步，為了保證后續(xù)搭scaffold和基因組補(bǔ)洞等工作的順利，我們先得挑選一個比較高質(zhì)量的contig。
組裝可視化和評估
理想條件下，我們希望一個物種有多少染色體，結(jié)果最好就只有多少個contig。當(dāng)然對于二代測序而言，這絕對屬于妄想，可以通過一款graph可視化工具bandage來感受一下最初得到的contig graph是多么復(fù)雜。
一般看這圖直觀感受就是怎么那么多節(jié)，這些節(jié)就是造成contig不連續(xù)的元兇。不同組裝工具在構(gòu)建de bruijn graph的差異不會那么大，contig的數(shù)量和大小和不同工具如何處理復(fù)雜節(jié)點(diǎn)有關(guān)。我們希望得到的contig文件中，每個contig都能足夠的長，能夠有一個完整的基因結(jié)構(gòu)，歸納一下就是3C原則:
- 連續(xù)性(Contiguity): 得到的contig要足夠的長
- 正確性(Correctness): 組裝的contig錯誤率要低
- 完整性(Completeness)：盡可能包含整個原始序列
但是這三條原則其實(shí)是相互矛盾的，連續(xù)性越高，就意味著要處理更多的模糊節(jié)點(diǎn)，會導(dǎo)致整體錯誤率上升，為了保證完全的正確，那么就會導(dǎo)致contig非常的零碎。此外，這三條原則也比較定性，我們需要更加定量的數(shù)值衡量，比如說contig數(shù), 組裝的總長度等, N50等。問題來了，什么叫做N50呢，
N50定義比較繞口，有一種只可意會不可言傳的感覺，所以索性看圖
假設(shè)一個基因組的大小為10，但是這個值只有神知道，你得到的信息就是組裝后有3個contig,長度分別為"3,4,1,1"，所以組裝總長度為9。為了計算N50，我們需要先把contig從大到小排列，也就是"4,3,1"。然后先看最大的contig，長度是4，他的長度是不是超過組裝總大小的一半了嗎？如果是，那么N50=4, 4 < 4.5, 不是。 那么在此基礎(chǔ)上加上第二長的contig,也就是4+3=7, 是不是超過一半了？7>4.5, 那么N50=3. 因此，N50的定義可以表述為"使得累加后長度超過組裝總長度一半的contig的長度就是N50"。為了方便管理和使用軟件，建議建立如下幾個文件夾
N50是基于一個未知的基因組得到得結(jié)果，如果基因組測序比較完整，那么就可以計算NG50，也就是"使得累加后長度超過基因組總長度一半的contig的長度就是NG50"。NA50比較稍微復(fù)雜，需要將組裝結(jié)果進(jìn)一步比對到參考基因組上，以contig實(shí)際和基因組匹配的長度進(jìn)行排序計算。
說完N50，我們介紹兩款工具，QUAST和BUSCO。
QUAST使用質(zhì)量標(biāo)準(zhǔn)(quality metrics)來評估不同組裝工具和不同參數(shù)的組裝效果，無論是否有基因組都可以使用。我們分別以有參和無參兩種模式比較Minia,IDBA和SPAdes三個組裝的運(yùn)行結(jié)果
# without reference
quast.py -o compare idba_ud/contig.fa minia/minia.contigs.fa spades/contigs.fasta
# with reference
quast.py -R ../genome/Saureus.fna -o compare idba_ud/contig.fa minia/minia.contigs.fa spades/contigs.fasta
這個結(jié)果非常直觀的告訴我們一個事實(shí)就是spades組裝的contigs`各方面表現(xiàn)都很優(yōu)秀，minia由于內(nèi)存使用率最低，所以組裝效果一般也是可以理解。
BUSCO通過同源基因數(shù)據(jù)庫從基因完整度來評價基因組組裝結(jié)果。BUSCO首先構(gòu)建了不同物種的最小基因集，然后使用HMMER，BLAST,Augustus等工具分析組裝結(jié)果中的同源基因，從而定量評估組裝是否完整。
busco -i assembly/spades/contigs.fasta -o result -l /home/wangjw/db/busco/bacteria_odb9 -m genome -f
運(yùn)行結(jié)果會在當(dāng)前目錄下的run_result生成一些列文件，其中的short_summary_result.txt內(nèi)容如下
# Summarized benchmarking in BUSCO notation for file assembly/spades/contigs.fasta
# BUSCO was run in mode: genome
C:98.6%[S:98.6%,D:0.0%],F:0.0%,M:1.4%,n:148
146 Complete BUSCOs (C)
146 Complete and single-copy BUSCOs (S)
0 Complete and duplicated BUSCOs (D)
0 Fragmented BUSCOs (F)
2 Missing BUSCOs (M)
C值表示和BUSCO集相比的完整度，M值表示可能缺少的基因數(shù)，D則是重復(fù)數(shù)。正所謂沒有比較，就沒有傷害，我們拿之前QUAST對比中表現(xiàn)比較差的minia結(jié)果作為對比。
C:85.1%[S:85.1%,D:0.0%],F:2.7%,M:12.2%,n:148
126 Complete BUSCOs (C)
126 Complete and single-copy BUSCOs (S)
0 Complete and duplicated BUSCOs (D)
4 Fragmented BUSCOs (F)
18 Missing BUSCOs (M)
98% vs 85%, 一下子對比就出來了。綜上，從兩個維度上證明的SPAdes不但組裝效果好，而且基因完整度也高，當(dāng)然它的內(nèi)存消耗也是很嚴(yán)重。這都是取舍的過程。
附錄
參考資料
軟件安裝
由于不同軟件對不同的基因組的適合度不同，一般都需要參數(shù)多個工具的不同參數(shù)，根據(jù)N50和BUSCO等衡量標(biāo)準(zhǔn)選擇比較好的結(jié)果。為了避免后續(xù)花篇幅在工具安裝上，因此先準(zhǔn)備后續(xù)的分析環(huán)境。對于組裝而言，我們需要安裝如下工具:
- 質(zhì)量控制:
- FastQC
- fastp
- BFC
- 主流組裝工具:
- ABySS
- IDBA
- SOAPdenovo2
- Velvet
- Sapdes
- Minia
- Ray
- MasuRCA
- 基因組組裝評價工具
- BUSCO
- Quast
- 基因結(jié)構(gòu)預(yù)測和功能注釋暫時不在考慮范圍內(nèi)
以下操作所用服務(wù)器的基本信息為：Linux的內(nèi)核為3.10.0-693.el7.x86_64， GCC版本為4.8.5。為了方便管理和使用軟件，建議建立如下幾個文件夾, 分門別類的存放不同工具及其源碼。
# 普通用戶
mkdir -p ~/opt/{sysoft,biosoft}
mkdir -p ~/src
# 管理員
sudo mkdir -p /opt/{sysoft,biosoft}
sudo mkdir -p /src
sudo chmod 1777 /opt/biosoft /opt/sysoft /src
系統(tǒng)自帶的GCC版本是4.8，而BLESS2要求4.9+, ABySS要求6.0+，直接編譯這些工具可能會出錯，但直接升級系統(tǒng)的GCC版本可能會影響整體穩(wěn)定性，因此推薦將在opt/sysoft下安裝高版本的GCC。當(dāng)然GCC的版本也不是越高越好，最好和作者開發(fā)的版本一致，也就是他們要求的最低版本。
# gcc,mpfr,gmp,mpc,isl
cd ~/src
wget -4 https://mirrors.tuna.tsinghua.edu.cn/gnu/gcc/gcc-6.4.0/gcc-6.4.0.tar.xz
tar xf gcc-6.4.0.tar.xz
cd gcc-6.4.0
./contrib/download_prerequisites
mkdir build && cd build
../configure --prefix=$HOME/opt/sysoft/gcc-6.4.0 --enable-threads=posix --disable-multilib --with-system-zlib
make -j 8 && make install
根據(jù)我之前關(guān)于GCC編譯的文章，程序編譯不成功大多是因為找不到頭文件(存放在include目錄下)和鏈接庫文件(存放在lib目錄下)，默認(rèn)編譯頭文件只會搜索/usr/include,/usr/local/include, 而鏈接庫文件只會搜索/lib,/usr/lib[64],/usr/local/lib[64]. 為了讓編譯完成的GCC的頭文件和鏈接庫文件能被搜索到，需要在~/.bashrc文件中添加幾個環(huán)境變量:
-
PKG_CONFIG_PATH: 同時添加搜索頭文件和鏈接頭文件的路徑 -
C_INCLUDE_PATH: 編譯時搜索頭文件的路徑 -
LIBRARY_PATH: 編譯時搜索鏈接文件的路徑 -
LD_LIBRARY_PATH: 運(yùn)行時搜索鏈接文件的路徑
即添加如下幾行內(nèi)容到~/.bashrc文件中，并執(zhí)行source ~/.bashrc更新環(huán)境變量。
export PKG_CONFIG_PATH=~/opt/sysoft/gcc-6.4.0/lib64/pkgconfig:$PKG_CONFIG_PATH
export C_INCLUDE_PATH=~/opt/sysoft/gcc-6.4.0/include:$C_INCLUDE_PATH
export LIBRARY_PATH=~/opt/sysoft/gcc-6.4.0/lib64:$LIBRARY_PATH
export LD_LIBRARY_PATH=~/opt/sysoft/gcc-6.4.0/lib64:$LD_LIBRARY_PATH
export PATH=~/opt/sysoft/gcc-6.4.0/bin:$PATH
genome survey工具: 功能都類似，GCE安裝最方便勝出
cd ~/src
wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.0.tar.gz
tar xf gce-1.0.0.tar.gz -C ~/opt/biosoft
組裝軟件種類很多，對于小基因組(<100Mb)而言SPAdes是很好的選擇， 但是對于大基因組就得多試試幾個，比如說MaSuRCA, Discover de novo, Abyss,SOAPdenovo2, IDBA。內(nèi)存不太夠的話可以嘗試Minia。
組裝軟件一：ABySS的安裝依賴boost1.62, OpenMPI, Google/sparsehash, SQLite，且GCC支持OpenMP，因此也就是一個個下載，一個個安裝的過程。
# boost1.62
cd ~/src
wget -4 https://sourceforge.net/projects/boost/files/boost/1.62.0/boost_1_62_0.tar.bz2
tar xf boost_1_62_0.tar.bz2
cd boost_1_62_0
./bootstrap.sh --prefix=$HOME/opt/sysoft/boost-1.62
./b2
# 引入頭文件的路徑為~/src/boost_1_62_0, 引入鏈接庫的路徑為~/src/boost_1_62_0/stage/lib
# openmpi
wget https://www.open-mpi.org/software/ompi/v3.0/downloads/openmpi-3.0.0.tar.gz
tar xf openmpi-3.0.0.tar.gz
cd openmpi-3.0.0
./configure --prefix=$HOME/opt/sysoft/openmpi-3.0.0
make -j 8 && make install
# 在.bashrc中添加環(huán)境變量或手動修改也行
echo 'export PKG_CONFIG_PATH=~/opt/sysoft/openmpi-3.0.0/lib/pkgconfig:$PKG_CONFIG_PATH' >> ~/.bashrc
echo 'export PATH=~/opt/sysoft/openmpi-3.0.0/bin:$PATH' >> ~/.bashrc
# sparsehash
cd ~/src
git clone https://github.com/sparsehash/sparsehash.git
cd sparsehash
./configure --prefix=$HOME/opt/sysoft/sparsehash
make && make install
# sqlite
cd ~/src
wget -4 http://www.sqlite.org/2018/sqlite-tools-linux-x86-3220000.zip
unzip sqlite-tools-linux-x86-3220000.zip
mv sqlite-tools-linux-x86-3220000 ~/opt/sysoft/sqlite3
最后在安裝ABySS時要以--with-PACKAGE[=ARG]形式指定依賴軟件的路徑
cd ~/src
wget -4 http://www.bcgsc.ca/platform/bioinfo/software/abyss/releases/2.0.2/abyss-2.0.2.tar.gz
tar xf abyss-2.0.2.tar.gz
cd abyss-2.0.2
./configure --prefix=$HOME/opt/biosoft/abyss-2.0.2--with-boost=$HOME/src/boost_1_62_0 --with-sparsehash=$HOME/opt/sysoft/sparsehash --with-sqlite=$HOME/opt/sysoft/sqlite3
make && make install
組裝軟件二：SOAPdenovo2，華大出品，目前使用率最高的工具
cd ~/src
git clone https://github.com/aquaskyline/SOAPdenovo2.git
cd SOAPdenovo2
mkdir -p ~/opt/biosoft/SOAPdenovo2
mv SOAPdenovo-* ~/opt/biosoft/SOAPdenovo2/
組裝軟件三: IDBA. de Brujin圖依賴于K-mers的k的選擇，IDBA能夠自動化遞歸使用不同的k進(jìn)行組裝，從而確定最優(yōu)的K。
cd ~/src
git clone https://github.com/loneknightpy/idba.git
idba/build.sh
mv idba ~/opt/biosoft/
組裝軟件四：MaSuRCA，能夠純用二代，也能二代三代測序混合使用，先用 de bruijn 圖構(gòu)建長reads，然后再用OLC算法進(jìn)行組裝
cd src
wget ftp://ftp.genome.umd.edu/pub/MaSuRCA/latest/MaSuRCA-3.2.4.tar.gz
tar xf MaSuRCA-3.2.4.tar.gz
cd MaSuRCA-3.2.4
export DEST=$HOME/opt/biosoft/MasuRCA
./install.sh
質(zhì)控軟件一: 原本是要推薦BLESS2,但是這個軟件在編譯完成后出現(xiàn)各種核心轉(zhuǎn)移的毛病，和我的系統(tǒng)相性太差，于是改用Li Heng的BFC
cd ~/src
git clone https://github.com/lh3/bfc.git
cd bfc
make
mkdir -p ~/opt/biosoft/bfc
mv bcf hash2cnt ~/opt/biosoft/bfc
質(zhì)控軟件二: fastp是一款基于C/C++編寫的工具，速度會比較塊，而且運(yùn)行之后會有比較好看的圖哦
mkdir -p ~/opt/biosoft/fastp
cd ~/opt/biosoft/fastp
wget http://opengene.org/fastp/fastp
chmod a+x ./fastp
評估工具一：Quast, 它通過比較N50,N G50等參數(shù)來評價基因組組裝質(zhì)量.Quast由Python編寫，推薦使用bioconda安裝
conda create --name assembly python=2.7
source activate assembly
conda install quast
評估工具二: BUSCO，這是一個利用進(jìn)化信息從基因完整性角度評估組裝準(zhǔn)確性的工具，推薦使用biconda安裝。
source activate assembly
conda install busco
盡管conda安裝了busco，但是離實(shí)際運(yùn)行還需要添加幾個環(huán)境變量和不同物種的基因數(shù)據(jù)集，請使用printenv確保如下如下幾個路徑都已經(jīng)添加到環(huán)境變量中。
export PATH="/path/to/AUGUSTUS/augustus-3.2.3/bin:$PATH"
export PATH="/path/to/AUGUSTUS/augustus-3.2.3/scripts:$PATH"
export AUGUSTUS_CONFIG_PATH="/path/to/AUGUSTUS/augustus-3.2.3/config/"
之后，根照自己研究的物種在http://busco.ezlab.org/選擇進(jìn)化上接近的評估數(shù)據(jù)集，比如說你如果研究魚，那么"actinopterygii(輻鰭魚類)"就比"metazoa(多細(xì)胞動物)"更加合適.
實(shí)際運(yùn)行時可能還存在鏈接庫無法找尋以至于程序出錯，解決方法就是將相對應(yīng)或著接近的庫拷貝或軟鏈接到
~/miniconda3/env/assembly/lib下。
