腫瘤突變印跡的求解算法比較①

何小雨, 韓鑫胤, 牛北方

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

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

癌細(xì)胞的基因組攜帶體細(xì)胞突變(somatic mutation).體細(xì)胞突變指發(fā)生在除生殖細(xì)胞之外的人體細(xì)胞中的DNA結(jié)構(gòu)上堿基對的改變, 其過程貫穿個體的整個生命周期. 體細(xì)胞突變的誘因主要有物理因素(如射線)、化學(xué)因素(如抗生素、煙草)以及生物因素(如細(xì)菌或病毒基因的融合)[1]. 如圖1所示, 從受精卵第一次分裂開始, 細(xì)胞受到某種致癌因子的刺激而發(fā)生突變.突變過程可能會引起DNA損傷和修復(fù), 并產(chǎn)生單核苷酸替換(single nucleotide substitutions)、短插入和刪除(short insertions and deletions)、結(jié)構(gòu)變異(structure variants)和染色體拷貝數(shù)變化(copy number variation)等[2]. 各致癌因子之間相互獨立, 在體細(xì)胞基因組上留下獨特的記號, 隨著個體生命的延續(xù), 體細(xì)胞突變不斷疊加累積, 最終使細(xì)胞脫離控制, 形成生長增殖不受限制的癌細(xì)胞. 生物醫(yī)學(xué)研究中普遍認(rèn)為, 細(xì)胞在正常生長分裂過程中產(chǎn)生的突變在基因組上是隨機分布的,而由特定致癌因子導(dǎo)致的突變則具有一定的模式. 這種模式能夠反映細(xì)胞癌變過程中曾暴露于哪些致癌因子以及在其中的暴露程度. 研究中將基因組中由致癌因子引起的這種特有突變模式叫做突變印跡(mutational signature)[3].

圖1 突變過程在癌癥基因組中留下的特征性印記

早在DNA雙螺旋結(jié)構(gòu)被發(fā)現(xiàn)之前, 已有臨床研究發(fā)現(xiàn)連續(xù)的紫外線輻射會加速細(xì)胞增殖的相對速度, 從而明確了紫外線的過量輻射是皮膚癌的一大誘因[4]. 測序技術(shù)的出現(xiàn)尤其是下一代測序技術(shù)(Next Generation Sequencing, NGS)的蓬勃發(fā)展, 迅速推動了分子生物學(xué)和生物信息學(xué)的研究進展[5]. 研究顯示, 人的基因組序列中存在大約400萬個突變位點, 全外顯子組序列中可以檢測到大約6萬-8萬個突變位點[6].大多數(shù)癌癥中的體細(xì)胞突變屬于“乘客”突變(passenger mutations), 即不導(dǎo)致細(xì)胞癌變, 在癌癥發(fā)展的過程中也不會被正向選擇. 僅有少數(shù)的突變是“驅(qū)動”突變(driver mutations), 即賦予細(xì)胞生長優(yōu)勢并且?guī)椭┘?xì)胞增殖[7].

長久以來, 對體細(xì)胞突變的研究方向局限在驅(qū)動基因(driver genes)上, 如TP53[8]. 然而, “乘客”事件可能不是癌癥發(fā)展的原因, 卻一定是細(xì)胞癌變過程中的產(chǎn)物. 因此“乘客”突變攜帶并記錄了DNA損傷和修復(fù)的豐富歷史證據(jù), 這些證據(jù)在癌癥的發(fā)生發(fā)展中起關(guān)鍵作用. 因此, 探尋突變過程時應(yīng)充分考慮“乘客”突變.突變印跡分析將發(fā)生在基因組上的所有突變位點納入考慮, 是對癌癥基因組研究的重要補充.

