關(guān)于更多生物醫(yī)療大數(shù)據(jù)分析工具和軟件的介紹和使用請看六點(diǎn)了官網(wǎng)[1]。

目錄

1、基因組組裝

2、基于De-Bruijn Graph的組裝算法

3、SOAPdenovo的安裝和使用說明：安裝、說明、配置、運(yùn)行

4、SOAPdenovo案例實(shí)戰(zhàn)：數(shù)據(jù)下載、配置、運(yùn)行、輸出

寫在前面

大家好，這是我們六點(diǎn)了給大家介紹生物信息大數(shù)據(jù)分析 基因組數(shù)據(jù)分析系列文章第一篇。我們會持續(xù)為大家分享關(guān)于生物醫(yī)療大數(shù)據(jù)處理相關(guān)的知識和案例，希望幫助大家更好地進(jìn)行自己項(xiàng)目中生物醫(yī)療健康大數(shù)據(jù)處理工作。本篇文章主要四部分來為大家介紹基因組的denovo的知識和以及詳細(xì)應(yīng)用案例。①基因組組裝、 ②基于De-Bruijn Graph的組裝算法、 ③SOAPdenovo的安裝和使用說明：安裝、說明、配置、運(yùn)行，以及 ④SOAPdenovo案例實(shí)戰(zhàn)：數(shù)據(jù)下載、配置、運(yùn)行、輸出。

1. 基因組組裝

基因組組裝 (Genome assembly)是生物信息學(xué)領(lǐng)域的核心問題, 想要深入研究一個生物體，獲得參考基因組是第一步也是必須的一步。基因組組裝是將原始的下機(jī)序列還原成DNA序列片段、以至于整個物種全基因組序列的過程。

基因組組裝是基因組分析的關(guān)鍵，對物種起源與進(jìn)化，挖掘功能基因進(jìn)而研究疾病發(fā)生和發(fā)展具有重大意義。

然而由于目前市面上廣為應(yīng)用的二代測序技術(shù)獲得的測序序列一般都較短，因此如何通過短片段組裝成完整的基因組成了亟待解決的問題。

基因組組裝可分為基于參考基因組的組裝（Mapping assembly）和從頭組裝（denovo assembly）。兩者主要的區(qū)別在于是否存在已知的基因組參考序列作為參照。本文我們主要介紹的是denovo組裝，即不依賴任何基因組參考序列相關(guān)信息而進(jìn)行的序列組裝。目前，應(yīng)用于主流的基因組denovo組裝的算法主要有兩個[1]：OLC方法 (Overlap-Layout-Consensus)和 DBG方法 (De-Bruijn Graph)[2]。

而DBG方法的核心思想是將序列拼接問題轉(zhuǎn)化為人們所熟知的歐拉圖（Euler Graph）問題[3]。

DBG方法內(nèi)存消耗相對較低，運(yùn)算速度快，且準(zhǔn)確率高。

目前主流的基因組裝算法都是基于DBG方法改進(jìn)設(shè)計(jì)的。

2.基于De-Bruijn Graph的組裝算法

前面我們說到基因組denovo組裝兩種方法，下面主要展開說說基于De-Bruijn Graph的組裝算法的基本原理。此處，就以目前使用比較廣泛，由華大基因團(tuán)隊(duì)開發(fā)的SOAPdenovo[4]為例。軟件的參考文獻(xiàn)[5]有興趣可以在參考資料看一下讀讀。

A：基因組DNA打斷成小的片段，進(jìn)行建庫和雙端測序。150~500bp的進(jìn)行直接雙端測序，長的片段2-10kb的則先進(jìn)行環(huán)化再進(jìn)行雙端測序。

B：組裝的核心部分，進(jìn)行De-Bruijn Graph的構(gòu)建。構(gòu)建De-Bruijn圖的第一步是將測序read k-mer化，而所謂的k-mer是指將reads分成包含k個堿基的字符串，即拿一個k長度的窗口在整個read上1個堿基一個堿基的滑動，每次滑動窗口內(nèi)部都會產(chǎn)生一個k大小的序列,即為一個k-mer，因此一般長短為m的reads可以分成m-k+1個k-mers。其中k一定是奇數(shù)，如果是偶數(shù)遇到回文序列可能會產(chǎn)生完全相同的k-mers。我們將k-mers作為圖的節(jié)點(diǎn)，如果兩個節(jié)點(diǎn)有 K-1個共同重疊子集，就把兩個節(jié)點(diǎn)連接在一起，這樣就會形成De-Bruijn Graph，可以看到該圖可以很好地展現(xiàn)出序列的順序信息。

