單腫瘤組織微衛(wèi)星不穩(wěn)定探測方法①

趙 丹,尚秋明,韓鑫胤,李瑞琳,何小雨,祝海棟,牛北方

1(中國科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心,北京 100190)

2(中國科學(xué)院大學(xué),北京 100190)

3(中國互聯(lián)網(wǎng)絡(luò)信息中心,北京 100190)

近些年,基因組測序技術(shù)飛速發(fā)展,其發(fā)展經(jīng)歷包含三個(gè)階段：第一代主要通過凝膠電泳和放射自顯影技術(shù)確定待測的DNA 序列,其單次測序讀長長、準(zhǔn)確率高,但測序過程緩慢、吞吐量低；第二代核心思想為邊合成邊測序,即通過捕捉新合成的末端的標(biāo)記來確定DNA 的序列,其特點(diǎn)是讀長短、測序速度快、吞吐量大；第三代是單分子測序,成本較高,且目前測序技術(shù)尚未成熟[1].高通量測序技術(shù)作為第二代測序技術(shù),滿足了生命科學(xué)諸多領(lǐng)域的大規(guī)模測序需求,尤其在MSI 探測領(lǐng)域中,對(duì)于腫瘤的早期診斷、化療敏感性判斷[2]以及高危人群的圈定[3]等研究發(fā)揮了重要作用.

另外,作為腫瘤遺傳不穩(wěn)定的敏感指標(biāo),MSI 探測在化學(xué)治療劑的選擇[4,5]、預(yù)后判斷[6,7]及家族性癌癥風(fēng)險(xiǎn)評(píng)估[8,9]等方面均具有重要的臨床意義.有效地利用MSI 探測結(jié)果可達(dá)到減少醫(yī)療盲目性,提高M(jìn)SI 患者的治愈率,實(shí)現(xiàn)精準(zhǔn)醫(yī)療的目標(biāo).目前,MSI 探測主要分為生物學(xué)實(shí)驗(yàn)方法和基于高通量測序數(shù)據(jù)的計(jì)算方法.通過分析高通量測序數(shù)據(jù)的MSI 探測方法,無需額外的臨床測試或患者樣本處理,大大減少由于實(shí)驗(yàn)方法依賴少量微衛(wèi)星位點(diǎn)的探測或錯(cuò)配修復(fù)蛋白是否表達(dá)的局限性帶來的誤差與人力檢測成本,因此被廣大研究者所采用.但是,隨著測序成本的急劇下降,測序通量呈指數(shù)提高,眾多探測算法軟件已無法滿足使用者在運(yùn)行速度、存儲(chǔ)空間及用戶體驗(yàn)等多方面的需求.基于計(jì)算的MSI 探測方法或軟件主要存在兩個(gè)限制條件：(1) 當(dāng)使用樣本的腫瘤與正常組織成對(duì)測序數(shù)據(jù)探測MSI 時(shí),需要正常組織測序數(shù)據(jù)作為輸入?yún)⒄眨?2) 僅使用腫瘤組織測序數(shù)據(jù)時(shí),使用者需要針對(duì)不同的癌種或不同的測序平臺(tái)得到的測序數(shù)據(jù),分別使用大量MSS 樣本的正常組織測序數(shù)據(jù)構(gòu)建基準(zhǔn)線,這不僅會(huì)在一定程度上造成特異性MSI 位點(diǎn)的丟失,而且給使用者帶來了大量的額外工作,對(duì)存儲(chǔ)空間和運(yùn)行時(shí)間要求也非常高,無形之中給使用者提出了一個(gè)難題.因此,構(gòu)建基于單腫瘤組織測序數(shù)據(jù)探測樣本MSI 狀態(tài)的模型十分必要.

本文主要分為四個(gè)部分.第1 部分,介紹已有的相關(guān)研究工作；第2 部分,在前期MSI 探測領(lǐng)域的研究基礎(chǔ)上,針對(duì)MSIsensor1.1[10]無法實(shí)現(xiàn)只基于樣本腫瘤組織測序數(shù)據(jù)探測MSI 狀態(tài)的問題,提出一種采用信息熵理論的MSI 探測模型；第3 部分,對(duì)本文提出的模型進(jìn)行性能評(píng)估測試；第4 部分,總結(jié)本文提出的模型的意義和展望未來MSI 探測研究方向.