突變印跡的概念由Alexandrov等人提出的[3], 他們在21個乳腺癌的全基因組單堿基替換突變類型中計算并提取了5個特殊的突變模式, 這些單基因替換突變并未進行“驅(qū)動”和“乘客”的區(qū)別, 而是將他們?nèi)孔鳛榉治鰧ο? 研究中還嘗試解釋了其生物學(xué)機制, 如以TpCpX三核苷酸處的C＞T, C＞G和C＞A替換為主要特征的突變印跡是由5-甲基胞嘧啶作用自發(fā)的內(nèi)源性突變過程導(dǎo)致的, 是癌癥樣本中絕大多數(shù)的突變類型, 并且以這種顯性的突變形式存在于約10%的雌激素受體(Estrogen Receptor, ER)陽性乳腺癌患者中[2].雖然當(dāng)時并未對計算出的5個印跡進行完全的生物學(xué)解釋, 但該研究已經(jīng)對潛在的突變機制形成了一些見解并提供了研究的思路. 2013年, Alexandrov等人又開展了一項涉及30個癌種, 7042個樣本的研究, 在4 938 362個單堿基替換中解析出21個經(jīng)過生物學(xué)驗證的突變印跡[3]. 2018年, 該研究團隊將研究數(shù)據(jù)集繼續(xù)擴大到23 829個癌癥樣本的4729 690個體細(xì)胞突變集合上, 同時納入了更多的突變類型: 單堿基替換, 雙堿基替換, 短插入刪除和結(jié)構(gòu)變異. 更大規(guī)模的數(shù)據(jù)集中發(fā)現(xiàn)了新的突變印跡: 49種單堿基替換印跡, 11種雙堿基替換印跡, 17個短插入刪除印跡以及11個結(jié)構(gòu)變異的聚類[9].這些印跡既驗證了紫外線、煙草、酒精等與癌癥發(fā)生相關(guān)聯(lián)的外部因素, 一些新的印跡也被證實與獨特的臨床特征相關(guān), 這說明突變印跡的分析很可能成為靶向治療中新的潛在生物標(biāo)志物[10].

雖然一些突變印跡反映的致癌因子(如紫外線照射、煙草)可以通過跟蹤調(diào)查和統(tǒng)計來識別, 但是一種正式的數(shù)學(xué)方法可以提取人類難以察覺的更微妙的元素和較弱的信號, 同時還需要評估該元素或信號在癌變過程中的比例. 并且隨著基因組計劃的完成, 大規(guī)模的基因組測序數(shù)據(jù)也對突變印跡問題的求解提出了考驗. 目前, 已有多個基于NGS的識別突變印跡的算法和軟件, 但對該問題的求解在參數(shù)設(shè)置等方面尚未達成共識. 據(jù)調(diào)研, 目前缺乏關(guān)于體細(xì)胞突變印跡分析算法和軟件的比較, 并且隨著更多癌癥突變位點的檢出,突變印跡分析將獲得很大的挖掘空間. 因此, 開展該領(lǐng)域的相關(guān)研究, 清晰詳細(xì)地介紹和討論是十分必要的.

總體來講, 本文的主要研究內(nèi)容包括:

(1) 突變印跡問題及其數(shù)學(xué)模型闡述.

(2) 突變印跡提取算法及評價.

(3) 突變印跡提取軟件的其他功能介紹.

(4) 比較總結(jié)實驗結(jié)果并提出新的解決方案.

1 突變印跡的數(shù)據(jù)模型

突變印跡問題可以描述為: 從復(fù)雜的體細(xì)胞突變信號中尋找獨立致癌因子使基因組發(fā)生的特有的改變模式. 癌細(xì)胞基因組的改變可能是多個致癌因子引發(fā)突變的累積, 原始致癌因子使基因組發(fā)生改變的程度也不相同. 因此, 可將突變印跡問題抽象為盲源分離問題(Blind Source Separation, BSB)[11].

盲源分離問題是研究在系統(tǒng)的傳遞函數(shù)、源信號的混合系數(shù)及概率分布未知的情況下, 利用源信號之間相互獨立這一微弱的已知條件, 如何從一組復(fù)雜的混合信號中分離出獨立的不可觀測的源信號. 盲源分離作為陣列信號處理的一種新技術(shù), 允許有意義地學(xué)習(xí)對象的不同部分, 近幾年來受到廣泛關(guān)注.

盲源分離問題在突變印跡分析中的應(yīng)用可以表述為算法1.

算法1. 突變印跡分離算法1. 計算體細(xì)胞突變信號中的最優(yōu)突變印跡的組合, 以表示在癌癥的發(fā)展過程中每個獨立突變過程的累積；2. 計算每個印跡在每個獨立癌癥基因組的體細(xì)胞突變中的比例, 表示印跡對應(yīng)的致癌因子對癌變過程的貢獻度.

