基因組組裝

基因組組裝一般分為三個層次，contig, scaffold和chromosomes. contig表示從大規(guī)模測序得到的短讀(reads)中找到的一致性序列。組裝的第一步就是從短片段(pair-end)文庫中組裝出contig。進一步基于不同長度的大片段(mate-pair)文庫妻率，將原本孤立的contig按序前后連接弓千，其中會調(diào)整contig方向以及contig可能會存在開口(gap,用N表示)兵扬，這一步會得到scaffolds,就相當于supercontigs和meatacontigs指蚜。最后基于遺傳圖譜或光學圖譜將scaffold合并調(diào)整乞巧，形成染色體級別的組裝(chromosome).

目前基于二代測序的組裝存在挑戰(zhàn)：

• 全基因組測序得到的短讀遠遠小于原來的分子長度
• 高通量測序得到海量數(shù)據(jù)會增加組裝的計算復雜性，消耗更高的計算資源
• 測序錯誤會導致組裝錯誤姚炕，會明顯影響contig的長度
• 短讀難以區(qū)分基因組的重復序列
• 測序覆蓋度不均一摊欠，會影響統(tǒng)計檢驗和結果結果診斷

上述的問題可以嘗試從如下角度進行解決

• 短讀長度：可以通過提供更多樣本丢烘，并且建庫時保證位置足夠隨機
• 數(shù)據(jù)集大小: 使用K-mers算法對數(shù)據(jù)進行組裝柱宦。assembler不再搜尋overlap些椒，而是搜索具有相同k-mers的reads。但是k-mer算法相比較overlap-based算法掸刊，靈敏度有所欠缺免糕，容易丟失一些true overlaps。關鍵在于定義K忧侧。注: K-mer表示一條序列中長度為k的連續(xù)子序列,如ABC的2-mer為AB,BC
• 測序錯誤: 必須保證測序結果足夠正確, 如提高質量控制的標準
• 基因組重復區(qū)： 測序深度要高石窑，結果要正確。如果repeat短于read長度蚓炬，只要保證有足夠多且特異的read松逊。如果repeat長于read，就需要paired ends or “mate-pairs”
• 覆蓋度不均一： 提高深度肯夏，保證隨機
• 組裝結果比較：contig N50, scaffold N50, BUSCO

什么叫做N50

二代數(shù)據(jù)組裝的算法和工具

基因組組裝的組裝工具主要分為三類：基于貪婪算法的拼接方法经宏，基于讀序之間的重疊序列(overlapped sequence)進行拼接的OLC(Overlap-Layout-Consensus)拼接方法和基于德布魯因圖(de bruijn graph)的方法，這三種方法或多或少基于圖論驯击。第一種是最早期的方法烁兰，目前已被淘汰，第二種適用于一代測序產(chǎn)生長片段序列徊都，可以稱之為字符串圖(string graph),第三種是目前二代測序組裝基因組的工具的核心基礎沪斟，也就是要繼續(xù)介紹的de bruijn圖。

示意圖

de bruijn圖由兩部分組成暇矫，節(jié)點(Nodes)和邊(Edges)主之，節(jié)點由k-mers組成，節(jié)點之間要想形成邊就需要是兩個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。

上面是理想情況冀续，實際序列中的測序錯誤琼讽，序列之間的SNP以及基因組低復雜度(重復序列)就會出現(xiàn)如下de brujin圖

測序錯誤引起分支

SNPs的氣泡結構

重復區(qū)域引起的氣泡

用圖的方式表示就是下面情況

幾種比較復雜的圖

組裝軟件的任務就是從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這篇文章有比較好的總結钻蹬，

• ALLPaths-LG是公認比較優(yōu)秀的組裝工具，但消耗內(nèi)存大凭需，并且要提供至少兩個不同大小文庫的數(shù)據(jù)
• SPAdes是小基因組(<100Mb)組裝時的首選
• SOAPdenovo是目前使用率最高的工具(華大組裝了大量的動植物基因組)问欠，效率也挺好肝匆，就是錯誤率也高
• Minia是內(nèi)存資源最省的工具，組裝人類基因組contig居然只要5.7G的RAM顺献，運行23小時旗国，簡直難以相信。