C：進(jìn)行圖結(jié)構(gòu)的精簡。盡管前面步驟已經(jīng)初步構(gòu)建出圖形，但是實(shí)際上由于測序錯誤，重復(fù)，雜合等原因，圖上會出現(xiàn)很多類似翼尖（tips）、氣泡（bubbles）等問題，因此還需要進(jìn)一步簡化。此處簡化主要包含四個方面：1）去除tips（可能為測序錯誤導(dǎo)致的）；2）去除低覆蓋度的路徑；3）解開微小重復(fù)的區(qū)域（可以通過read穿過來解決）4）合并bubbles氣泡區(qū)（可能為測序錯誤，重復(fù)或者雜合導(dǎo)致的）。

D: 拆分出contig。在重復(fù)的節(jié)點(diǎn)處剪斷，輸出contigs。

E: 構(gòu)建scaffolds。重新用reads和contigs進(jìn)行比對，使用paired-end信息來把單一的contigs連接成scaffolds。1）paired reads 比對到contigs上，使臨近的contig建立連接；3）paired-end信息的不同插入片段被用來一步步從短到長的建立scaffold.

F: 最終是把多個scaffold組裝成無GAP的基因組序列。

3.SOAPdenovo的安裝和使用說明

3.1安裝

SOAPdenovo目前已更新到SOAPdenovo2,github[6]鏈接:https://github.com/aquaskyline/SOAPdenovo2。

直接下載二進(jìn)制[7]（https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/）

mac鏈接: <https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/bin/r240/SOAPdenovo2-bin-r240-mac.tgz/download>linux鏈接: <https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/bin/r240/SOAPdenovo2-bin-LINUX-generic-r240.tgz/download>

源代碼安裝:

wget <https://github.com/aquaskyline/SOAPdenovo2/archive/refs/tags/r242.zip>unzip r242.zipcd SOAPdenovo2-r242make

安裝完可以看到SOAPdenovo-127mer，SOAPdenovo-63mer兩個執(zhí)行文件。63mer代表支持的kmer最大長度為63，127mer代表支持的kmer最大長度為127，除了支持的kmer長度不同外，其他用法完全相同。

3.2使用說明

SOAPdenovo由于計(jì)算量相對較大，對電腦的配置有一定的要求，官網(wǎng)對運(yùn)行配置的說明：SOAPdenovo 的適用目標(biāo)是大型植物和動物基因組，盡管它也適用于細(xì)菌和真菌基因組。它運(yùn)行在至少 5G 物理內(nèi)存的 64 位 Linux 系統(tǒng)上。對于像人類這樣的大基因組，大約需要 150 GB 的內(nèi)存。運(yùn)行SOAPdenovo-63mer即可看到SOAPdenovo主要包含了以下6個子命令：

1. pregraph? ? ? ? construct kmer-graph # 2. sparse_pregraph construct sparse kmer-graph 3. contig? ? ? ? ? eliminate errors and output contigs 4. map? ? ? ? ? ? map reads to contigs 5. scaff? ? ? ? ? construct scaffolds 6. all? ? ? ? ? ? do pregraph-contig-map-scaff in tur

其中，1-5分別表示組裝的4個步驟（1,2是兩種構(gòu)圖方式，二選一），all則用于一次執(zhí)行以上的4個步驟。實(shí)際應(yīng)用中，可以使用SOAPdenovo all 一步式跑完，也可以分成4步單獨(dú)去跑。

3.3. 軟件運(yùn)行前準(zhǔn)備工作-配置文件

soapdenovo需要一個配置文件config_file,里面給定輸入文件和一些參數(shù)設(shè)置。

下面是配置文件的示例和說明：

