基因組組裝

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

目前基於二代測序的組裝存在挑戰：

• 全基因組測序得到的短讀遠遠小於原來的分子長度
• 高通量測序得到海量資料會增加組裝的計算複雜性，消耗更高的計算資源
• 測序錯誤會導致組裝錯誤，會明顯影響contig的長度
• 短讀難以區分基因組的重複序列
• 測序覆蓋度不均一，會影響統計檢驗和結果結果診斷

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

• 短讀長度：可以通過提供更多樣本，並且建庫時保證位置足夠隨機
• 資料集大小: 使用K-mers演算法對資料進行組裝。assembler不再搜尋overlap，而是搜尋具有相同k-mers的reads。但是k-mer演算法相比較overlap-based演算法，靈敏度有所欠缺，容易丟失一些true overlaps。關鍵在於定義K。注: K-mer表示一條序列中長度為k的連續子序列,如ABC的2-mer為AB,BC
• 測序錯誤: 必須保證測序結果足夠正確, 如提高質量控制的標準
• 基因組重複區： 測序深度要高，結果要正確。如果repeat短於read長度，只要保證有足夠多且特異的read。如果repeat長於read，就需要paired ends or “mate-pairs”
• 覆蓋度不均一： 提高深度，保證隨機
• 組裝結果比較：contig N50, scaffold N50, BUSCO

二代資料組裝的演算法和工具

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

de bruijn圖由兩部分組成，節點(Nodes)和邊(Edges)，節點由k-mers組成，節點之間要想形成邊就需要是兩個k-mers存在K-1個完全匹配。比如說，ACTG, CTGC, TGCC在K=3時的k-mers為ACT,CTG,TGC,GCC，可以表示為ACT -> CTG -> TGC -> GC.

對於de brujin圖而言，冗餘序列不會影響k-mers的數量，比如說ACTG,ACTG,CTGC,CTGC,CTGC,TGCC,TGCC在K=3時依舊錶示為ACT -> CTG -> TGC -> GCC。

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

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

組裝軟體的任務就是從k-mers形成的圖按照一定的演算法組裝出可能的序列，根據"GAGE: A critical evaluation of genome assemblies and assembly algorithms"以及自己的經驗，目前二代資料比較常用的工具有Velvet, ABySS, AllPaths/AllPaths-LG, Discovar, SOAPdenovo, Minia, spades,Genomic Assemblers這篇文章有比較好的總結，

• ALLPaths-LG是公認比較優秀的組裝工具，但消耗記憶體大，並且要提供至少兩個不同大小文庫的資料
• SPAdes是小基因組(<100Mb)組裝時的首選
• SOAPdenovo是目前使用率最高的工具(華大組裝了大量的動植物基因組)，效率也挺好，就是錯誤率也高
• Minia是記憶體資源最省的工具，組裝人類基因組contig居然只要5.7G的RAM，執行23小時，簡直難以相信。

當然工具之間的差別並沒有想象的那麼大，也沒有想象中那麼小，可能在物種A表現一般的工具可能在物種B裡就非常好用，因此要多用幾個工具，選擇其中最好的結果。

資料準備

這裡使用來自於GAGE的金黃色葡萄球菌 Staphylococcus aureusa 資料進行練習。一方面資料量小，伺服器能承受並且跑得快，另一方面本身基因組就組裝的不錯，等於是考完試能夠自己對答案。

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

在正式組裝之前，需要先根據50X左右的illumina測序結果對基因組進行評估，瞭解基因組的大小，重複序列含量和複雜度。基於這些資訊，確定後續策略以及是否真的需要對該物種進行測序。

基因組survery的核心就是使用k-mers對整體進行評估，k-mers時基因組裡長度為k的子序列，當k=17時，ATCG的組合數就有170億種，也就說理想條件下基因組大小隻有超過17Gb才會出現2條一摸一樣的k-mers。比如說有一個長度為14的序列，給定k-mers的k為8，於是能產生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除以拷貝數10，那麼就得到了和之前一樣的預測值7。當然上述都是理想條件，實際上測序不均一，低複雜區域，重複序列等都會影響預測結果。舉個例子，"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所佔比例最高，表示絕大多數的k-mers僅僅出現了幾次，這可能是測序錯誤造成。後續在depth=20逐漸形成一個峰，說明測序測度大概是20x附近，實際上是19x有極大值。於是基因組的大小就是"52,143,337,243/19=2744386170", 差不多就是2.74Gb

k-mers一般選擇17即可，對於高度重複基因組或者基因組過大，可以選擇19甚至31也行。但不是越大越好，因為如果一條reads裡有一個錯誤位點，越大的k-mers就會導致包含這個錯誤位點的k-mers個數增多。