當然工具之間的差別并沒有想象的那么大注整，也沒有想象中那么小能曾，可能在物種A表現(xiàn)一般的工具可能在物種B里就非常好用，因此要多用幾個工具肿轨，選擇其中最好的結果寿冕。

數(shù)據(jù)準備

這里使用來自于GAGE的金黃色葡萄球菌 Staphylococcus aureusa 數(shù)據(jù)進行練習。一方面數(shù)據(jù)量小椒袍，服務器能承受并且跑得快驼唱，另一方面本身基因組就組裝的不錯，等于是考完試能夠自己對答案驹暑。

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測序結果對基因組進行評估，了解基因組的大小岗钩，重復序列含量和復雜度纽窟。基于這些信息兼吓，確定后續(xù)策略以及是否真的需要對該物種進行測序臂港。

基因組survery的核心就是使用k-mers對整體進行評估，k-mers時基因組里長度為k的子序列视搏，當k=17時审孽，ATCG的組合數(shù)就有170億種，也就說理想條件下基因組大小只有超過17Gb才會出現(xiàn)2條一摸一樣的k-mers浑娜。比如說有一個長度為14的序列佑力，給定k-mers的k為8，于是能產(chǎn)生7條長度為8的子序列筋遭，于是推測基因組大小為7bp打颤，但是這似乎和實際的14bp偏離有點遠.

GATCCTACTGATGC (L=14), k-mers for k=8
n = (L-k) + 1 = 14 - 8 + = 7
GATCCTAC,     ATCCTACT,     TCCTACTG,     CCTACTGA,     CTACTGAT,     TACTGATG,     ACTGATGC

如果基因組大小為1MB, 那么當k-mers的k=18時，會得到(1000000-18)+1=999983個不同的k-mers漓滔，與實際大小偏差僅僅只有0.0017%编饺，也就說基因組越大，預測就越接近响驴。這是對單條基因組的估計結果透且，實際上高通量測序會得到基因組30X到50X深度的測序結果，比如說10個拷貝(C)的“GATCCTACTGATGC”在k-mers=8時會有70條子序列豁鲤，

n = [(L-K) + 1] * C = [(14-8)+1]*10 = 70

為了得到實際的基因組大小秽誊，既需要將70除以拷貝數(shù)10鲸沮，那么就得到了和之前一樣的預測值7。當然上述都是理想條件锅论，實際上測序不均一讼溺，低復雜區(qū)域，重復序列等都會影響預測結果棍厌。舉個例子肾胯，"Genome sequencing reveals insights into physiology and longevity of the naked mole rat"的k=17, k_num=52,143,337,243竖席，測序程度可以通過k-mers深度分布曲線來估計

17-k-mer預測基因組大小

圖中耘纱，深度為1的k-mers所占比例最高，表示絕大多數(shù)的k-mers僅僅出現(xiàn)了幾次毕荐，這可能是測序錯誤造成束析。后續(xù)在depth=20逐漸形成一個峰，說明測序測度大概是20x附近憎亚，實際上是19x有極大值员寇。于是基因組的大小就是"52,143,337,243/19=2744386170", 差不多就是2.74Gb

k-mers一般選擇17即可，對于高度重復基因組或者基因組過大第美，可以選擇19甚至31也行蝶锋。但不是越大越好，因為如果一條reads里有一個錯誤位點什往，越大的k-mers就會導致包含這個錯誤位點的k-mers個數(shù)增多扳缕。