1 相關(guān)工作

微衛(wèi)星,又稱為短串聯(lián)重復(fù)序列(short tandem repeats,STRs),由長度在1～6 堿基對(duì)(base pair,bp)之間,重復(fù)次數(shù)為10～60 次的重復(fù)短核苷酸序列組成,其在真核生物的基因組中廣泛存在[11].MSI 是指在DNA 復(fù)制期間,由于錯(cuò)配修復(fù)系統(tǒng)的缺陷(Mismatch Repair Deficiency,MMRD),導(dǎo)致重復(fù)短核苷酸序列的插入或缺失,出現(xiàn)微衛(wèi)星新等位基因的現(xiàn)象.根據(jù)微衛(wèi)星不穩(wěn)定程度,MSI 分為MSI-H (Microsatellite Instability High),MSI-L (Microsatellite Instability Low)以及MSS,由于MSI-L 與MSS 沒有顯著性差異,本研究將MSI-L 視為MSS 處理.

隨著下一代測序技術(shù)的發(fā)展,基于高通量測序數(shù)據(jù)的MSI 檢測方法及軟件逐漸涌現(xiàn).目前,無需輸入腫瘤與正常組織成對(duì)測序的方法和軟件有mSINGS[12]、Lu 等人提出的基于插入和缺失變異(Insertion and Deletion,Indel)的MSI 檢測方法[13]、MSIseq[14]和MSIpred[15].對(duì)BAM (Binary sequence Alignment/Map format)文件直接進(jìn)行分析的mSINGS 無需配對(duì)的正常組織測序數(shù)據(jù)做參照,但在對(duì)樣本進(jìn)行分析之前,需要足夠的MSS樣本正常組織測序數(shù)據(jù)建立各個(gè)位點(diǎn)的基準(zhǔn)線,根據(jù)3σ法則判斷各個(gè)位點(diǎn)的穩(wěn)定性,最后通過統(tǒng)計(jì)不穩(wěn)定位點(diǎn)的比例獲得樣本的MSI 狀態(tài).Lu 等人提出的方法使用PI/PD 作為判別樣本MSI 狀態(tài)的指標(biāo),并通過實(shí)驗(yàn)證實(shí)該方法的可行性.其中,PD 是指微衛(wèi)星區(qū)域Deletion 占所有Deletion 的比例,PI 是指微衛(wèi)星區(qū)域Insertion 占所有Insertion 的比例.MSIseq 基于決策樹算法的分類器在單核苷酸替代(Single Nucleotide Substitution,SNS)和Indel 短突變的高級(jí)別數(shù)據(jù)上進(jìn)行模型訓(xùn)練和測試,該方法在測試時(shí)僅使用微衛(wèi)星區(qū)域小片段堿基插入與刪除的比例這一特征.MSIpred 利用支持向量機(jī)對(duì)從突變數(shù)據(jù)中提取的12 個(gè)特征進(jìn)行模型訓(xùn)練和測試.MSIseq、MSIpred 等是從突變注釋格式(Mutation Annotation Format,MAF)中提取特征用于模型訓(xùn)練和測試,但MAF 是由腫瘤與配對(duì)正常組織組織的測序數(shù)據(jù)產(chǎn)生.因此,提供一個(gè)直接對(duì)腫瘤組織測序數(shù)據(jù)的BAM 文件進(jìn)行分析的檢測工具十分必要.

