基于BERT和CNN的致病剪接突變預(yù)測(cè)方法

宋程程 趙依然 李曉艷 夏俊峰

剪接是真核生物一種重要的轉(zhuǎn)錄后調(diào)控過(guò)程,通過(guò)切除內(nèi)含子和連接外顯子可提高RNA(Ribo-nucleic Acid)和蛋白質(zhì)水平的多樣性和復(fù)雜性[1-3].多項(xiàng)研究表明,90%以上具有多個(gè)外顯子的基因會(huì)進(jìn)行可變剪接[4-5].但剪接過(guò)程可能受到基因組序列變異的影響或破壞,導(dǎo)致原本的剪接位點(diǎn)失去功能或通過(guò)破壞剪接體組裝以改變剪接過(guò)程[6-9].這種剪接異?？赡軙?huì)引起蛋白質(zhì)的合成異常,進(jìn)而導(dǎo)致疾病發(fā)生[10-11].

研究發(fā)現(xiàn)影響剪接的突變與多種疾病密切相關(guān)[12-15].然而,目前還沒(méi)有通用的規(guī)則解釋剪接突變的致病性.美國(guó)醫(yī)學(xué)遺傳學(xué)與基因組學(xué)會(huì)(Ameri-can College of Medical Genetics and Genomics,ACMG)在指南中指出,在功能喪失的已知機(jī)制的基因中,位于經(jīng)典剪接位點(diǎn)(外顯子和內(nèi)含子邊界的GT與AG二核苷酸位置)上的突變往往具有致病性[16].對(duì)于非經(jīng)典剪接突變,由于該區(qū)域無(wú)明顯的保守序列特征,且影響剪接的方式多種多樣[17-18],確定其致病性仍然具有挑戰(zhàn)性[19].因此,迫切需要開(kāi)發(fā)剪接突變致病性預(yù)測(cè)算法,尤其是針對(duì)非經(jīng)典剪接突變的算法,促進(jìn)在下一代疾病診斷測(cè)序中有效地對(duì)致病剪接突變進(jìn)行優(yōu)先排序.

現(xiàn)有的廣譜性致病突變預(yù)測(cè)方法可預(yù)測(cè)包含剪接突變?cè)趦?nèi)的多種突變的致病性,但由于該類(lèi)方法納入的剪接生物學(xué)特征有限,限制其對(duì)剪接突變致病性的預(yù)測(cè)[20-21].預(yù)測(cè)剪接分子效應(yīng)的工具也可對(duì)剪接突變致病性進(jìn)行預(yù)測(cè).例如:SpliceAI[22]是一個(gè)深度殘差神經(jīng)網(wǎng)絡(luò),可通過(guò)野生型和突變序列得分的差異預(yù)測(cè)剪接突變的致病性.類(lèi)似的預(yù)測(cè)剪接分子效應(yīng)的工具還有很多,如Maximum Entropy Mo-dels[23]、MMsplice(Modular Modeling of Splicing)[24]、CADD-Splice[25]、SPiP(Splicing Prediction Pipeline)[26]等.然而,預(yù)測(cè)剪接的分子效應(yīng)與預(yù)測(cè)剪接突變致病性是兩個(gè)截然不同的任務(wù).

目前,針對(duì)剪接突變致病性預(yù)測(cè),學(xué)者們已經(jīng)提出多種方法.S-CAP(Splicing Clinically Applicable Pathogenicity Prediction)[27]使用梯度增強(qiáng)樹(shù)分類(lèi)器構(gòu)建模型,定義與剪接位點(diǎn)相關(guān)的6個(gè)區(qū)域,包括外顯子及其邊界50 bp,用于預(yù)測(cè)影響RNA剪接的遺傳突變的致病性.SQUIRLS(Super Quick Information-Content Random-Forest Learning of Splice Variants)[19]為剪接供體和受體訓(xùn)練兩個(gè)隨機(jī)森林分類(lèi)器,并通過(guò)邏輯回歸組合它們的輸出,得出最終的剪接突變致病性預(yù)測(cè)分?jǐn)?shù).