max_rd_len=100# 全局配置參數(shù)：如果序列大于該長度，會被切成該長度，然后在分析[LIB]#每個文庫的配置以[LIB]開頭avg_ins=200#文庫插入片段的平均長度，在實(shí)際設(shè)置時(shí)，可以參考文庫size分布圖，取峰值即可reverse_seq=0#是否需要將序列反向互補(bǔ)，對于pair-end數(shù)據(jù)，不需要反向互補(bǔ)，設(shè)置為0；對于mate-pair數(shù)據(jù)，需要反向互補(bǔ)，設(shè)置為1asm_flags=3#1表示只組裝contig. 2表示只組裝scaffold,3表示同時(shí)組裝contig和scaffold,4表示只補(bǔ)gaprd_len_cutoff=100#序列長度閾值，作用和max_rd_len相同，大于該長度的序列會被切除到該長度rank=1#設(shè)置不同文庫數(shù)據(jù)的優(yōu)先級順序，取值范圍為整數(shù)，rank值相同的多個文庫，在組裝scaffold時(shí)，會同時(shí)使用。pair_num_cutoff=3# contig或者scaffold之前的最小overlap個數(shù)，對于pair-end數(shù)據(jù)，默認(rèn)值為3；對于mate-paird數(shù)據(jù)，默認(rèn)值為5map_len=32#比對長度的最小閾值，對于pair-end數(shù)據(jù)，默認(rèn)值為32；對于mate-pair數(shù)據(jù)，默認(rèn)值為35q1=fastq1_read_1.fq# read1的序列文件q2=fastq1_read_2.fq# read2的序列文件

3.4. 軟件運(yùn)行

拆分式：

#step1:SOAPdenovo-63mer pregraph -s config_file -K 63 -R -o graph_prefix 1>pregraph.log 2>pregraph.err#ORSOAPdenovo-63mer sparse_pregraph -s config_file -K 63 -z 5000000000 -R -o graph_prefix 1>pregraph.log 2>pregraph.err#step2:SOAPdenovo-63mer contig -g graph_prefix -R 1>contig.log 2>contig.err#step3:SOAPdenovo-63mer map -s config_file -g graph_prefix 1>map.log 2>map.err#step4:SOAPdenovo-63mer scaff -g graph_prefix -F 1>scaff.log 2>scaff.err

一步式：

SOAPdenovo-63mer all -s config_file -K 63 -R -o graph_prefix 1>ass.log 2>ass.err

輸出文件：運(yùn)行完會有不少的文件生成，其中后綴分別為contig和scafSeq即為對應(yīng)組裝結(jié)果，分別對應(yīng)contig和scaffold的結(jié)果。

4. SOAPdenovo實(shí)戰(zhàn)

下面我們找個NA12878樣本的測序數(shù)據(jù)，具體來實(shí)踐一下吧。

4.1 數(shù)據(jù)下載

下載測序數(shù)據(jù)：

wget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_1.fastq.gzwget http://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_2.fastq.gz

4.2 準(zhǔn)備配置文件

vi config_file， 填入以下內(nèi)容：

max_rd_len=100[LIB]avg_ins=200reverse_seq=0asm_flags=3rd_len_cutoff=100rank=1q1=NA12878.WGS-100K_1.fastq.gzq2=NA12878.WGS-100K_2.fastq.gz

4.3 運(yùn)行命令

此處我們將程序運(yùn)行的標(biāo)準(zhǔn)輸出和標(biāo)準(zhǔn)錯誤都分別重定向到對應(yīng)的log和err文件中了。一步式運(yùn)行：

SOAPdenovo-63mer all -s config_file -K 63 -R -o graph_prefix 1>ass.log 2>ass.err

四步單獨(dú)運(yùn)行：

SOAPdenovo-63mer pregraph -s config_file -K 63 -R -o graph_prefix 1>pregraph.log 2>pregraph.errSOAPdenovo-63mer contig -g graph_prefix -R 1>contig.log 2>contig.errSOAPdenovo2-r242/SOAPdenovo-63mer map -s config_file -g graph_prefix 1>map.log 2>map.errSOAPdenovo2-r242/SOAPdenovo-63mer scaff -g graph_prefix -F 1>scaff.log 2>scaff.err

4.4 輸出結(jié)果

此處我們的測試數(shù)據(jù)做了截取，因此可以非?？焖俚呐芡辏唧w的結(jié)果如下圖所示，可以看到生成了不少的中間結(jié)果文件，其中組裝出來的contig和scafford結(jié)果即圖上圈出來的兩個文件：

*.contig：contig序列文件，fasta格式；

*.scafSeq：scaffold序列文件，contig之間的gap用N填充 。

*.log和*.err是運(yùn)行的日志，里面包含很多的統(tǒng)計(jì)信息，如N50，N90，contig/Scaffold等信息。