算法中原始體細(xì)胞突變信號常采用結(jié)合上下文的三堿基結(jié)構(gòu)表示. 基因組是由腺嘌呤(A)、胸腺嘧啶(T)、鳥嘌呤(G)、胞嘧啶(C)組成的序列. 由于基因組正負(fù)鏈上的堿基以互補配對原則, 即A與T配對, G與C配對形成堿基對, 則基因組上可能出現(xiàn)的單核苷酸替換情況有C24種(圖2), 其中:

(1) C＞A: 代表C＞A和G＞T兩種單核苷酸替換方式；

(2) C＞G: 代表C＞G和G＞C兩種單核苷酸替換方式；

(3) C＞T: 代表C＞T和G＞A兩種單核苷酸替換方式；

(4) T＞A: 代表T＞A和A＞T兩種單核苷酸替換方式；

(5) T＞C: 代表T＞C和A＞G兩種單核苷酸替換方式；

(6) T＞G: 代表T＞G和A＞C兩種單核苷酸替換方式.

將這6種替換方式表示為如式(1)所示的單堿基替換字典V:

然后, 將其上游(5′端)、下游(3′端)各一個堿基作為其上下文(如圖2中標(biāo)識), 三堿基結(jié)構(gòu)表示為式(2):

圖2 從體細(xì)胞突變列表中構(gòu)建96種突變模式

顯然, 單堿基上下文結(jié)構(gòu)的表示方法有4×C24×4=96種. 則單個基因組上的體細(xì)胞突變譜可表示為如式(3)所示的向量:

其中,Tkn表示與Tn關(guān)聯(lián)的突變過程(第n個致癌因子)引起的第k個突變符號的頻率. 因此:

將一個癌癥隊列中每個患者的體細(xì)胞突變表示為如式(3)所示向量, 則該隊列的突變目錄即可表示為:

其中, 突變目錄的每個元素可以近似地認(rèn)為是使正常細(xì)胞發(fā)展為腫瘤細(xì)胞的潛在突變過程的特征的線性疊加, 且每個特征通過在相應(yīng)過程中的暴露程度來加權(quán),如式(6):

因此突變印跡可通過將M表示為兩個較小的矩陣實現(xiàn). 即:

其中,M是研究隊列中g(shù)個基因組的96個突變類型的突變頻率:

而Sk×n是該隊列基因組由n個致癌因子導(dǎo)致的特定的突變模式:

En×g表示g個基因組在n個致癌因子中的暴露程度:

此外, 突變印跡因其代表的自然意義而具有兩個特征: (1)研究對象的非負(fù)性. 一個癌癥隊列中所有基因組體細(xì)胞突變目錄M表示的是每個樣本的每種突變類型的突變頻率, 分解矩陣S表示突變過程的特征,系數(shù)矩陣E表示突變過程在基因組上的強度； (2)研究目標(biāo)是期望從體細(xì)胞突變目錄中提取具有生物學(xué)意義的特征, 即致癌因子在基因組上留下的特殊記號.

2 算法分類和比較

根據(jù)對式(7)的求解方法, 可以將突變印跡問題分為3類: 一是非負(fù)矩陣分解(Nonnegative Matrix Factorization, NMF)的方法； 二是期望最大化(Expectation Maximization, EM)方法； 三是線性回歸(Linear Regression, LR)方法. 按照是否能夠發(fā)現(xiàn)新印跡可以將突變印跡問題分為允許新印跡發(fā)現(xiàn)的“提取”方法和對已知印跡的“擬合”方法.

2.1 求解方法的比較

2.1.1 NMF方法

NMF算法是由Lee和Seung于1999年提出的一種矩陣分解方法, 它使分解后的所有分量均為非負(fù)值,并且同時實現(xiàn)非線性的維數(shù)約減[12]. NMF作為從各種類型的高維生物數(shù)據(jù)中提取有意義成分的一種強大技術(shù)屢次脫穎而出. 此外在其他領(lǐng)域也有NMF方法的成功應(yīng)用, 如“雞尾酒會”問題, 人臉識別問題等[13]. 由于上述數(shù)學(xué)模型中矩陣的內(nèi)在非負(fù)性, 使得NMF特別適合于突變特征推斷問題. NMF也是第一個用來嘗試分析突變印跡問題的算法. 如表1所示, 目前已有如SigProfiler、SomaticSignatures、sigfit和Mutational-Patterns等多個軟件基于NMF算法來解決印跡分解問題.

表1 突變印跡分析軟件匯總