根據上述的介紹，便可以使用jellyfish統計k-mer，然後用R作圖對基因組進行評估。當然這類工具其實已經有人開發，比如說ALLPATHS-LG/FindErrors，它不但能夠修正低質量的短讀，還能初步評估基因組，還有GCE(genome characteristics Estimation)，由華大基因開發出來的一款基因組評估工具等。為了避免重複造輪子，簡單就用這些工具即可。

使用GCE評估基因組: 先用kmer_freq_hash統計k-mer頻數

# Staphylococcus_aureus專案根目錄下
mkdir genome_survey && cd genome_survey
## 提供用於read的位置資訊
ls raw-data/lib1/frag_*.fastq.gz > genome_survey/reads.list
## k-mer_freq_hash統計
~/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 資料用於gce的輸入引數。

隨後用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 選擇估算模型，真實資料選擇1，表示連續型

在這次的日誌檔案中有預測後的結果4.34Mb，但是根據NCBI的資料，這個物種的基因組大小是2.8M左右。因此使用k-mers通過數學方法預測存在一定的侷限性，需要結合流式細胞儀和粗組裝的結果。

雖然也可以使用FindErrors對基因組進行評估，但是我實際使用時出現了各種問題，這裡不做介紹。其他的工具也是大同小異，不做額外推薦。

基因組正式組裝

當你拿到測序資料後，就可以按照如下幾步處理資料。第一步是資料質控控制，這一步對於組裝而言非常重要，處理前和處理後的組裝結果可能會天差地別；第二步，根據經驗確定起始引數，如K-mer和覆蓋率；第三步，使用不同軟體進行組裝；第四步，評估組裝結果，如contig N50, scaffold N50, 判斷是否需要修改引數重新組裝。

原始資料質量控制

儘管目前的測序技術已經非常成熟，公司提供的資料一般都可以直接用於普通的專案(特殊專案如miRNA-seq除外)。但由於 de novo 組裝對資料質量比較敏感，因此需要通過質控來降低偏差。原始資料質量控制分為四個部分：

• 瞭解資料質量: 瞭解質量這一步可以暫時忽略，基本上基因組測序的結果都能通過FastQC的標準。
• 去接頭和低質量reads過濾: 去接頭和低質量reads過濾可供選擇的軟體非常之多，如NGSQCToolkit, Trimmomatic, cutadapter, 似乎都是國外開發的軟體，但其實國內也有一款很優秀的工具叫做fastp
• 去除PCR重複: 去重一般都是在比對後根據位置資訊進行，沒有基因組的話只能根據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，一度讓我懷疑軟體是否出現了問題，直到我用同樣的程式碼處理現在Illumina的測序結果以及看了FastQC的結果才打消了我的疑慮，沒錯，以前的資料質量就是那麼差。注，除非是去接頭，否則不建議通過刪除序列的方式提高質量。

質控另一個策略是對短讀中一些可能的錯誤鹼基進行糾正，測序錯誤會引入大量無意義的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

總之，質控的目標是在不引入的錯誤的情況下儘量提高整體質量，這一步對後續的組裝影響很大，所以儘量做這一步，除非組裝軟體要求你別做，那你就不要手賤了。

使用不同工具和引數進行組裝

二代組裝可供選擇的工具很多, 但是主流其實就那麼幾個, 所以組裝的時候選擇3~5個工具執行比較結果即可，比如說MaSuRCA
, IDBA-UD, SOAPdenovo2, Abyss, velvet和Spades。當然一旦你選擇一個軟體準備執行的時候，你就會遇到引數選擇問題，比如怎麼確定k-mers，組裝軟體最基礎也是最核心的引數。這裡有幾條原則值得借鑑：

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

但是說了那麼多，你依舊不知道應該選擇什麼樣的K，如果你的計算資源無限，那麼窮舉法最簡單粗暴。如果窮舉法不行，那麼建議先用k=21, 55,77 組裝一下contig, 對不同引數的contig N50有一個大致的瞭解，然後繼續調整。此外還有一個工具叫做KmerGenie可以預測一個初始值。總之，讓我們先執行第一個工具--SPAdes，可通過bioconda安裝。

SPAdes全稱是聖彼得堡基因組組裝工具，包含了一系列組裝工具處理不同的專案，如高雜合度的dipSPAdes，巨集基因組的metaSPAdes。官方文件中以大腸桿菌為例執行整個流程，花了將近1個小時。我們的資料集比較小，速度會更快

# 專案根資料夾下
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/

你會發現之前說的k-mers在這裡根本沒出現，而且用的也是原始資料，這是因為spades.py有一個元件BayesHammer處理測序錯誤，並且它是多K類組裝工具(multi-k assembly), 也就是說它會自動選擇不同的K執行，從而挑選比較合適的k值，當然你還可以自己設定，比如說-k 21,55,77。最後結果為糾正後的短讀資料，組裝後的contig, 組裝後的scaffold, 不同格式的組裝graph。