基于腫瘤與正常組織成對(duì)測序數(shù)據(jù)的MSIsensor 1.1和MANTIS[16]可以直接對(duì)BAM 文件進(jìn)行分析.MANTIS將樣本的腫瘤組織與配對(duì)的正常組織在每個(gè)微衛(wèi)星位點(diǎn)的等位基因分布看作兩個(gè)向量,定義這兩個(gè)向量的L1范數(shù)作為衡量該位點(diǎn)穩(wěn)定程度的度量,最后通過統(tǒng)計(jì)該樣本所有位點(diǎn)的L1范數(shù)的平均值獲得該樣本的MSI 分?jǐn)?shù).MSIsensor 通過卡方檢驗(yàn)比較腫瘤組織和配對(duì)的正常組織的相同微衛(wèi)星位點(diǎn)的等位基因分布是否顯著不同判斷該位點(diǎn)是否穩(wěn)定,最后統(tǒng)計(jì)不穩(wěn)定的位點(diǎn)占總位點(diǎn)的比例判斷樣本的MSI 狀態(tài)[17,18].

通過對(duì)MSIsensor1.1和MANTIS 軟件之間的對(duì)比,本研究發(fā)現(xiàn)MSIsensor1.1 具有以下兩方面優(yōu)點(diǎn)：(1)設(shè)計(jì)基于C++實(shí)現(xiàn)、計(jì)算效率高；(2) 通過使用全基因組范圍內(nèi)的微衛(wèi)星位點(diǎn)信息的探測和卡方檢驗(yàn)方法判斷微衛(wèi)星狀態(tài),準(zhǔn)確性非常高.因此,本研究選擇該軟件進(jìn)行優(yōu)化.

MSIsensor1.1 基于卡方檢驗(yàn)的探測模型包含以下兩個(gè)功能模塊：

(1)獲取位點(diǎn)基本信息模塊.通過對(duì)參考基因組文件(Reference Genome.fa)執(zhí)行掃描(scan)命令,構(gòu)建一個(gè)由微衛(wèi)星位點(diǎn)基本信息組成的序列字典(microsatellites.list),其基本信息包括：微衛(wèi)星位點(diǎn)所在染色體號(hào)及其在染色體上的起始位置、微衛(wèi)星重復(fù)單元的大小及重復(fù)次數(shù)以及微衛(wèi)星位點(diǎn)的左右翼信息；

(2)探測樣本MSI 狀態(tài)模塊.通過對(duì)腫瘤組織測序數(shù)據(jù)文件(Sample.tumor.bam)、正常組織測序數(shù)據(jù)文件(Sample.normal.bam)和microsatellites.list 執(zhí)行msi 命令獲得待測樣本的MSI 分?jǐn)?shù)(MSIscore).該模塊共有四個(gè)結(jié)果文件,分別為：MSI 分?jǐn)?shù)文件(Output 文件)、微衛(wèi)星位點(diǎn)等位基因分布文件(Output_dis 文件)、germline 位點(diǎn)文件(Output_germline 文件)和somatic位點(diǎn)文件(Output_somatic 文件).該模塊具體步驟如下：① 計(jì)算每個(gè)位點(diǎn)在樣本腫瘤與正常組織成對(duì)測序數(shù)據(jù)上微衛(wèi)星的長度分布,并篩選在腫瘤和正常組織測序數(shù)據(jù)中測序深度均大于20 的微衛(wèi)星位點(diǎn)；② 利用卡方檢驗(yàn)比較符合條件①的微衛(wèi)星位點(diǎn)在腫瘤與正常組織中的分布差異,若顯著不同,則判定該位點(diǎn)為不穩(wěn)定的微衛(wèi)星位點(diǎn)；③ 通過統(tǒng)計(jì)不穩(wěn)定位點(diǎn)占總位點(diǎn)的比例獲得樣本的MSIscore,進(jìn)而判斷該樣本的MSI狀態(tài).

2 實(shí)驗(yàn)方法

目前,基于腫瘤與正常組織成對(duì)測序數(shù)據(jù)的MSI探測方法,利用配對(duì)的正常組織測序數(shù)據(jù)為參考,對(duì)腫瘤組織的微衛(wèi)星位點(diǎn)的不穩(wěn)定情況進(jìn)行評(píng)判,最終達(dá)到判斷樣本MSI 狀態(tài)的目標(biāo).對(duì)于缺少配對(duì)的正常組織測序數(shù)據(jù)做參照的情況,為了避免使用大量MSS 樣本的正常組織測序數(shù)據(jù)構(gòu)建基準(zhǔn)線而丟失特異性MSI 位點(diǎn),本文提出一種新模型,該模型根據(jù)腫瘤組織每個(gè)微衛(wèi)星位點(diǎn)等位基因分布的混亂情況,來判斷該位點(diǎn)及樣本的MSI 狀態(tài).