根據(jù)上述的介紹，便可以使用jellyfish統(tǒng)計k-mer别威，然后用R作圖對基因組進行評估躯舔。當然這類工具其實已經(jīng)有人開發(fā)，比如說ALLPATHS-LG/FindErrors省古，它不但能夠修正低質量的短讀粥庄，還能初步評估基因組，還有GCE(genome characteristics Estimation)豺妓，由華大基因開發(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運行結束后會有粗略估計基因組大小训堆，粗略估計為4.22Mb。注意臀栈，Kmer_individual_num 數(shù)據(jù)用于gce的輸入?yún)?shù)蔫慧。

隨后用gce程序基于前面的輸出結果進行估計

~/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為主峰對應depth
# -g使用的就是Kmer_individual_num對應值
# -m 選擇估算模型，真實數(shù)據(jù)選擇1权薯，表示連續(xù)型

在這次的日志文件中有預測后的結果4.34Mb姑躲，但是根據(jù)NCBI的數(shù)據(jù)睡扬，這個物種的基因組大小是2.8M左右。因此使用k-mers通過數(shù)學方法預測存在一定的局限性黍析，需要結合流式細胞儀和粗組裝的結果卖怜。

雖然也可以使用FindErrors對基因組進行評估，但是我實際使用時出現(xiàn)了各種問題阐枣，這里不做介紹马靠。其他的工具也是大同小異，不做額外推薦蔼两。

基因組正式組裝

當你拿到測序數(shù)據(jù)后甩鳄，就可以按照如下幾步處理數(shù)據(jù)。第一步是數(shù)據(jù)質控控制额划，這一步對于組裝而言非常重要妙啃，處理前和處理后的組裝結果可能會天差地別；第二步俊戳，根據(jù)經(jīng)驗確定起始參數(shù)揖赴，如K-mer和覆蓋率；第三步抑胎，使用不同軟件進行組裝燥滑；第四步，評估組裝結果阿逃，如contig N50, scaffold N50, 判斷是否需要修改參數(shù)重新組裝铭拧。

原始數(shù)據(jù)質量控制

盡管目前的測序技術已經(jīng)非常成熟，公司提供的數(shù)據(jù)一般都可以直接用于普通的項目(特殊項目如miRNA-seq除外)盆昙。但由于 de novo 組裝對數(shù)據(jù)質量比較敏感羽历，因此需要通過質控來降低偏差。原始數(shù)據(jù)質量控制分為四個部分：

• 了解數(shù)據(jù)質量: 了解質量這一步可以暫時忽略淡喜，基本上基因組測序的結果都能通過FastQC的標準秕磷。
• 去接頭和低質量reads過濾: 去接頭和低質量reads過濾可供選擇的軟件非常之多，如NGSQCToolkit, Trimmomatic, cutadapter, 似乎都是國外開發(fā)的軟件炼团，但其實國內(nèi)也有一款很優(yōu)秀的工具叫做fastp
• 去除PCR重復: 去重一般都是在比對后根據(jù)位置信息進行澎嚣，沒有基因組的話只能根據(jù)PE的reads是否完全一樣進行過濾。從理論上說瘟芝，測序相當于是從基因組上隨機抽樣易桃，不太可能存在完全一摸一樣的兩條序列。不過貌似只有FastUniq能做這件事情锌俱，后來有一個人寫了sequniq晤郑。
• reads修正: 除了過濾或修剪低質量的reads外，一般而言，還需要對reads中的錯誤堿基進行修正造寝。尤其當測序的覆蓋度比較高時磕洪，錯誤的reads也就越來越多，會對 de novo 組裝造成不良的影響诫龙。工具有BLESS2, BFC, Musket等析显，其中BLESS2的效率最高，效果也不錯签赃。

去接頭和低質量reads過濾這一步推薦fastp谷异，主要是因為它基于C/C++，運行速度塊锦聊。

# 使用, 項目文件夾下
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的測序結果以及看了FastQC的結果才打消了我的疑慮伶选，沒錯型雳，以前的數(shù)據(jù)質量就是那么差擎鸠。注，除非是去接頭，否則不建議通過刪除序列的方式提高質量哲嘲。

