群體基因組結(jié)構(gòu)變異檢測工作流

曹舒淇，劉詩琦，姜 濤

(哈爾濱工業(yè)大學(xué) 計算學(xué)部，哈爾濱150001)

基因組結(jié)構(gòu)變異(Structural Variation，SV)是基因組上大尺度的核苷酸序列重排性變化，它包括長度大于50 bp的插入(INS)、缺失(DEL)、倒位(INV)、重復(fù)(DUP)、易位(BND)[1]。相關(guān)研究表明，平均每個人類個體上存在大約兩萬個結(jié)構(gòu)變異[2]，結(jié)構(gòu)變異盡管相較于單核苷酸變異(SNV)、短插入缺失變異(INDEL)數(shù)量較少，但因其變異長度較大，因此對基因組上核苷酸序列的影響是最廣泛的[3]。結(jié)構(gòu)變異會改變基因序列信息，進而影響轉(zhuǎn)錄過程，改變蛋白質(zhì)空間結(jié)構(gòu)，從而引發(fā)性狀與表型的改變[4]。此外，結(jié)構(gòu)變異對基因表達調(diào)控[5]、種群多樣性[6]等方面有著重要影響，同時與以自閉癥[7]、阿爾茲海默癥[8]等為代表的許多疾病的引發(fā)有密切的關(guān)系。

結(jié)構(gòu)變異會對人類遺傳、進化產(chǎn)生影響，形成個體之間的差異，影響種群的發(fā)展與演進。對于同一個群體，相當(dāng)數(shù)量的結(jié)構(gòu)變異對于群體中大部分個體是共享的，這些共享的結(jié)構(gòu)變異可以有效對群體的特征與結(jié)構(gòu)進行刻畫[9]。此外，在群體中仍存在個體特有的結(jié)構(gòu)變異，這些個體特有的結(jié)構(gòu)變異反映了個體獨有的特性，通過對特有結(jié)構(gòu)變異以及個體表型的分析，能夠發(fā)掘結(jié)構(gòu)變異與表型、疾病之間的重要關(guān)系[10]。

隨著國際千人基因組計劃的實施與推動[11-13]，各國也紛紛啟動了本國的大規(guī)模人群基因組計劃[14-17]，希望通過分析和構(gòu)建本國、本民族的基因組變異圖譜，更加深入地解讀本國人群在遺傳、進化上的機理，為接下來開展的疾病診治、精準(zhǔn)健康發(fā)展提供支撐。結(jié)構(gòu)變異作為對基因序列影響最為廣泛的基因組變異類型，如何高效、精準(zhǔn)的檢測群體結(jié)構(gòu)變異已成為當(dāng)前群體基因組研究中的核心。因此我們基于多層過濾的質(zhì)量控制，多種算法的聯(lián)合檢測、多維度變異融合和校對，開發(fā)了一個高性能的群體結(jié)構(gòu)變異檢測工作流，實現(xiàn)了群體基因組結(jié)構(gòu)變異的全面、精準(zhǔn)檢測。該工作流總體分為四個環(huán)節(jié)：基因組測序片段比對，單樣本基因組結(jié)構(gòu)變異檢測，單樣本基因組結(jié)構(gòu)變異融合以及群體基因組結(jié)構(gòu)變異檢測(見圖1)。

圖1 群體基因組結(jié)構(gòu)變異檢測工作流Fig.1 Workflow of structural variation detection from population genomes

1 基因組測序片段比對

高通量基因組測序片段比對是基因組數(shù)據(jù)分析的首要環(huán)節(jié)，測序片段的比對的精度將對變異檢測、基因組拼接等下游分析產(chǎn)生重要的影響。因此，對基因組片段測序數(shù)據(jù)、片段比對數(shù)據(jù)等有效的質(zhì)量控制，是保障以測序片段比對為基礎(chǔ)的基因組數(shù)據(jù)分析的關(guān)鍵。為此，本研究設(shè)計了多重質(zhì)量控制與過濾的基因組測序片段比對工作流，該工作流主要包含以下步驟(見圖2)。

圖2 基因組測序片段比對流程Fig.2 Workflow of sequencing read alignment

