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

曹舒淇，劉詩琦，姜 濤

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

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

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

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

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

1 基因組測(cè)序片段比對(duì)

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

圖2 基因組測(cè)序片段比對(duì)流程Fig.2 Workflow of sequencing read alignment

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)果與分析

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

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

在該267個(gè)樣本構(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)變異中，罕見變異的占比基本是均超過總體檢驗(yàn)出變異的50%，相比之下，在MEI、BND兩種類型結(jié)構(gòu)變異中，檢測(cè)出的常見變異數(shù)量是總體可檢測(cè)變異數(shù)量的主要占比。這些結(jié)果與過去開展的基因組計(jì)劃發(fā)現(xiàn)的結(jié)果相一致[25-26]，說明本研究建立的群體基因組結(jié)構(gòu)變異檢測(cè)工作流具有良好的檢測(cè)能力。

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

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

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

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

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

6 結(jié) 論

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

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

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

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