上述預(yù)測(cè)致病剪接突變的方法大多基于傳統(tǒng)機(jī)器學(xué)習(xí)的方法,依賴(lài)人工提取的剪接生物學(xué)特征.這種依賴(lài)性可能會(huì)使某些剪接突變數(shù)據(jù)在預(yù)測(cè)過(guò)程中被忽略,導(dǎo)致方法預(yù)測(cè)的不夠全面.同時(shí),由于缺乏足夠的致病剪接突變正樣本數(shù)據(jù),現(xiàn)有深度學(xué)習(xí)工具中的傳統(tǒng)神經(jīng)網(wǎng)絡(luò)模型無(wú)法充分發(fā)揮優(yōu)勢(shì).此外,非經(jīng)典剪接突變往往難以識(shí)別,現(xiàn)有預(yù)測(cè)工具對(duì)非經(jīng)典剪接突變的預(yù)測(cè)精度仍有待提高.因此,亟需提出一種致病剪接突變預(yù)測(cè)方法,對(duì)剪接突變功能效應(yīng)進(jìn)行優(yōu)先排序,改善疾病的診斷和治療.

序列編碼內(nèi)部存在多義性和遙遠(yuǎn)的語(yǔ)義關(guān)系,這是自然語(yǔ)言的關(guān)鍵屬性.以往的研究表明,DNA(Deoxyribonucieic Acid)序列與人類(lèi)語(yǔ)言具有極大的相似性[28],也存在一詞多義的現(xiàn)象.例如:相同的順式調(diào)控元件(Cis-Regulatory Elements,CREs)在不同的上下游核苷酸序列中往往具有不同功能;間隔較遠(yuǎn)的多個(gè)CREs可能會(huì)協(xié)同,對(duì)啟動(dòng)子產(chǎn)生不同的功能作用.遺傳突變可能會(huì)通過(guò)影響CREs導(dǎo)致異常剪接[29-30].由此可見(jiàn),剪接突變的功能后果與其上下游序列密切相關(guān),尤其是那些發(fā)生在多義性的序列片段中的突變.

因此,本文提出基于BERT(Bidirectional Enco-der Representations from Transformer)[31]和CNN(Con-volutional Neural Networks)[32]的致病剪接突變預(yù)測(cè)方法(BERT and CNN-Based Deleterious Splicing Mu-tation Prediction Method,BCsplice).BCsplice中BERT模塊可全面提取序列的上下文信息,與提取局部特征的CNN結(jié)合后,可充分學(xué)習(xí)序列的語(yǔ)義信息,預(yù)測(cè)剪接突變致病性.非經(jīng)典剪接突變的影響往往更依賴(lài)序列上下文的深層語(yǔ)義信息,通過(guò)CNN將BERT的多級(jí)別語(yǔ)義信息進(jìn)行組合提取,可獲得豐富的信息表示,有助于識(shí)別非經(jīng)典剪接突變.對(duì)比實(shí)驗(yàn)表明BCsplice性能較優(yōu),尤其是在非經(jīng)典剪接區(qū)表現(xiàn)出一定性能優(yōu)勢(shì),有助于識(shí)別致病剪接突變和臨床遺傳診斷.

1 基于BERT和CNN的致病剪接

突變預(yù)測(cè)方法

基于Transformers 編碼器的語(yǔ)言模型BERT引入預(yù)訓(xùn)練和微調(diào)的范式,在許多自然語(yǔ)言處理任務(wù)中取得較優(yōu)性能.BERT首先從大量的無(wú)標(biāo)簽數(shù)據(jù)中發(fā)展通用的理解,再通過(guò)在特定任務(wù)上的微調(diào)取得較優(yōu)性能,即使使用少量的數(shù)據(jù)也能獲得出色的表現(xiàn).目前相關(guān)研究將BERT基于DNA序列進(jìn)行建模,稱(chēng)為DNABERT[33].因此,本文采用DNA-BERT的預(yù)訓(xùn)練模型作為BCsplice中提取序列上下文信息的一部分.此外,在使用BERT處理分類(lèi)任務(wù)時(shí),通常會(huì)使用BERT最后一層的信息.然而,Devlin等[31]對(duì)比BERT不同層的組合發(fā)現(xiàn),相比頂層的輸出,最后四個(gè)隱藏層的輸出組合可編碼更多的信息.因此,本文選擇提取BERT最后四個(gè)隱層的向量,并進(jìn)一步進(jìn)行組合處理.考慮到CNN在提取局部特征方面具有顯著優(yōu)勢(shì),以及經(jīng)典的TextCNN在文本分類(lèi)領(lǐng)域中的良好表現(xiàn),將BERT最后四個(gè)隱層的輸出與TextCNN結(jié)合,用于特征提取.