SigProfiler方法通過找到矩陣S和E來準(zhǔn)確地提取N個突變印跡, 同時解決由式(7)導(dǎo)出的非凸優(yōu)化問題. 該方法選擇矩陣范數(shù)作為Frobenius重構(gòu)誤差:

以SigProfiler為代表的NMF算法具體步驟如算法2.

算法2. 以SigProfiler為代表的NMF算法1. 初始化隨機非負(fù)矩陣S, E；2. 將初始突變目錄矩陣降維, 將所有突變類型中占比≤1%的突變類型刪除, 得到矩陣M′；3. 迭代:(1)對矩陣M′進行蒙特卡洛自舉重采樣(Monte Carlo bootstrap resampling), 得到矩陣M"；(2)乘法更新算法應(yīng)用于M", 得到使式(11)中的Frobenius范數(shù)最小的S和E；4. 對S劃分聚類, 得到N個簇；5. 將每個簇中的s歸一化, 得到S的N個向量；6. 求暴露矩陣E.

基于NMF方法提取突變印跡的軟件需要指定分解個數(shù)N作為程序輸入, SigProfiler根據(jù)大量實驗給出的建議值是:

在實際操作中, SigProfiler根據(jù)每個N值計算模型的總體再現(xiàn)性和Frobenius范數(shù)誤差. 最終需要人工干預(yù)選擇N值, 使得分解的印跡矩陣S具有高度的再現(xiàn)性同時顯示出低的總體重建誤差.

2.1.2 EM 方法

Fischer等人基于NMF方法從基因組內(nèi)在特性出發(fā), 應(yīng)用概率模型解決突變印跡問題[14]. 基因組的內(nèi)在的特性(如CpG雙核苷酸的不均一分布, 拷貝數(shù)變化)會影響三堿基序列結(jié)構(gòu)發(fā)生突變的可能性, 繼而使模型在推斷的突變模式上產(chǎn)生偏倚. 理論上, 使用概率模型可在求解過程中充分考慮突變發(fā)生的可能性, 能夠更準(zhǔn)確地分離出真實突變過程的印跡.

將三堿基序列結(jié)構(gòu)的突變可能性表示為非0的k元組, 其中Okg表示基因組g上第k個突變類型發(fā)生突變的可能性:

EM算法的目標(biāo)是使包含隱變量的數(shù)據(jù)集的后驗概率或似然函數(shù)最大化, 進而得到最優(yōu)的參數(shù)估計. 文獻[14]將突變印跡分解問題重構(gòu)為一個如式(14)所示的概率模型. 其中突變目錄矩陣(M)分布為獨立的泊松隨機變量, 其元素由印跡矩陣(S)與暴露矩陣(E)乘積確定, 通過期望最大化算法對S和E進行估計.

算法的具體執(zhí)行過程如算法3.

算法3. 以EMu為代表的EM算法1. 猜想模型參數(shù)S(0)；2. 迭代?E(1)給定當(dāng)前猜測參數(shù)S(k), 得到曝光估計值 ；?E (2)使用 更新下一次迭代的參數(shù)估計值S(k+1)；(3)當(dāng)P(M|S)收斂到極大值時, 迭代結(jié)束；3. 比較不同N下的數(shù)據(jù)可能性, 確定突變過程數(shù)量.

值得注意的是, 雖然EMu是建立在對NMF的有效替代解釋的基礎(chǔ)上, 該解釋將NMF視為對特定問題的EM應(yīng)用, 但EMu的新概念和優(yōu)點并不是固有的EM范例特性, 也可以通過其他方法進行同化的顯式增強. 另一方面, EMu對初始條件的敏感度與常規(guī)NMF相同. 盡管如此, EMu成功地利用了突變印跡推斷的概率形式來解決以前未曾探索過的方向, 即結(jié)合了基因組的內(nèi)在特性和腫瘤特定的突變可能性并確定了突變印跡的個數(shù).

2.1.3 LR方法

線性回歸(Linear Regression, LR)指通過對大量的觀測數(shù)據(jù)進行處理, 從而得到比較符合事物內(nèi)部規(guī)律的數(shù)學(xué)表達式[15]. 在NMF方法發(fā)現(xiàn)了一些可解釋的突變印跡的基礎(chǔ)上, 印跡提取問題可擴展為新的描述: 對某個癌癥隊列的研究不再需要發(fā)現(xiàn)新的印跡, 而重點在于得到腫瘤中存在的、可解釋的印跡以及其對癌癥發(fā)展的作用程度. 因此, 突變印跡提取轉(zhuǎn)變?yōu)閷σ阎蛔冇≯E分布的擬合, 即對線性方程(15)的求解.