(1)使用測序片段質(zhì)量評價算法FastQC(https://github.com/s-andrews/FastQC)(v0.11.9，默認(rèn)參數(shù))對各樣本基因組測序數(shù)據(jù)進行質(zhì)量控制，通過對測序片段中GC含量、重復(fù)性、堿基質(zhì)量、片段長度分布等指標(biāo)進行統(tǒng)計和閾值判定，若任意滿足：測序片段GC含量與理論分布偏差30%；測序片段重復(fù)度超過理論重復(fù)總量的50%；測序片段任意位置的堿基質(zhì)量下四分位數(shù)低于5或中位數(shù)低于20；任意測序片段長度不足或長于150 bp，則將其認(rèn)定為低質(zhì)量測序樣本數(shù)據(jù)，并進行過濾處理。

(2)使用高通量測序片段比對算法BWA (https://github.com/lh3/bwa)(v0.7.17，默認(rèn)參數(shù))，完成各測序樣本向參考基因組序列的比對。使用比對格式轉(zhuǎn)換算法Sambamba[18](v0.8.0，默認(rèn)參數(shù))對各樣本比對結(jié)果進行格式轉(zhuǎn)換和排序；使用測序重復(fù)片段標(biāo)記算法Samblaster[19](v0.1.2，默認(rèn)參數(shù))對各樣本轉(zhuǎn)換后的比對文件進行重復(fù)標(biāo)記；使用GATK (https://github.com/broadinstitute/gatk)(v4.2.0.0，默認(rèn)參數(shù))對測序片段比對中堿基質(zhì)量校正形成最終的片段比對數(shù)據(jù)。

(3)使用測序片段比對質(zhì)量評價算法Qualimap[20](v2.2.1，默認(rèn)參數(shù))對各樣本測序片段比對結(jié)果進行質(zhì)量控制，通過對片段比對中測序覆蓋度(不低于30×測序深度)、片段重復(fù)性(不高于5%)、片段比對率(不低于95%)等指標(biāo)進行統(tǒng)計和閾值判定，進一步過濾低質(zhì)量測序樣本數(shù)據(jù)。

(4)使用DNA污染估計算法Verifybamid[21](v2.0.1，默認(rèn)參數(shù))計算各樣本中DNA污染程度，過濾高污染率(高于3%)測序樣本數(shù)據(jù)，形成最終用于下游分析的群體樣本集合。

多重質(zhì)量控制與過濾的基因組測序片段比對工作流在完成各樣本測序片段比對任務(wù)的同時，將有效監(jiān)控測序數(shù)據(jù)質(zhì)量、比對數(shù)據(jù)質(zhì)量、樣本污染情況等多重指標(biāo)，為高質(zhì)量基因組結(jié)構(gòu)變異檢測奠定基礎(chǔ)。

2 單樣本基因組結(jié)構(gòu)變異檢測

受限于高通量測序數(shù)據(jù)讀長與系統(tǒng)性測序誤差的限制，采用單一結(jié)構(gòu)變異檢測工具識別各樣本基因組中的結(jié)構(gòu)變異往往存在敏感性與準(zhǔn)確性較低的問題，這將制約結(jié)構(gòu)變異的檢測能力和向下游研究轉(zhuǎn)化的水平。針對這一問題，本研究采用三款當(dāng)前性能最好的個體基因組結(jié)構(gòu)變異檢測算法，全面挖掘各樣本基因組結(jié)構(gòu)變異，主要步驟如下(見圖3)。

圖3 單樣本基因組結(jié)構(gòu)變異檢測流程Fig.3 Workflow of structural variation detection from individual sample

(1)使用快速檢測基因組結(jié)構(gòu)變異檢測算法Manta[22](v1.6.0，默認(rèn)參數(shù))識別各基因組中DEL變異、INS變異、INV變異、BND變異、DUP變異。

(2)使用簡化集成基因組結(jié)構(gòu)變異檢測與基因分型算法Smoove(https://github.com/brentp/smoove)(v0.2.7，默認(rèn)參數(shù))識別各基因組中DEL變異、INS變異、INV變異、BND變異、DUP變異。

(3)使用拷貝數(shù)變異檢測算法CNVNator[23](v0.4.1，默認(rèn)參數(shù))識別各基因組中的拷貝數(shù)變異(CNV)，并計算各基因組區(qū)域上測序覆蓋度信息。

分別使用Manta、Smoove、CNVNator三種結(jié)構(gòu)變異檢測算法，將有效挖掘每個樣本基因組中多種類型結(jié)構(gòu)變異，為融合形成群體基因組結(jié)構(gòu)變異提供支撐。

3 單樣本基因組結(jié)構(gòu)變異融合

群體基因組結(jié)構(gòu)變異主要由各單樣本基因組結(jié)構(gòu)變異融合產(chǎn)生。如何對來自不同樣本、不同檢測算法形成的結(jié)構(gòu)變異進行融合，是當(dāng)前產(chǎn)生高精度群體基因組結(jié)構(gòu)變異的核心。為此，本研究分別對群體樣本中由同種檢測算法與不同檢測算法預(yù)測的結(jié)構(gòu)變異分層次整合，從而產(chǎn)生最終群體結(jié)構(gòu)變異候選位點，主要步驟如下(見圖4)。

圖4 單樣本基因組結(jié)構(gòu)變異融合流程Fig.4 Workflow of integration of individual structural variation

(1)分別對由Manta、Smoove檢測算法產(chǎn)生的結(jié)構(gòu)變異按照基因組坐標(biāo)進行排序，完成在相同檢測算法上不同樣本結(jié)構(gòu)變異的融合。若相鄰兩個結(jié)構(gòu)變異存在交疊，則將兩個變異合并為一個變異，直至所有變異均不存在交疊性。對于合并后的結(jié)構(gòu)變異，分別記錄來源樣本標(biāo)號，同時以累加方式累積各來源樣本在此變異上的變異質(zhì)量數(shù)、測序片段支持度等信息。

(2)再次對由Manta、Smoove檢測算法產(chǎn)生的基因組結(jié)構(gòu)變異融合結(jié)果進行排序，完成對不同檢測算法產(chǎn)生的融合結(jié)構(gòu)變異數(shù)據(jù)的二次融合。

通過使用不同工具對各樣本基因組中的結(jié)構(gòu)變異雙重融合，在充分保留各樣本基因組中潛在的結(jié)構(gòu)變異的同時，在融合過程中記錄每個結(jié)構(gòu)變異的樣本支持情況、變異質(zhì)量情況等信息，為過濾低質(zhì)量群體結(jié)構(gòu)變異提供了保障。

4 群體基因組結(jié)構(gòu)變異檢測

在完成單樣本基因組結(jié)構(gòu)變異融合后，進行質(zhì)量控制和多重變異屬性過濾，是完成高質(zhì)量群體基因組結(jié)構(gòu)變異檢測，繪制高精度人群基因組結(jié)構(gòu)變異圖譜的核心。本研究實現(xiàn)這一目標(biāo)主要采用如下3個步驟(見圖5)。

圖5 群體基因組結(jié)構(gòu)變異檢測流程Fig.5 Workflow of structural variation detection from population genomes

(1)將融合后的基因組結(jié)構(gòu)變異按照變異類型分別拆分為DEL、INS、INV、DUP、BND五種類型。使用基因組結(jié)構(gòu)變異斷點基因型計算算法SVTyper[24](v0.7.1，默認(rèn)參數(shù))，分別從單樣本層面對以上融合后的五種類型結(jié)構(gòu)變異重新計算基因型。使用CNVNator計算的單樣本層面的測序覆蓋度信息對重新校準(zhǔn)基因型信息的各結(jié)構(gòu)變異進行注釋。

(2)對重新校正基因型和測序覆蓋度信息的各類型結(jié)構(gòu)變異合并，計算各變異在人群中的變異頻率。將合并后群體結(jié)構(gòu)變異檢測結(jié)果轉(zhuǎn)換為bedpe文件格式并排序，對存在變異區(qū)域交疊的結(jié)構(gòu)變異進行聚類，保留聚類中具有最大變異頻率的結(jié)構(gòu)變異，將聚類中其余結(jié)構(gòu)變異修剪刪除。

(3)將經(jīng)過修剪的結(jié)構(gòu)變異集合重新轉(zhuǎn)換為vcf格式，依據(jù)測序覆蓋度信息和基因型一致性信息對各結(jié)構(gòu)變異重新校對變異類型，新增移動元件變異(MEI)類型，形成最終的群體結(jié)構(gòu)變異檢測結(jié)果。

經(jīng)過對不同類型結(jié)構(gòu)變異基因型的重新校正、過濾和變異類型校對，有效消減檢測形成的假陽性結(jié)構(gòu)變異預(yù)測結(jié)果，最大限度反映群體中真實的結(jié)構(gòu)變異位點、類型和變異頻率，為最終繪制群體基因組變異圖譜提供了堅實的保障。

5 結(jié)果與分析

為了驗證群體基因組結(jié)構(gòu)變異檢測工作流的真實效果，本研究構(gòu)建了由267個樣本組成的人群，使用Illumina高通量測序平臺對該人群樣本進行了30×高深度全基因組測序，并使用本研究提出的群體基因組結(jié)構(gòu)變異檢測工作流對此267個樣本進行群體結(jié)構(gòu)變異檢測(見表1)，合計檢測出了96 202個結(jié)構(gòu)變異，其中包括：11 697個DEL變異、18 385個INS變異、3 563個DUP變異、1 278個INV變異、2 007個MEI變異、59 272個BND變異。

表1 267樣本人群中結(jié)構(gòu)變異檢測結(jié)果統(tǒng)計Table 1 Results of structural variation detection of 267 samples

在該267個樣本構(gòu)成的人群中(見圖6)，常見變異(AF≥0.05)占總體檢出變異的41%(39 086/96 202)，低頻變異(0.05>AF≥0.01)占總體檢出變異的18%(17 554/96 202)，罕見變異(0.01>AF)占總體檢出變異的41%(39 562/96 202)。值得關(guān)注的是，在DEL、DUP、INS、INV四種類型結(jié)構(gòu)變異中，罕見變異的占比基本是均超過總體檢驗出變異的50%，相比之下，在MEI、BND兩種類型結(jié)構(gòu)變異中，檢測出的常見變異數(shù)量是總體可檢測變異數(shù)量的主要占比。這些結(jié)果與過去開展的基因組計劃發(fā)現(xiàn)的結(jié)果相一致[25-26]，說明本研究建立的群體基因組結(jié)構(gòu)變異檢測工作流具有良好的檢測能力。

圖6 不同變異頻率中結(jié)構(gòu)變異分布統(tǒng)計Fig.6 The distribution of structural variation among various allele frequencies

就每個樣本可檢測的結(jié)構(gòu)變異而言，平均每個樣本可以檢測出18 388個結(jié)構(gòu)變異，其中包含1 634個DEL變異、657個DUP變異、1 216個MEI變異、13 155個BND變異、1 521個INS變異、206個INV變異(見圖7)。受限于高通量測序技術(shù)中讀長的限制，基因組重復(fù)片段區(qū)域中的結(jié)構(gòu)變異難以檢測和精確分型，其中僅可檢測到斷點連接關(guān)系的結(jié)構(gòu)變異均歸結(jié)為BND變異，因此導(dǎo)致了每個樣本中包含了相當(dāng)數(shù)量的BND變異。然而，僅獲取變異斷點連接關(guān)系，無法解析結(jié)構(gòu)變異精準(zhǔn)結(jié)構(gòu)(如：是否為平衡變異，DNA變化方向等)將嚴(yán)重影響B(tài)ND變異的可信度和準(zhǔn)確性。經(jīng)過對BND變異按照置信度進行過濾(見圖7)，總計移除55 492個BND變異，僅保留3 780個高置信度BND變異(移除率93.62%)。平均每個樣本移除11 703個BND變異，僅保留1 451個高置信度BND變異(移除率88.97%)。

圖7 各樣本不同類型結(jié)構(gòu)變異檢測數(shù)量分布和統(tǒng)計Fig.7 The quantitative distribution of various structural variation types among samples

此外，本研究還對群體基因組結(jié)構(gòu)變異檢測工作流中變異融合與群體變異檢測兩個關(guān)鍵環(huán)節(jié)的計算開銷和內(nèi)存使用進行了統(tǒng)計(見表2)。該工作流完成群體基因組結(jié)構(gòu)變異檢測和融合兼容串行分析和并行分析兩種方式，其中串行計算方式需要約173.4 h，最大內(nèi)存開銷30 GB，而采用并行計算方式僅需不足3 h，并維持最大30 GB的內(nèi)存開銷。這一結(jié)果表明，對于大規(guī)模人群基因組結(jié)構(gòu)變異檢測分析，在保持有限內(nèi)存消耗的前提下，采用并行方式運行該工作流將顯著提升計算速度，為高效、快速的群體基因組結(jié)構(gòu)變異檢測提供了保證。

表2 267樣本人群中結(jié)構(gòu)變異檢測運行時間及內(nèi)存統(tǒng)計Table 2 Time and memory cost of structural variation detection of 267 samples

6 結(jié) 論

1)本研究構(gòu)建了一套高效、精準(zhǔn)的群體基因組結(jié)構(gòu)變異檢測工作流，該工作流通過多層過濾的質(zhì)量控制，為高質(zhì)量群體基因組結(jié)構(gòu)變異檢測提供支撐。

2)該工作流通過使用多種高性能結(jié)構(gòu)變異檢測算法，提高了結(jié)構(gòu)變異檢測的準(zhǔn)確性與敏感性，并通過雙重融合實現(xiàn)了群體結(jié)構(gòu)變異候選位點的精準(zhǔn)定位。

3)該工作流通過多維度重新校正結(jié)構(gòu)變異候選位點的基因型與變異類別，進一步保障群體結(jié)構(gòu)變異圖譜的高質(zhì)量構(gòu)建。

4)利用該工作流對由267個樣本組成的人群進行基因組結(jié)構(gòu)變異檢測，結(jié)果表明該工作流具有良好、快速、高效的檢測能力。通過并行分析策略在控制內(nèi)存消耗的基礎(chǔ)上，提高了工作流的計算速度，為大規(guī)模群體基因組研究提供了可能。