因此本文提出基于BERT和CNN的致病剪接突變預(yù)測(cè)方法(BCsplice),整體架構(gòu)如圖1所示.BCsplice首先處理數(shù)據(jù)集,再輸入序列.由于遺傳突變?cè)诓煌纳舷掠魏塑账嵝蛄兄锌赡軙?huì)產(chǎn)生不同的功能后果,突變位點(diǎn)的語(yǔ)義變化與上下文信息密切相關(guān).因此,首先提取突變位點(diǎn)的上下游基因組序列數(shù)據(jù),并將突變前后的基因組序列轉(zhuǎn)換為一組使用k-mer標(biāo)記表示的序列作為模型的輸入.序列中還包括一個(gè)[CLS]標(biāo)記(表示整個(gè)句子含義的標(biāo)簽)和兩個(gè)[SEP]標(biāo)記(表示句子分隔符,在突變前后的序列結(jié)尾分別添加).輸入序列首先經(jīng)過(guò)嵌入層,再送入BERT模塊.然后,將BERT最后四個(gè)隱層的[CLS]向量輸出進(jìn)行拼接,拼接后的向量通過(guò)TextCNN進(jìn)行卷積處理.最后,通過(guò)一個(gè)分類(lèi)層輸出剪接突變是否致病.

圖1 BCsplice框架圖

BCsplice主要包括4個(gè)模塊:嵌入模塊、BERT模塊、CNN模塊、分類(lèi)模塊.各模塊的細(xì)節(jié)描述如下.

1)嵌入模塊.輸入序列通過(guò)k-mer標(biāo)記表示后,每k個(gè)堿基相當(dāng)于一個(gè)詞(這里以k=6為例,即每6個(gè)堿基作為一個(gè)詞).嵌入模塊包含Token Embe-dding與Positional Embedding.Token Embedding將每個(gè)詞轉(zhuǎn)換為768維的向量,作為詞的向量表示.Positional Embedding對(duì)序列中每個(gè)詞的位置信息進(jìn)行編碼,使模型能學(xué)習(xí)輸入序列中每個(gè)詞的順序?qū)傩?最后將每個(gè)詞對(duì)應(yīng)位置的Token Embedding與Positional Embedding相加,得到嵌入向量矩陣.

2)BERT模塊.BERT的主要結(jié)構(gòu)是Transformer中的編碼器單元,這里共有12層編碼器.每個(gè)編碼器單元依次由多頭注意力機(jī)制、層歸一化、前饋神經(jīng)網(wǎng)絡(luò)和層歸一化疊加產(chǎn)生.

3)CNN模塊.使用3個(gè)相同大小的卷積核對(duì)BERT的輸出矩陣進(jìn)行卷積操作.卷積核的大小設(shè)置為2×768.卷積操作提取輸入矩陣中的局部特征.卷積后得到的特征向量通過(guò)ReLU激活函數(shù)進(jìn)行非線(xiàn)性映射,增強(qiáng)表示能力.然后,通過(guò)一個(gè)最大池化層對(duì)特征向量進(jìn)行降維,選擇每個(gè)卷積核輸出的最大值作為該卷積核的池化輸出.最后,將3個(gè)池化輸出進(jìn)行拼接,形成一個(gè)更豐富的特征表示.

4)分類(lèi)模塊.通過(guò)一個(gè)線(xiàn)性層預(yù)測(cè)輸入序列是否致病,得到最終的輸出結(jié)果.