由于線性回歸方法對先驗知識的依賴, 尤其依賴對已知突變印跡的個數(shù)選擇和組合, 導(dǎo)致其在實際應(yīng)用中十分受限, 如目前已知的突變印跡集尚不能完全地解釋癌變過程, 個別突變印跡沒有得到完備的生物學(xué)解釋等.同時研究中也指出該方法的準(zhǔn)確性偏低, Maura等人在2019年的研究中發(fā)現(xiàn), 使用同樣的突變目錄(M)情況下, 非負(fù)矩陣分解算法提取的突變印跡大部分能夠被驗證其所代表的生物學(xué)意義, 而使用線性回歸方法擬合的突變印跡大多為“無特征”或“平坦”印跡, 即6種突變模式的頻率分布相對均勻, 無明顯差異[16].

2.2 “提取”算法和“擬合”算法

表1中的軟件和算法分別從不同角度出發(fā)來解決印跡問題, 根據(jù)算法特征及必需的輸入數(shù)據(jù), 可將其分為印跡“提取”算法和“擬合”算法. 基于NMF的“提取”算法以突變目錄矩陣M和分解個數(shù)N為輸入, 求解印跡和暴露矩陣, 因此分解出的印跡矩陣S中可能會出現(xiàn)新的印跡. 與提取方法不同, 擬合方法以目錄矩陣M和已知的印跡矩陣S為輸入, 將M中潛在的印跡擬合為S的線性表達.

“提取”的方法的優(yōu)勢是不依賴先驗知識(已知的突變印跡), 同時允許提取出新的突變印跡. 該算法也存在局限性. 首先, 同時發(fā)生的多個獨立的突變印跡可能會被合并為一個印跡. 其次, 對于非常復(fù)雜的印跡可能會因其較小的貢獻度而拆分為兩個或多個印跡.

反之, “擬合”的方法依賴大量先驗知識, 如分析的癌種有哪些致癌因子, 這些致癌因子的突變印跡分別是什么等. 現(xiàn)有公開發(fā)布的突變印跡總計有81種, 且50%尚未得到生物學(xué)解釋, 突變目錄矩陣本身的大小以及疾病類型乃至分型都會對印跡個數(shù)和組合的選擇產(chǎn)生影響. 由于其對已知突變印跡的依賴, 因而不能發(fā)現(xiàn)新的印跡. 此外, 當(dāng)主觀性較強的先驗知識被輸入時,很可能導(dǎo)致過擬合現(xiàn)象, 即夸大某個印跡在該癌種發(fā)生發(fā)展中的權(quán)重. 反之, 拋棄先驗知識的限制, 將全部已知的突變印跡作為輸入, 則會導(dǎo)致特異性突變印跡的滲透, 即少量樣本的突變印跡分配至整個隊列的樣本中或擬合到并未在癌種的發(fā)生發(fā)展中起作用的印跡上.

同時, 兩類算法也存在共性問題, 當(dāng)不同的突變印跡的組合可解釋同一個突變目錄矩陣時, 印跡提取就會變得不明確； 當(dāng)隊列中存在少量異質(zhì)性較高樣本時,其突變印跡因貢獻度不高而被過濾掉, 從而掩蓋隊列的異質(zhì)性.

3 突變印跡分析軟件的其他功能

在實際應(yīng)用中, 除了對式(7)的求解外, 突變目錄矩陣的構(gòu)建、突變印跡及暴露矩陣的可視化、運行環(huán)境等問題在生物信息學(xué)領(lǐng)域也是廣受關(guān)注的需求.