質控另一個策略是對短讀中一些可能的錯誤堿基進行糾正，測序錯誤會引入大量無意義的K-mers高氮，從而增加運算復雜度彼哼。此處使用BFC對測序質量：

~/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

處理前后

總之，質控的目標是在不引入的錯誤的情況下盡量提高整體質量河绽，這一步對后續(xù)的組裝影響很大己单，所以盡量做這一步，除非組裝軟件要求你別做耙饰，那你就不要手賤了纹笼。

使用不同工具和參數(shù)進行組裝

二代組裝可供選擇的工具很多, 但是主流其實就那么幾個, 所以組裝的時候選擇3~5個工具運行比較結果即可，比如說MaSuRCA
, IDBA-UD, SOAPdenovo2, Abyss, velvet和Spades苟跪。當然一旦你選擇一個軟件準備運行的時候廷痘，你就會遇到參數(shù)選擇問題，比如怎么確定k-mers件已，組裝軟件最基礎也是最核心的參數(shù)笋额。這里有幾條原則值得借鑒：

• k要大于log4(基因組大小)，如果數(shù)學不好篷扩，無腦選擇20以上
• 盡量減少測序錯誤形成的k-mers, 因為這是無意義的噪音, 也就是要求k不能過大
• 當然k也不能太小兄猩，否則會導致重復壓縮,比如說ATATATA，在2kmers的情況下，就只有AT了
• 測序深度越高枢冤，K值也就可以選擇的越大

但是說了那么多援岩，你依舊不知道應該選擇什么樣的K，如果你的計算資源無限掏导，那么窮舉法最簡單粗暴享怀。如果窮舉法不行，那么建議先用k=21, 55,77 組裝一下contig, 對不同參數(shù)的contig N50有一個大致的了解趟咆，然后繼續(xù)調(diào)整添瓷。此外還有一個工具叫做KmerGenie可以預測一個初始值≈瞪矗總之鳞贷，讓我們先運行第一個工具--SPAdes，可通過bioconda安裝虐唠。

SPAdes全稱是圣彼得堡基因組組裝工具搀愧，包含了一系列組裝工具處理不同的項目，如高雜合度的dipSPAdes疆偿，宏基因組的metaSPAdes咱筛。官方文檔中以大腸桿菌為例運行整個流程，花了將近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運行，從而挑選比較合適的k值撤蟆，當然你還可以自己設置奕塑，比如說-k 21,55,77。最后結果為糾正后的短讀數(shù)據(jù)家肯，組裝后的contig, 組裝后的scaffold, 不同格式的組裝graph龄砰。

同樣運行多k-mers運行后比較的工具還有IDBA，它也有一系列的工具息楔。IDBA是基礎版寝贡，IDBA-UD適用于宏基因組和單細胞測序的數(shù)據(jù)組裝，IDBA-Hybrid則是基于相似的基因組提高組裝結果值依，IDBA-Tran是專門處理轉錄組數(shù)據(jù)圃泡。對于無參考基因組組裝，作者推薦使用IDBA-UD愿险。

IDBA-UD工具要求將兩個配對的短讀文件合并成一個颇蜡，我們的原始數(shù)據(jù)需要先用它提供的fq2fa先轉換格式