2 實(shí)驗(yàn)及結(jié)果分析

2.1 數(shù)據(jù)集構(gòu)建

為了構(gòu)建一個(gè)高質(zhì)量的與剪接突變相關(guān)(Splicing-Associated Variant,SAV)的數(shù)據(jù)集,首先整合人類(lèi)基因突變數(shù)據(jù)庫(kù)(Human Gene Mutation Database,HGMD專(zhuān)業(yè)版2022.4)[34]上與剪接相關(guān)的突變數(shù)據(jù),共得到34 852條剪接突變數(shù)據(jù),其中包含突變類(lèi)型為splice、exonic-splice、canonical-splice的數(shù)據(jù).同時(shí)移除在Y染色體上的唯一一條剪接突變數(shù)據(jù),最終僅保留染色體1～22和X上的突變數(shù)據(jù).然后,根據(jù)人類(lèi)參考基因組版本為GRCh37和致病性標(biāo)簽為DM(Disease Causing Mutation),對(duì)剪接突變數(shù)據(jù)進(jìn)行篩選,最終獲得26 550條致病剪接突變數(shù)據(jù).另外,還從ClinVar(https://www.ncbi.nlm.nih.gov/clinvar,發(fā)布于2022年7月9日的VCF文件)[35]中獲取402 763條臨床意義為Benign的突變數(shù)據(jù).

由于本文的致病剪接突變數(shù)據(jù)僅包含單核苷酸突變,因此進(jìn)一步篩選突變類(lèi)型為單核苷酸的良性突變數(shù)據(jù)(GRCh37版本).同時(shí),為了保證負(fù)樣本數(shù)據(jù)的置信度較高,選擇ReviewStatus為criteria provi-ded,multiple submitters,no conflicts的突變數(shù)據(jù),并限定染色體為1～22和X.

經(jīng)過(guò)上述篩選,最終得到32 064條良性突變數(shù)據(jù).根據(jù)文獻(xiàn)[27],將剪接區(qū)域定義為外顯子及其周?chē)?0 bp,并篩選位于剪接區(qū)的突變數(shù)據(jù),從而獲得30 772條良性剪接突變數(shù)據(jù).

致病剪接突變與良性剪接突變數(shù)據(jù)的來(lái)源不同,導(dǎo)致存在少量數(shù)據(jù)標(biāo)簽矛盾.因此,移除正負(fù)樣本標(biāo)簽不一致的數(shù)據(jù),最終得到26 543條致病剪接突變數(shù)據(jù)和30 765條良性剪接突變數(shù)據(jù).另外,鑒于非經(jīng)典剪接突變的臨床意義較難解釋,根據(jù)剪接區(qū)域的不同(如表1所示),進(jìn)一步將剪接突變數(shù)據(jù)劃分到經(jīng)典剪接區(qū)(Canonical Splicing Region)與非經(jīng)典剪接區(qū)(Non-Canonical Splicing Region).

表1 數(shù)據(jù)集在不同剪接區(qū)的分布情況

為了減少來(lái)自同一基因的不同突變?cè)斐傻钠?采用Cytogenetic Band劃分訓(xùn)練集和測(cè)試集.首先,從總共721個(gè)染色體條帶中隨機(jī)選擇20%(144個(gè))的染色體條帶,將這些條帶中包含的突變作為測(cè)試集.然后,將剩下的80%(577個(gè))染色體條帶中包含的突變作為訓(xùn)練集.染色體條帶是指突變位點(diǎn)所處的染色體區(qū)段,每個(gè)突變都與特定的染色體區(qū)段關(guān)聯(lián).該方法的目的是減少位于同一基因中的不同突變同時(shí)被用于訓(xùn)練和測(cè)試的幾率,因?yàn)橄噜彽耐蛔兛赡芫哂邢嗨频奶卣?可能導(dǎo)致結(jié)果偏差[19].通過(guò)Cytogenetic Band,可確保同一染色體區(qū)段內(nèi)的突變要么用于訓(xùn)練,要么用于測(cè)試,而不會(huì)同時(shí)用于訓(xùn)練和測(cè)試.劃分后得到包含45 999條數(shù)據(jù)的訓(xùn)練集和包含11 309條數(shù)據(jù)的測(cè)試集,具體如表2所示.