獲取突變上下文是構(gòu)建突變目錄矩陣的關(guān)鍵. 基于NGS的突變識別軟件如VarScan2[17]、Strelka2[18]、SomaticSniper[19]等給出的突變并不包含上下文信息.有些軟件如VarDict[20]在INFO列使用“LSEQ”和“RSEQ”分別給出突變上下文, 但大多數(shù)的軟件不能提供此信息, 需要重新計算. 可解決的方案是利用突變位置從參考基因組中獲取上下文. 但是人類的參考基因組有超過30億個堿基, 分別計算一個基因組中上百萬個突變位點的上下文需要一定的時間消耗. 此外, 通過突變注釋軟件(如Oncotator[21])能夠獲得上下文堿基序列, 然而這就需要在計算突變印跡前先對突變列表進行注釋, 增加了構(gòu)建突變印跡分析流程的復(fù)雜性. 因此很多軟件在內(nèi)部增加了計算上下文的過程, Pedersen等人開發(fā)了對VCF快速處理的Python程序Cyvcf2,可實現(xiàn)快速的VCF文件處理[22].

突變印跡和暴露矩陣可視化也是基因組數(shù)據(jù)挖掘的重要內(nèi)容. 除了EMu外, 目前用于計算突變印跡的軟件都提供了模塊化的可視化方法. 這些可視化方法都直接使用各自的分析結(jié)果作為輸入, 在實際操作中無需分析人員具備專業(yè)的繪圖知識.

最后, 基因組數(shù)據(jù)挖掘軟件大多依賴Linux環(huán)境,需要操作人員具備在Linux環(huán)境下編譯、安裝和運行軟件的能力. 因此, 除了發(fā)布軟件包, 多個軟件如PMSignature、MutaGene等提供了門戶網(wǎng)站的計算方式, 這在很大程度上方便了缺少計算機知識的生物醫(yī)學(xué)研究人員, 同時這也是整個生物信息學(xué)領(lǐng)域軟件的發(fā)展趨勢.

4 實驗分析

4.1 數(shù)據(jù)說明

本文選擇了21個宮頸癌的全基因組測序數(shù)據(jù), 分別使用基于非負(fù)矩陣分解算法的軟件SigProfiler[3]和基于線性回歸方法的軟件Mutalisk[23]分析突變印跡. 由于基于期望最大化的方法需要關(guān)鍵的先驗知識(基因組序列的復(fù)雜性), 本文實驗中不包括對EM算法的實驗.

4.2 實驗過程

SigProfiler的實驗過程如下:

(1) 構(gòu)建體細(xì)胞突變目錄矩陣；

(2) 將最大提取印跡數(shù)設(shè)置為10, 即N的最大取值為10；

(3) 運行SigProfiler；

(4) 選擇結(jié)果穩(wěn)定性高且殘差相對較小的突變印跡的個數(shù)的最大值N′；

(5) 提取N′個突變印跡, 構(gòu)成集合S；

(6) 計算S中每個印跡與公開發(fā)表印跡的相似性.

Mutalisk的實驗過程如下:

(1) 合并21個基因組的突變識別格式文件；

(2) 選擇COSMIC的30個突變印跡作為先驗輸入S；

(3) 運行Mutalisk；

(4) 選擇文獻中記載的與宮頸癌相關(guān)的印跡作為先驗輸入S′；

(5) 再次運行Mutalisk.

4.3 實驗結(jié)果

圖3表示了對21個癌癥患者的體細(xì)胞突變目錄矩陣進行NMF分解時, 不同的突變印個數(shù)對應(yīng)的突變印跡穩(wěn)定性(左)和原始突變目錄矩陣的重建誤差(右).

圖3 突變印跡數(shù)對應(yīng)的突變印跡穩(wěn)定性和重建誤差

實驗中, 將分解印跡數(shù)N設(shè)置為[1,10], 對每個N的取值分別計算兩個指標(biāo), 最終選擇再現(xiàn)性高且重建誤差小的印跡個數(shù). 通過圖3可以看出, 當(dāng)突變印跡數(shù)N＜7時, 隨著N的增大, 再現(xiàn)性逐漸降低(N=3)除外, 重構(gòu)誤差逐漸降低, 而在N＞7時, 再現(xiàn)性顯著降低,重構(gòu)誤差穩(wěn)定. 因此, 為了同時保證印跡矩陣的穩(wěn)定性高和原始矩陣的重構(gòu)誤差低, 選擇印跡個數(shù)為7. 圖4表示了具體7個突變印跡.

圖4 21個癌癥基因組中提取的7個突變印記

為了說明提取的7個印跡與該癌種發(fā)病之間的關(guān)系, 計算了7個突變印跡與國際公開數(shù)據(jù)庫COSMIC(Catalogue Of Somatic Mutations In Cancer)中2017年收錄的已知30種突變印跡的余弦相似性以及可能的生物學(xué)解釋(表2).