在正常細(xì)胞中,錯(cuò)配修復(fù)系統(tǒng)(Mismatch Repair,MMR)提供了一種高效的機(jī)制來糾正DNA 復(fù)制過程中發(fā)生的錯(cuò)誤.當(dāng)MMR 受損時(shí),例如錯(cuò)配修復(fù)基因MLH1、MSH2和MSH3 的失活,將導(dǎo)致基因序列——特別是微衛(wèi)星位點(diǎn)重復(fù)單元的錯(cuò)誤插入或刪除得不到糾正,微衛(wèi)星位點(diǎn)的重復(fù)單元的重復(fù)次數(shù)出現(xiàn)波動(dòng),即MSI 發(fā)生[7].從酵母的研究中發(fā)現(xiàn),在DNA 錯(cuò)配修復(fù)蛋白突變后,微衛(wèi)星中短缺失的數(shù)量明顯增加,而微衛(wèi)星中插入的數(shù)量沒有明顯變化[8].后續(xù)研究也證實(shí)了,由于MMRD,微衛(wèi)星區(qū)域會(huì)發(fā)生小片段堿基的插入與缺失.受已有研究成果的啟發(fā),我們考慮到,與MSS 樣本相比,MSI 樣本腫瘤組織的微衛(wèi)星位點(diǎn)的重復(fù)單元的重復(fù)次數(shù)的分布會(huì)因?yàn)镸MRD 呈現(xiàn)比較混亂的狀態(tài).圖1為MSS 樣本與MSI 樣本的腫瘤組織在同一微衛(wèi)星位點(diǎn)的等位基因分布.從圖中可以看出,MSI 樣本的微衛(wèi)星位點(diǎn)等位基因分布與MSS 樣本相比較混亂.鑒于MSI 樣本腫瘤組織微衛(wèi)星位點(diǎn)測序數(shù)據(jù)這一特性,本研究提出使用信息熵理論來探測樣本的MSI 狀態(tài).

圖1 MSI 樣本與MSS 樣本的同一微衛(wèi)星位點(diǎn)等位基因分布比較

信息熵是描述“混亂”程度的量度.系統(tǒng)越有序,信息數(shù)據(jù)越集中的地方,熵值越小；系統(tǒng)越混亂,信息數(shù)據(jù)越分散的地方,熵值越大[19].擴(kuò)展后的MSIsensor 工作流程如圖2所示.

圖2 MSIsensor 工作流程圖

圖中虛線所指是MSIsensor1.1 已有模塊,實(shí)線所指功能是本文在實(shí)現(xiàn)探測樣本MSI 狀態(tài)功能中提出的基于信息熵理論在樣本單腫瘤組織測序數(shù)據(jù)進(jìn)行探測的新模型.該模型共有三個(gè)輸出,分別為：Output 文件、Output_dis 文件、Output_somatic 文件.新模型主要包含以下子模塊：

(1)目標(biāo)位點(diǎn)篩選功能.通過計(jì)算每個(gè)微衛(wèi)星位點(diǎn)在腫瘤組織測序數(shù)據(jù)上的分布,可篩選出測序深度大于20(默認(rèn))的微衛(wèi)星位點(diǎn).

(2)目標(biāo)位點(diǎn)過濾功能.由于測序過程中可能會(huì)帶來背景噪聲信息,因此,對(duì)(1)篩選出的目標(biāo)位點(diǎn)集進(jìn)行過濾.過濾規(guī)則為：過濾掉腫瘤組織在該位點(diǎn)支持reads 數(shù)小于3 的等位基因(被認(rèn)為是噪聲數(shù)據(jù)),接下來以剩下的等位基因數(shù)作為樣本在該位點(diǎn)的等位基因數(shù).