表2 數(shù)據(jù)集劃分情況

2.2 參數(shù)設(shè)置及性能指標(biāo)

本文使用PyTorch框架構(gòu)建BCsplice.模型以DNABERT為基礎(chǔ)架構(gòu),采用其默認(rèn)參數(shù),詳細(xì)參數(shù)見(jiàn)表3.對(duì)于輸入的k-mer序列,文獻(xiàn)[35]研究表明,當(dāng)k=6時(shí),模型具有最優(yōu)性能,因此設(shè)置k=6,即每6個(gè)堿基作為一個(gè)詞.

表3 BCsplice參數(shù)信息

本文使用7個(gè)常用指標(biāo)評(píng)估BCsplice分別在全基因組區(qū)域和非經(jīng)典剪接區(qū)的性能表現(xiàn).具體指標(biāo)如下:Accuracy(ACC)、ROC(Receiver Operating Cha-racteristic Curve)曲線(xiàn)下面積(Area under the ROC Curve,AUC),精確率-召回率曲線(xiàn)下面積(Area under the Precision-Recall Curve,AUPR)、F1 score(F1)、馬修斯相關(guān)系數(shù)(Matthews Correlation Coeffi-cient,MCC)、精確率(Precision)、召回率(Recall).這些評(píng)價(jià)指標(biāo)是根據(jù)混淆矩陣計(jì)算的,參數(shù)包括真正例(True Positive,TP)、真反例(True Negative,TN)、假正例(False Positive,FP)和假反例(False Nega-tive,FN),具體計(jì)算方式定義如下:

一般來(lái)說(shuō),各指標(biāo)值越高,模型性能越優(yōu).在進(jìn)行性能對(duì)比時(shí),主要考慮AUC與AUPR,因?yàn)樵谶@7個(gè)評(píng)價(jià)指標(biāo)中,AUC與AUPR是最重要的2個(gè)綜合評(píng)價(jià)指標(biāo).

2.3 不同剪接區(qū)的突變數(shù)據(jù)分布

經(jīng)典剪接區(qū)的突變位于外顯子與內(nèi)含子邊界±2 bp位置,具有明顯的基序特征和位置特征,通常容易解釋.非經(jīng)典剪接區(qū)的突變相對(duì)更隱蔽,致病性較難確定.因此,本文將突變分為2個(gè)區(qū)域:經(jīng)典剪接區(qū)和非經(jīng)典剪接區(qū).外顯子和內(nèi)含子邊界的GT與AG二核苷酸位置的突變是發(fā)生在經(jīng)典剪接區(qū)的突變,這個(gè)位置的突變絕大多數(shù)都是致病的[16].在本文設(shè)計(jì)的數(shù)據(jù)集上,經(jīng)典剪接區(qū)的致病突變有21 965條,而良性突變只有25條.由此可見(jiàn),經(jīng)典剪接區(qū)上的良性突變極少,僅占0.08%,良性突變基本上都位于非經(jīng)典剪接區(qū).此外,根據(jù)本文的分類(lèi)結(jié)果,經(jīng)典剪接區(qū)的正負(fù)樣本非常不平衡.盡管非經(jīng)典剪接區(qū)的正負(fù)樣本也存在不平衡,但相比之下,在非經(jīng)典剪接區(qū)發(fā)生的突變并不像經(jīng)典剪接區(qū)那樣直觀,其致病性需要進(jìn)一步研究.基于不同區(qū)域的致病性情況和突變數(shù)據(jù)的分布情況,在接下來(lái)的性能對(duì)比中主要關(guān)注非經(jīng)典剪接區(qū)和全基因組區(qū)域(包括經(jīng)典剪接區(qū)和非經(jīng)典剪接區(qū)).測(cè)試集在非經(jīng)典剪接區(qū)和全基因組區(qū)域上的具體分布情況如表4所示.