*.scaf：包括scaffold中contig的詳細(xì)信息；在scaffold行中包括scaffold名字、contig長度和該scaffold長度。在contig行包括contig名字、contig在scaffold上的起始位置、正反鏈、長度和contig間的鏈接信息

*.links：contig間的pair-end連接信息

*.readOnContig：reads在contig上的位置。

4.5?從sixoclock下載soapdenovo2

此外，六點(diǎn)了官網(wǎng)基于CWL (common workflow language) 對SOAPdenovo2軟件進(jìn)行了封裝，通過我們開發(fā)的`sixbox` 軟件可以快速進(jìn)行軟件的運(yùn)行。對sixbox不了解可以通過六點(diǎn)了官網(wǎng)了解下。下面是具體的運(yùn)行步驟如下：

1)下載cwl 源碼

sixbox pull cad377c5-1a22-4a60-b761-d6e95e0d806b 或 在六點(diǎn)了官網(wǎng)上下載soapdenovo2.cwl

2)下載數(shù)據(jù)

wgethttp://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_1.fastq.gzwgethttp://www.sixoclock.net/resources/data/NGS/Homo_sapiens/WGS/NA12878.WGS-100K_2.fastq.gz

3)使用sixbox生成參數(shù)模板文件(YAML) , 并配置yaml文件

sixboxrun --make-template ./soapdenovo2.cwl > soapdenovo2.job.yamlvim soapdenovo2.job.yaml# 編輯參數(shù)配置文件，替換或設(shè)置參數(shù)以實(shí)現(xiàn)個性化分析

不熟悉的，可以直接粘貼下方示例內(nèi)容到soapdenovo2.job.yaml

q1:? # array oftype"File"(optional)-class: Filepath: NA12878.WGS-100K_1.fastq.gzq2:? # array oftype"File"-class: Filepath: NA12878.WGS-100K_2.fastq.gzo: graph_prefix? #defaultvalue oftype"string".K:63#defaultvalue oftype"int". (optional)R:true#defaultvalue oftype"boolean".rank:-1# array oftype"int"(optional)rd_len_cutoff:-100# array oftype"int"(optional)asm_flags:-3# array oftype"int"(optional)reverse_seq:-0# array oftype"int"(optional)avg_ins:-200# array oftype"int"(optional)max_rd_len:200#defaultvalue oftype"int". (optional)all:true#type"boolean"

4）使用sixbox運(yùn)行

sixboxrun ./soapdenovo2.cwl ./soapdenovo2.job.yaml#或者指定輸出目錄sixboxrun --outdir /home/path ./soapdenovo2.cwl ./soapdenovo2.job.yaml

運(yùn)行結(jié)束即可看到當(dāng)前目錄或者指定的輸出目錄輸出對應(yīng)的SOAPdenovo 組裝的結(jié)果文件。

至此，SOAPdenovo的實(shí)戰(zhàn)體驗(yàn)基本就結(jié)束了。

以上為我們給大家?guī)淼幕蚪Mdenovo的基本原理知識，以及在平臺上運(yùn)行經(jīng)典的SOAPdenovo的詳細(xì)操作過程。也歡迎大家去我們六點(diǎn)了官網(wǎng)看我們放上去的SOAPdenovo2的CWL流程工具。

如果對生物醫(yī)療健康大數(shù)據(jù)相關(guān)內(nèi)容感興趣也可以持續(xù)關(guān)注我們。想要探索更多的軟件流程或者知識文檔，可以到六點(diǎn)了官網(wǎng)查看。

References

[1]六點(diǎn)了官網(wǎng):http://www.sixoclock.net

[2]OLC方法 (Overlap-Layout-Consensus)和 DBG方法 (De-Bruijn Graph):https://zh.wikipedia.org/wiki/%E5%BA%8F%E5%88%97%E7%B5%84%E8%A3%9D

[3]歐拉圖（Euler Graph）問題:https://baike.baidu.com/item/歐拉圖/2587300

[4]SOAPdenovo:https://github.com/aquaskyline/SOAPdenovo2

[5]參考文獻(xiàn):http://www.genome.org/cgi/doi/10.1101/gr.097261.109

[6]github:https://github.com/aquaskyline/SOAPdenovo2

[7]二進(jìn)制:https://sourceforge.net/projects/soapdenovo2/files/SOAPdenovo2/