表2 SigProfiler提取突變印跡與已知印跡相似性

表3是Mutalisk采用30個COSMIC[24]的印跡為輸入時的結(jié)果, 其將21個宮頸癌基因組的突變目錄擬合到S5, S2, S13, S8, S1, S28六個已知印跡中.

當(dāng)選擇已知相關(guān)印跡(S1, S2, S5, S6, S8, S10, S13,S26)為先驗知識時, Mutalisk給出5種擬合結(jié)果. 其中選擇貝葉斯信息準(zhǔn)則(Bayesian Information Criterion,BIC)最小的擬合結(jié)果, 即S5, S2, S13, S1, S6 (表4).

表4 Mutalisk與已知6種印跡擬合情況及相似性

4.4 實驗總結(jié)

通過對比表2至表4中的余弦相似性可知, 基于非負(fù)矩陣分解的軟件SigProfiler發(fā)現(xiàn)的印跡與國際研究公開的印跡相似性均高于基于線性回歸的擬合方法.雖然SigProfiler提取的2, 6, 7并未得到相應(yīng)的生物學(xué)解釋, 但其余4個印跡分別與APOBEC基因家族和吸煙有關(guān). 既往研究顯示, APOBEC基因家族的突變印跡存在于乳腺癌、宮頸癌、肺癌等多個癌種, 其在慢性炎癥條件下的異常表達可能誤傷人類本身的基因組[25].

此外, 基于線性回歸的擬合方法雖然也擬合到了APOBEC基因家族, 但其相似性極低. 當(dāng)使用先驗知識干預(yù)時, 相似性有微弱的提升, 但并不具有統(tǒng)計效力.同時, 先驗知識的加入使得擬合結(jié)果中過濾掉了印跡8和28, 出現(xiàn)了印跡6, 而印跡6是一個“平坦”印跡, 所以先驗知識的加入也沒有使得擬合出現(xiàn)可接受的結(jié)果.最后, 該實驗驗證了在印跡分析問題中, NMF方法比LR方法適用性更強.

5 結(jié)束語

本文全面詳盡地探討了體細(xì)胞突變印跡分析的相關(guān)概念和模型, 并對算法進行了說明, 分類和比較； 闡述了基因組測序數(shù)據(jù)進行突變印跡提取的注意事項和缺陷. 此外, 使用真實的基因組數(shù)據(jù)全面而詳細(xì)地示范了如何在癌癥基因組的研究中應(yīng)用突變印跡分析.

本文介紹的分析框架使用了一些重要的限制和假設(shè)來描述基因組上的突變. 因此, 目前已經(jīng)提取到的突變印跡仍然是數(shù)學(xué)近似值, 其輪廓可能受到所用數(shù)學(xué)方法的影響. 從概念和使用的簡單性出發(fā), 研究中假設(shè)一個印跡與某個致癌因子引起的特定突變過程相關(guān)聯(lián),并以均一化的形式來表示. 然而, 不同的數(shù)學(xué)方法可以發(fā)現(xiàn)具有相似性和差異性的印跡, 并且這些印跡已通過多種方式得到證實. 隨著突變數(shù)量的增加和不同類型突變之間的數(shù)量級差異, 單一數(shù)學(xué)方法可能無法實現(xiàn)準(zhǔn)確的印跡分離, 因此結(jié)合充分的先驗知識, 進一步研究破譯和鑒定突變印跡的方法是避免產(chǎn)生在生物學(xué)上不可信或難以解釋結(jié)果的有效手段.

突變印跡的研究已經(jīng)發(fā)現(xiàn)了一些誘導(dǎo)癌癥發(fā)生發(fā)展的原因, 人類癌癥中自然發(fā)生的突變特征很可能有相當(dāng)一部分已經(jīng)被描述出來. 然而, 一些罕見的或者由治療導(dǎo)致的突變印跡可能還沒有被捕捉到, 需要進行徹底的探索. 目前許多最新發(fā)現(xiàn)的印跡背后的機制尚未明確, 還有待進一步地實驗和理解. 未來, 突變印跡的提取還應(yīng)該包含更多的突變類型, 無論是致病型還是繼發(fā)型, 同時也包括由遺傳導(dǎo)致的癌癥易感基因中的種系突變. 這些印跡背后的誘因?qū)τ诎┌Y預(yù)防和公共衛(wèi)生具有重要意義.