表4 測(cè)試集在不同剪接區(qū)的分布情況

2.4 消融實(shí)驗(yàn)

為了研究BCsplice中BERT的隱層向量經(jīng)過(guò)TextCNN處理后對(duì)模型性能提升的重要性,定義BERT-TextCNN和BERT,進(jìn)行針對(duì)TextCNN模塊的消融實(shí)驗(yàn).由于經(jīng)典剪接突變通常會(huì)對(duì)剪接產(chǎn)生有害影響,故預(yù)測(cè)非經(jīng)典剪接突變的致病性更具挑戰(zhàn)性.因此,除了在全基因組區(qū)域進(jìn)行分析以外,還對(duì)非經(jīng)典剪接區(qū)進(jìn)行單獨(dú)分析.在列出的評(píng)價(jià)指標(biāo)中,AUC與AUPR是最重要的兩個(gè)綜合評(píng)價(jià)指標(biāo),因此主要根據(jù)這兩個(gè)指標(biāo)進(jìn)行性能對(duì)比.

在全基因組區(qū)域,首先對(duì)模型進(jìn)行五折交叉驗(yàn)證和測(cè)試數(shù)據(jù)集上的性能評(píng)估,具體如表5所示.根據(jù)表中結(jié)果可見(jiàn),BERT-TextCNN在各個(gè)評(píng)價(jià)指標(biāo)上都有性能提升.而在非經(jīng)典剪接區(qū)域,BERT-TextCNN和BERT的指標(biāo)值對(duì)比如表6所示.由表可見(jiàn),BERT-TextCNN性能提升更顯著,BERT-Text-CNN的AUC值為0.947,AUPR值為0.828,相比BERT分別提高5.1%和15.2%.

表5 各模型在五折交叉驗(yàn)證和測(cè)試集上的性能

表6 各模型在非經(jīng)典剪接區(qū)測(cè)試集上的性能

每個(gè)編碼層的句子向量表示不同級(jí)別的語(yǔ)義信息.淺層的編碼層能更好地捕捉低級(jí)別的語(yǔ)義信息,而深層的編碼層能更好地捕捉高級(jí)別的語(yǔ)義信息.通過(guò)TextCNN將多個(gè)編碼層的語(yǔ)義信息進(jìn)行組合提取,既能獲取原始的全局語(yǔ)義信息,又能獲得序列信息中關(guān)鍵的特征.這可能是模型性能提升的重要原因之一.此外,非經(jīng)典剪接區(qū)突變的影響更隱晦,對(duì)于序列上下文的深層語(yǔ)義信息更依賴(lài).

因此,提取多級(jí)別的隱含信息對(duì)于預(yù)測(cè)非經(jīng)典剪接突變的致病性具有重要意義.

2.5 對(duì)比實(shí)驗(yàn)

為了進(jìn)一步驗(yàn)證BCsplice的優(yōu)越性,選擇如下剪接突變預(yù)測(cè)方法進(jìn)行對(duì)比:SQUIRLS[19]、Spli-ceAI[22]、CADD-Splice[25]、SPiP[26]、S-CAP[27]、MLC-splice(Machine Learning-Based Classification of Splice Sites Variants)[36].

SpliceAI是一種可從pre-mRNA轉(zhuǎn)錄本序列中預(yù)測(cè)剪接連接(Splice Junctions)的方法,通過(guò)對(duì)比突變前后序列的分?jǐn)?shù)差異以預(yù)測(cè)突變的致病性.CADD-splice整合其它工具的剪接分?jǐn)?shù),并優(yōu)化對(duì)剪接突變效應(yīng)的預(yù)測(cè).MLCsplice是一個(gè)集成的機(jī)器學(xué)習(xí)方法,用于預(yù)測(cè)突變對(duì)剪接的影響.SPiP也是一種機(jī)器學(xué)習(xí)工具,用于檢測(cè)外顯子和內(nèi)含子突變對(duì)mRNA剪接的影響.