(3)加權(quán)信息熵值計(jì)算功能.根據(jù)步驟(2)得到的每個(gè)位點(diǎn)等位基因分布情況,計(jì)算加權(quán)信息熵值.信息熵值計(jì)算公式如下：

其中,pi指的是針對(duì)位點(diǎn)的每個(gè)重復(fù)次數(shù),其不為零支持reads 數(shù)占該位點(diǎn)支持reads 數(shù)總和的比例.支持reads 數(shù)不為零的重復(fù)次數(shù)的個(gè)數(shù)會(huì)較大程度影響每個(gè)位點(diǎn)的信息熵值,即每個(gè)位點(diǎn)的信息熵值會(huì)因支持reads 數(shù)不為零的重復(fù)次數(shù)的個(gè)數(shù)變化而跨越幅度較大.為了方便后續(xù)信息熵值閾值(cutoff)的確定,對(duì)每個(gè)位點(diǎn)的熵值賦予一定的權(quán)重(weight).若該位點(diǎn)的H(U)*weight大于等于0.3(默認(rèn)),則判定該位點(diǎn)不穩(wěn)定.其中,weight等于該位點(diǎn)支持reads 數(shù)不為零的重復(fù)次數(shù)的個(gè)數(shù)的倒數(shù).在TCGA 數(shù)據(jù)集的不同癌種上使用不同的閾值進(jìn)行多次實(shí)驗(yàn),通過比較得到,當(dāng)使用閾值0.3 時(shí),在不同癌種上Accuracy 都比較高.圖3給出了在TCGA 數(shù)據(jù)集上隨著cutoff變化時(shí)的Accuracy 趨勢.

圖3 不同cutoff 下的Accuracy 比較

(4) MSIscore 計(jì)算功能.通過統(tǒng)計(jì)由(3)得到不穩(wěn)定位點(diǎn)數(shù)占滿足條件(1)的微衛(wèi)星總位點(diǎn)數(shù)的比例,獲得該樣本的MSIscore,進(jìn)而判斷樣本的MSI 狀態(tài).

3 實(shí)驗(yàn)結(jié)果與討論

本研究使用的數(shù)據(jù)來自國際癌癥基因組圖譜(The Cancer Genome Atlas,TCGA)計(jì)劃[20]、歐洲基因組-表型組檔案(European Genome-phenome Archive,EGA)[21,22]和北京腫瘤醫(yī)院(Beijing Cancer Hospital),包含的癌癥類型有胃癌(STomach ADenocarcinoma,STAD)、子宮內(nèi)膜癌(Uterine Corpus Endometrial Carcinoma,UCEC)、結(jié)腸癌(COlon ADenocarcinoma,COAD)和直腸腺癌(REctum ADenocarcinoma,READ),具體情況如表1所示.

表1 TCGA、EGA和Beijing Cancer Hospital 數(shù)據(jù)集

為了衡量MSIsensor 擴(kuò)增的使用單腫瘤組織測序數(shù)據(jù)探測MSI 模塊的表現(xiàn)性能,本研究將其與MSI sensor 基于卡方檢驗(yàn)的MSI 探測模塊在表1所示數(shù)據(jù)上的測試結(jié)果進(jìn)行對(duì)比.由于不同的癌癥類型或不同測序平臺(tái)得到的測序序列會(huì)稍有不同,本文不建議使用統(tǒng)一的MSIscore 閾值劃分MSS 樣本和MSI 樣本.針對(duì)基于卡方檢驗(yàn)在腫瘤與正常組織成對(duì)測序數(shù)據(jù)上探測MSI 的模塊,本研究在TCGA STAD、TCGA CRC、TCGA UCEC、EGA STAD、EGA COAD和Beijing Cancer Hospital 數(shù)據(jù)集上采用的MSIscore 閾值分別為3、10、3、3、10、和15,其各項(xiàng)性能指標(biāo)表現(xiàn)如表2所示,其中,準(zhǔn)確率(Accuracy)表示預(yù)測狀態(tài)正確的樣本數(shù)占待檢測樣本總數(shù)的百分比；精確率(Precision)表示軟件預(yù)測的真實(shí)MSI 樣本數(shù)占軟件預(yù)測出的MSI 樣本總數(shù)的百分比；靈敏度(Sensitivity)表示軟件預(yù)測的真實(shí)MSI 樣本數(shù)占MSI 總樣本數(shù)的百分比；特異性(Specificity)表示軟件預(yù)測的真實(shí)MSS樣本數(shù)占MSS 總樣本數(shù)的百分比；F-分?jǐn)?shù)(F-score)表示Sensitivity和Precision 兩個(gè)性能度量的調(diào)和平均數(shù),用來綜合評(píng)估軟件的性能.具體分類效果如圖4所示.針對(duì)基于信息熵在單腫瘤組織上探測MSI 的模塊,本研究在TCGA STAD、TCGA CRC、TCGA UCEC、EGA STAD、EGA COAD和Beijing Cancer Hospital數(shù)據(jù)集上采用的閾值分別為13、20、13.5、11、15和20.該模型的各項(xiàng)性能指標(biāo)表現(xiàn)如表3所示,具體分類效果如圖5所示.