# 項目文件夾下
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相關參數(shù)為--mink,--maxk,--step, 通過--read_level_x 傳入不同大小的文庫价说，也提供了短讀糾正的相關參數(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

運行結束后在assembly/idba_ud下會生成一系列的文件，其中結果文件為contig.fa和scaffold.fa风秤。

最后介紹一個要手動運行不同k-mers的工具鳖目，如ABySS, 它有一個亮點，就是能夠可以使用多個計算節(jié)點缤弦。我們使用k=31進行組裝

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結尾，其次第二個文庫才37bp, 所以比對軟件要選擇bowtie,否則你運行一定會遇到histogram xxx.hist is empty的報錯碍沐。當然到最后狸捅，這個問題我都沒有解決掉，所以我放棄了累提。

雖然看起來abyss用起來很簡單尘喝，但其實背后的工作流程還是比較復雜，如下是它的流程示意圖

flowchart

小結一下斋陪，這里用到了spades, idba,abyss三種工具對同一種物種進行組裝朽褪，得到對應的contig結果，重點在于k-mers的選擇无虚。contig是組裝的第一步缔赠，也是非常重要的一步，為了保證后續(xù)搭scaffold和基因組補洞等工作的順利骑科，我們先得挑選一個比較高質量的contig橡淑。

組裝可視化和評估

理想條件下，我們希望一個物種有多少染色體咆爽，結果最好就只有多少個contig。當然對于二代測序而言置森，這絕對屬于妄想斗埂，可以通過一款graph可視化工具bandage來感受一下最初得到的contig graph是多么復雜。

Bandage

一般看這圖直觀感受就是怎么那么多節(jié)凫海，這些節(jié)就是造成contig不連續(xù)的元兇呛凶。不同組裝工具在構建de bruijn graph的差異不會那么大，contig的數(shù)量和大小和不同工具如何處理復雜節(jié)點有關行贪。我們希望得到的contig文件中漾稀，每個contig都能足夠的長，能夠有一個完整的基因結構建瘫，歸納一下就是3C原則:

• 連續(xù)性(Contiguity): 得到的contig要足夠的長
• 正確性(Correctness): 組裝的contig錯誤率要低
• 完整性(Completeness)：盡可能包含整個原始序列

但是這三條原則其實是相互矛盾的崭捍，連續(xù)性越高，就意味著要處理更多的模糊節(jié)點啰脚，會導致整體錯誤率上升殷蛇，為了保證完全的正確，那么就會導致contig非常的零碎。此外粒梦，這三條原則也比較定性亮航，我們需要更加定量的數(shù)值衡量，比如說contig數(shù), 組裝的總長度等, N50等匀们。問題來了缴淋，什么叫做N50呢，

N50定義比較繞口泄朴，有一種只可意會不可言傳的感覺宴猾，所以索性看圖

N50和NG50

假設一個基因組的大小為10，但是這個值只有神知道叼旋，你得到的信息就是組裝后有3個contig,長度分別為"3,4,1,1"仇哆，所以組裝總長度為9。為了計算N50夫植，我們需要先把contig從大到小排列讹剔，也就是"4,3,1"。然后先看最大的contig详民，長度是4延欠，他的長度是不是超過組裝總大小的一半了嗎？如果是沈跨，那么N50=4, 4 < 4.5, 不是由捎。 那么在此基礎上加上第二長的contig,也就是4+3=7, 是不是超過一半了？7>4.5, 那么N50=3. 因此饿凛，N50的定義可以表述為"使得累加后長度超過組裝總長度一半的contig的長度就是N50"狞玛。為了方便管理和使用軟件，建議建立如下幾個文件夾

N50是基于一個未知的基因組得到得結果涧窒，如果基因組測序比較完整心肪，那么就可以計算NG50，也就是"使得累加后長度超過基因組總長度一半的contig的長度就是NG50"纠吴。NA50比較稍微復雜硬鞍，需要將組裝結果進一步比對到參考基因組上，以contig實際和基因組匹配的長度進行排序計算戴已。

說完N50固该，我們介紹兩款工具，QUAST和BUSCO糖儡。

QUAST使用質量標準(quality metrics)來評估不同組裝工具和不同參數(shù)的組裝效果伐坏，無論是否有基因組都可以使用。我們分別以有參和無參兩種模式比較Minia,IDBA和SPAdes三個組裝的運行結果

# 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

quast運行結果

這個結果非常直觀的告訴我們一個事實就是spades組裝的contigs`各方面表現(xiàn)都很優(yōu)秀休玩，minia由于內(nèi)存使用率最低著淆，所以組裝效果一般也是可以理解劫狠。

BUSCO通過同源基因數(shù)據(jù)庫從基因完整度來評價基因組組裝結果。BUSCO首先構建了不同物種的最小基因集永部，然后使用HMMER独泞，BLAST,Augustus等工具分析組裝結果中的同源基因，從而定量評估組裝是否完整苔埋。

busco -i assembly/spades/contigs.fasta -o result -l /home/wangjw/db/busco/bacteria_odb9 -m genome -f

運行結果會在當前目錄下的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則是重復數(shù)荞膘。正所謂沒有比較，就沒有傷害玉工，我們拿之前QUAST對比中表現(xiàn)比較差的minia結果作為對比羽资。

    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不但組裝效果好屠升，而且基因完整度也高，當然它的內(nèi)存消耗也是很嚴重狭郑。這都是取舍的過程腹暖。

附錄

參考資料

• Bandage: https://github.com/rrwick/Bandage/wiki
• QUAST: http://quast.bioinf.spbau.ru/manual.html

軟件安裝

由于不同軟件對不同的基因組的適合度不同，一般都需要參數(shù)多個工具的不同參數(shù)翰萨，根據(jù)N50和BUSCO等衡量標準選擇比較好的結果脏答。為了避免后續(xù)花篇幅在工具安裝上，因此先準備后續(xù)的分析環(huán)境亩鬼。對于組裝而言殖告，我們需要安裝如下工具:

• 質量控制:
  • FastQC
  • fastp
  • BFC
• 主流組裝工具:
  • ABySS
  • IDBA
  • SOAPdenovo2
  • Velvet
  • Sapdes
  • Minia
  • Ray
  • MasuRCA
• 基因組組裝評價工具
  • BUSCO
  • Quast
• 基因結構預測和功能注釋暫時不在考慮范圍內(nèi)

更多相關工具見https://biosphere.france-bioinformatique.fr/wikia2/index.php/Tools_directory_in_Assembly_and_Annotation_(Lexicographic_ordering)

以下操作所用服務器的基本信息為：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。當然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ù)我之前關于GCC編譯的文章袒餐，程序編譯不成功大多是因為找不到頭文件(存放在include目錄下)和鏈接庫文件(存放在lib目錄下)飞蛹，默認編譯頭文件只會搜索/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: 運行時搜索鏈接文件的路徑

即添加如下幾行內(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進行組裝，從而確定最優(yōu)的K库物。

cd ~/src
git clone https://github.com/loneknightpy/idba.git
idba/build.sh
mv idba ~/opt/biosoft/

組裝軟件四：MaSuRCA霸旗，能夠純用二代，也能二代三代測序混合使用戚揭，先用 de bruijn 圖構建長reads诱告，然后再用OLC算法進行組裝

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

質控軟件一: 原本是要推薦BLESS2,但是這個軟件在編譯完成后出現(xià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

質控軟件二: fastp是一款基于C/C++編寫的工具精居，速度會比較塊，而且運行之后會有比較好看的圖哦

mkdir -p ~/opt/biosoft/fastp
cd ~/opt/biosoft/fastp
wget http://opengene.org/fastp/fastp
chmod a+x ./fastp

評估工具一：Quast, 它通過比較N50,N G50等參數(shù)來評價基因組組裝質量.Quast由Python編寫潜必，推薦使用bioconda安裝

conda create --name assembly python=2.7
source activate assembly
conda install quast

評估工具二: BUSCO靴姿，這是一個利用進化信息從基因完整性角度評估組裝準確性的工具，推薦使用biconda安裝磁滚。

source activate assembly
conda install busco

盡管conda安裝了busco佛吓，但是離實際運行還需要添加幾個環(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/選擇進化上接近的評估數(shù)據(jù)集，比如說你如果研究魚晒他，那么"actinopterygii(輻鰭魚類)"就比"metazoa(多細胞動物)"更加合適.

實際運行時可能還存在鏈接庫無法找尋以至于程序出錯吱型，解決方法就是將相對應或著接近的庫拷貝或軟鏈接到~/miniconda3/env/assembly/lib下。