為了對(duì)比方法的性能,使用11 309條測(cè)試數(shù)據(jù),并在全基因組區(qū)域和非經(jīng)典剪接區(qū)域進(jìn)行實(shí)驗(yàn),結(jié)果如表7和表8所示.

表7 各方法在全基因組區(qū)上的性能對(duì)比

表8 各方法在非經(jīng)典剪接區(qū)上的性能對(duì)比

在本文的測(cè)試數(shù)據(jù)中,并不是每個(gè)工具都能對(duì)所有突變給出預(yù)測(cè)結(jié)果.因此將對(duì)比預(yù)測(cè)工具分為兩類(lèi):1)不能對(duì)所有突變給出預(yù)測(cè)結(jié)果的工具,即存在缺失突變(CADD-splice、SpliceAI、S-CAP、MLC-splice);2)可對(duì)所有測(cè)試突變給出預(yù)測(cè)結(jié)果,即沒(méi)有缺失突變(SPiP、SQUIRLS).在表中使用“類(lèi)別”進(jìn)行標(biāo)注以區(qū)分這兩類(lèi)工具.需要特別注意的是,BCsplice也屬于沒(méi)有缺失突變的這一類(lèi)工具.因此BCsplice與SPiP、SQUIRLS可在同一水平上進(jìn)行對(duì)比.

在全基因區(qū)域,使用由5 319個(gè)致病剪接突變和5 990個(gè)良性突變組成的測(cè)試集.根據(jù)表7的評(píng)估結(jié)果顯示,BCsplice在AUC和AUPR指標(biāo)上取得最高值(AUC=0.987,AUPR=0.988).相比存在缺失突變的工具,BC-splice表現(xiàn)出明顯的優(yōu)勢(shì).BCsplice能對(duì)所有突變數(shù)據(jù)進(jìn)行預(yù)測(cè),而其余4種方法都存在一定的缺失突變.具體缺失突變的數(shù)量見(jiàn)表9.MLCsplice缺失的突變數(shù)量最多,達(dá)到10 092條(10092÷11309≈89%).其次是S-CAP,缺失5 076條突變(5076÷11309≈45%).相比沒(méi)有缺失突變的工具,BCsplice在AUC和AUPR值上的提升約為1%～2%(見(jiàn)圖2),也顯示出一定的優(yōu)勢(shì).

表9 各方法在全基因組區(qū)上的混淆矩陣值

(a)AUC

綜合來(lái)看,相比對(duì)比方法,BCsplice性能更優(yōu)或相當(dāng).這是因?yàn)閷?duì)于使用者來(lái)說(shuō),一個(gè)理想的方法不僅需要高精確度,還需要能分析更廣泛范圍的突變.BCsplice可對(duì)全基因組上任何區(qū)域的突變進(jìn)行預(yù)測(cè),而當(dāng)前的預(yù)測(cè)工具大多存在剪接突變預(yù)測(cè)分?jǐn)?shù)缺失的情況.那些能對(duì)所有剪接突變進(jìn)行預(yù)測(cè)的工具,在性能上并未高于BCsplice,特別是在非經(jīng)典剪接區(qū)域中.

另外需要說(shuō)明的是,由于獲取MLCsplice和S-CAP的預(yù)測(cè)分?jǐn)?shù)需要先計(jì)算其它多種工具的特征分?jǐn)?shù),過(guò)程較耗時(shí),因此本文直接使用預(yù)先計(jì)算好的MLCsplice和S-CAP的人類(lèi)基因組突變分?jǐn)?shù),這可能是兩個(gè)預(yù)測(cè)工具缺失值較多的原因.對(duì)于CADD-Splice,本文使用線(xiàn)上預(yù)測(cè)功能,下載包含PHRED標(biāo)注突變致病性評(píng)分的tsv文件.為了評(píng)估SpliceAI,根據(jù)其文獻(xiàn)中提供的網(wǎng)站資源,在服務(wù)器上下載SpliceAI工具包,SpliceAI為每個(gè)突變提供4個(gè)Δ分?jǐn)?shù),這4個(gè)分?jǐn)?shù)中的最大值代表該突變的ΔScore.然而,對(duì)于CADD-Splice和SpliceAI,目前無(wú)法明確它們存在缺失值的具體原因.可能是由于它們依賴(lài)于某些剪接特征或受到基因組區(qū)域的限制,導(dǎo)致無(wú)法預(yù)測(cè)所有的剪接突變.對(duì)于所有對(duì)比方法,在本文的測(cè)試集上,使用find_optimal_cutoff函數(shù)生成具有高AUC的最佳閾值,確保性能對(duì)比的準(zhǔn)確性與一致性.