同樣執行多k-mers執行後比較的工具還有IDBA，它也有一系列的工具。IDBA是基礎版，IDBA-UD適用於巨集基因組和單細胞測序的資料組裝，IDBA-Hybrid則是基於相似的基因組提高組裝結果，IDBA-Tran是專門處理轉錄組資料。對於無參考基因組組裝，作者推薦使用IDBA-UD。

IDBA-UD工具要求將兩個配對的短讀檔案合併成一個，我們的原始資料需要先用它提供的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相關引數為--mink,--maxk,--step, 通過--read_level_x 傳入不同大小的文庫，也提供了短讀糾正的相關引數--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, 它有一個亮點，就是能夠可以使用多個計算節點。我們使用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用起來很簡單，但其實背後的工作流程還是比較複雜，如下是它的流程示意圖

小結一下，這裡用到了spades, idba,abyss三種工具對同一種物種進行組裝，得到對應的contig結果，重點在於k-mers的選擇。contig是組裝的第一步，也是非常重要的一步，為了保證後續搭scaffold和基因組補洞等工作的順利，我們先得挑選一個比較高質量的contig。

組裝視覺化和評估

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

一般看這圖直觀感受就是怎麼那麼多節，這些節就是造成contig不連續的元凶。不同組裝工具在構建de bruijn graph的差異不會那麼大，contig的數量和大小和不同工具如何處理複雜節點有關。我們希望得到的contig檔案中，每個contig都能足夠的長，能夠有一個完整的基因結構，歸納一下就是3C原則:

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

但是這三條原則其實是相互矛盾的，連續性越高，就意味著要處理更多的模糊節點，會導致整體錯誤率上升，為了保證完全的正確，那麼就會導致contig非常的零碎。此外，這三條原則也比較定性，我們需要更加定量的數值衡量，比如說contig數, 組裝的總長度等, N50等。問題來了，什麼叫做N50呢，

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

假設一個基因組的大小為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)來評估不同組裝工具和不同引數的組裝效果，無論是否有基因組都可以使用。我們分別以有參和無參兩種模式比較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

這個結果非常直觀的告訴我們一個事實就是spades組裝的contigs`各方面表現都很優秀，minia由於記憶體使用率最低，所以組裝效果一般也是可以理解。

BUSCO通過同源基因資料庫從基因完整度來評價基因組組裝結果。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內容如下

# 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值表示可能缺少的基因數，D則是重複數。正所謂沒有比較，就沒有傷害，我們拿之前QUAST對比中表現比較差的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不但組裝效果好，而且基因完整度也高，當然它的記憶體消耗也是很嚴重。這都是取捨的過程。

附錄

參考資料

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

軟體安裝

由於不同軟體對不同的基因組的適合度不同，一般都需要引數多個工具的不同引數，根據N50和BUSCO等衡量標準選擇比較好的結果。為了避免後續花篇幅在工具安裝上，因此先準備後續的分析環境。對於組裝而言，我們需要安裝如下工具:

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

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

以下操作所用伺服器的基本資訊為：Linux的核心為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

系統自帶的GCC版本是4.8，而BLESS2要求4.9+, ABySS要求6.0+，直接編譯這些工具可能會出錯，但直接升級系統的GCC版本可能會影響整體穩定性，因此推薦將在opt/sysoft下安裝高版本的GCC。當然GCC的版本也不是越高越好，最好和作者開發的版本一致，也就是他們要求的最低版本。

# 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

根據我之前關於GCC編譯的文章，程式編譯不成功大多是因為找不到標頭檔案(存放在include目錄下)和連結庫檔案(存放在lib目錄下)，預設編譯標頭檔案只會搜尋/usr/include,/usr/local/include, 而連結庫檔案只會搜尋/lib,/usr/lib[64],/usr/local/lib[64]. 為了讓編譯完成的GCC的標頭檔案和連結庫檔案能被搜尋到，需要在~/.bashrc檔案中新增幾個環境變數:

• PKG_CONFIG_PATH: 同時新增搜尋標頭檔案和連結標頭檔案的路徑
• C_INCLUDE_PATH: 編譯時搜尋標頭檔案的路徑
• LIBRARY_PATH: 編譯時搜尋連結檔案的路徑
• LD_LIBRARY_PATH: 執行時搜尋連結檔案的路徑

即新增如下幾行內容到~/.bashrc檔案中，並執行source ~/.bashrc更新環境變數。

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。記憶體不太夠的話可以嘗試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中新增環境變數或手動修改也行
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進行組裝，從而確定最優的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,但是這個軟體在編譯完成後出現各種核心轉移的毛病，和我的系統相性太差，於是改用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等引數來評價基因組組裝質量.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，但是離實際執行還需要新增幾個環境變數和不同物種的基因資料集，請使用printenv確保如下如下幾個路徑都已經新增到環境變數中。

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/選擇進化上接近的評估資料集，比如說你如果研究魚，那麼"actinopterygii(輻鰭魚類)"就比"metazoa(多細胞動物)"更加合適.

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