表2 基于卡方檢驗(yàn)探測模塊的測試結(jié)果

以上結(jié)果證明,在TCGA 上確定的對(duì)每個(gè)位點(diǎn)判定是否穩(wěn)定的加權(quán)信息熵閾值0.3,在EGA、Beijing Cancer Hospital 數(shù)據(jù)集上同樣適用.同時(shí),該結(jié)果也表明本文提出的基于信息熵理論在單腫瘤組織測序數(shù)據(jù)上探測樣本MSI 狀態(tài)的方法是可行的.

基于信息熵的MSI 探測方法只需要對(duì)單腫瘤組織測序數(shù)據(jù)進(jìn)行分析,因此,在探測MSI 的前期準(zhǔn)備工作中不需要對(duì)正常組織進(jìn)行測序,可以節(jié)省測序成本；且與MSIsensor1.1 在腫瘤與正常組織測序數(shù)據(jù)上探測MSI 的方法相比較,其運(yùn)行效率提高了一倍.

4 結(jié)論與展望

本研究針對(duì)樣本單腫瘤組織的外顯子或全基因組測序數(shù)據(jù)的MSI 探測問題,選擇從高通量測序數(shù)據(jù)入手,基于已有軟件MSIsensor1.1,增加采用信息熵理論探測MSI 狀態(tài)的模塊.擴(kuò)增后的軟件不僅可針對(duì)腫瘤與正常組織成對(duì)測序數(shù)據(jù)探測樣本MSI 狀態(tài),而且當(dāng)無配對(duì)的正常組織測序數(shù)據(jù)做輸入?yún)⒄瘴锖筒粐?yán)重影響精度的情況下,可滿足樣本單腫瘤組織測序數(shù)據(jù)的MSI 探測需求,是一套功能完備、兼容多種測序數(shù)據(jù)模式、具備實(shí)用性和前瞻性的MSI 探測軟件.

后續(xù)我們將圍繞以下兩方面工作展開：(1)基于此工作深入挖掘每個(gè)位點(diǎn)信息并篩選與腫瘤發(fā)生發(fā)展相關(guān)的基因,這樣不但可以節(jié)省測序的成本,而且可以提高檢測的正確率；(2)探究基于血漿的MSI 探測方法.癌細(xì)胞具有擴(kuò)散性,當(dāng)癌細(xì)胞侵犯血管進(jìn)入血液時(shí),血液中也會(huì)有MSI 信號(hào),但此時(shí)信號(hào)很微弱,不易檢測.目前大多數(shù)MSI 檢測局限于組織測序數(shù)據(jù),但是對(duì)于不方便提取組織的患者,開發(fā)基于血漿的MSI 探測方法具有重要的意義.下一步,我們將圍繞此工作展開,探求更強(qiáng)、更符合實(shí)際應(yīng)用需求的MSI 探測軟件.

圖4 基于卡方檢驗(yàn)探測模塊的測試結(jié)果

表3 基于信息熵理論探測模塊的測試結(jié)果

圖5 基于信息熵理論探測模塊的測試結(jié)果