非經(jīng)典剪接突變的臨床意義較難解釋,發(fā)生在非經(jīng)典剪接區(qū)的突變也難以識(shí)別.大量的預(yù)測(cè)方法專(zhuān)注于經(jīng)典剪接區(qū),對(duì)非經(jīng)典剪接區(qū)突變的預(yù)測(cè)能力較弱.因此本文在非經(jīng)典剪接區(qū)上也進(jìn)行性能對(duì)比,具體如表10所示.在本文的測(cè)試集上,非經(jīng)典剪接區(qū)的致病剪接突變有916條,良性突變有5 985條.

表10 各方法在非經(jīng)典剪接區(qū)上的混淆矩陣值

由表10可見(jiàn),相比預(yù)測(cè)區(qū)域在全基因組區(qū)域,所有方法在非經(jīng)典剪接區(qū)上都顯示出一定程度的性能下降.其中CADD-Splice、SPiP與SQUIRLS性能下降最明顯,這可能與非經(jīng)典剪接區(qū)的樣本數(shù)據(jù)不平衡有關(guān).在非經(jīng)典剪接區(qū)域,本文依舊將對(duì)比方法分為兩類(lèi).相比存在缺失突變的方法,MLCsplice在AUC和AUPR指標(biāo)上具有最高值.然而,在本文的測(cè)試集上,MLCsplice存在5 994個(gè)突變的缺失情況(占總數(shù)的87%),具體缺失情況見(jiàn)表10.相比沒(méi)有缺失突變的方法,BCsplice具有顯著優(yōu)勢(shì)(AUC=0.947,AUPR=0.828).相比SPiP(AUC=0.855,AUPR=0.735)與SQUIRLS(AUC=0.864,AUPR=0.705),BCsplice的綜合性能提升8%～12%,具體如圖3所示.

因此,BCsplice在不同情況下都表現(xiàn)出獨(dú)特的優(yōu)勢(shì).非經(jīng)典剪接突變通常更依賴(lài)序列上下文的深層語(yǔ)義信息.BCsplice通過(guò)結(jié)合CNN和BERT以捕捉多級(jí)別的語(yǔ)義信息,這可能是其性能較優(yōu)的原因之一.

3 結(jié) 束 語(yǔ)

本文提出基于BERT和CNN的致病剪接突變預(yù)測(cè)方法(BCsplice),用于預(yù)測(cè)剪接突變的致病性.BCsplice利用基因組序列數(shù)據(jù)全面捕捉序列上下文信息,鑒別有害的剪接突變.該方法顯示自身的優(yōu)越性與實(shí)用性,尤其是對(duì)于非經(jīng)典剪接突變的預(yù)測(cè),可支持改善和擴(kuò)展對(duì)剪接突變的診斷能力,并對(duì)突變產(chǎn)生的剪接影響進(jìn)行評(píng)估,為識(shí)別剪接突變做出一定的貢獻(xiàn).此外,BCsplice是一個(gè)深度學(xué)習(xí)算法,通過(guò)自動(dòng)學(xué)習(xí)序列信息以提取特征,但是這些特征并不是容易理解的有明確生物學(xué)含義的特征.今后的研究方向之一是進(jìn)一步探索方法的可解釋性,并對(duì)突變導(dǎo)致的剪接缺陷進(jìn)行解釋說(shuō)明,這有助于理解剪接突變?nèi)绾螌?dǎo)致人類(lèi)疾病,并更深入了解剪接這一精密過(guò)程的調(diào)節(jié)機(jī)制.