陳念華 袁細國

(西安電子科技大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，陜西 西安 710071)

0 引言

癌癥對人類的健康和生命威脅極大，從基因分子水平上研究癌癥的預(yù)防和治療策略是當代醫(yī)學(xué)急需解決的問題.近年來國際生物醫(yī)學(xué)界廣泛關(guān)注的一種新的基因組變異形式：拷貝數(shù)變異(copy number variation，CNV)，為此提供了新的線索和思路.CNV 是一種基因組結(jié)構(gòu)性變異，主要表現(xiàn)為長度從幾 Kb 至 Mb 的染色體片段的擴增或缺失[1, 2]，是促使人類個體間基因差異的重要因素之一，也是引發(fā)癌細胞產(chǎn)生和發(fā)展的重要現(xiàn)象.CNV擴增是指基因組區(qū)域的拷貝數(shù)從正常細胞二倍體到多倍體的變化，CNV缺失是基因組區(qū)域中拷貝數(shù)減少的變異.雖然CNV發(fā)生的頻率較低，但累積的堿基數(shù)量卻大大超過了單核苷酸多態(tài).在癌細胞中，CNV變異通常會引起相應(yīng)區(qū)域中包含的基因的劑量變化，這會影響基因的正常功能[3, 4].因此，在癌癥基因組中CNV的準確檢測對于癌細胞發(fā)展機理研究及癌癥診斷具有重要的現(xiàn)實意義[5, 6].

以多樣本數(shù)據(jù)為背景的CNV檢測與分析，其過程不僅涉及到癌癥樣本與正常樣本信號的比較，而且涉及到癌癥樣本本身之間的比較，那么依據(jù)CNV在樣本中出現(xiàn)的頻率，可將其分為復(fù)發(fā)性和個體性CNV模式.復(fù)發(fā)CNV指在多數(shù)樣本中共同發(fā)生的CNV區(qū)域，即CNV在多樣本中表現(xiàn)的頻率較高，目前相關(guān)檢測方法的研究非常豐富[7, 8].個體CNV指在少部分樣本中共同發(fā)生或個體特異性的CNV，即CNV在多樣本中表現(xiàn)的頻率較低[9].而目前為止，針對個體性CNV檢測的研究方法較少，但這種CNV模式同樣非常重要.通過研究個體CNV與癌癥的關(guān)系，不僅可以發(fā)現(xiàn)更多與癌癥發(fā)生發(fā)展密切相關(guān)的變異，還對在醫(yī)學(xué)上進行個體化的有針對性的藥物開發(fā)和治療有極大的幫助.

因此，本文提出一種名為IndivCNV(An individual copy number variation detection algorithm based on hierarchical matrix energy spectrum)的算法，與現(xiàn)有方法相比，該算法主要具有3個特點：(1) 可以從原始數(shù)據(jù)中實現(xiàn)個體性CNV模式的檢測；(2) 通過全變分將觀察到的信號進行平滑處理，利用潛變量模型將其重建為特征與權(quán)重的乘積，以應(yīng)對噪聲較高情況下CNV的檢測；(3) 對信號進行分層，根據(jù)分層矩陣能量譜在每層的占比，將能量高的復(fù)發(fā)CNV信號層剔除，以更準確鑒別個體性CNV.

1 相關(guān)工作

基于陣列的比較基因組雜交技術(shù)(array-based comparative genomic hybridization, aCGH)是一種高通量、高分辨率的方法，可以用于測量數(shù)千個DNA區(qū)域中拷貝數(shù)的變化.要從aCGH數(shù)據(jù)中檢測CNV，就必須定位信號數(shù)據(jù)中CNV區(qū)域與非CNV區(qū)域間的變化點，這些變化點會將染色體分成多個離散的片段，進一步便可以檢測出CNV.多樣本CNV的檢測涉及多個樣本，以期發(fā)現(xiàn)那些單樣本檢測無法發(fā)現(xiàn)的模式.目前有許多相關(guān)方法可以對aCGH數(shù)據(jù)進行多樣本CNV檢測，例如PLA(Piecewise-constant and low-rank approximation for identification of recurrent copy number variations)[10]、fastRPCA(A fused lasso latent feature model for analyzing multi-sample aCGH data)[11]、FLLat(A variational approach to stable principal component pursuit )[12]等.

PLA將多樣本CNV檢測問題轉(zhuǎn)化為矩陣分解問題，其中原始數(shù)據(jù)矩陣被分解為低秩分量、稀疏分量和噪聲分量.這三個成分分別對應(yīng)于復(fù)發(fā)CNV、個體CNV和隨機噪聲.通過主成分分析，也就是計算出輸入矩陣的奇異值分解，并使用前幾個奇異向量形成一個新的低秩矩陣，可以很容易地從低秩分量中識別出復(fù)發(fā)性CNV，從稀疏分量中識別出個體CNV.

類似地，fastRPCA采用線性疊加的模型，為穩(wěn)定主成分跟蹤(stable principal component pursuit, SPCP)引入了新的凸公式，將原始信號分解為低秩分量和稀疏分量.fastRPCA首先建立了一個凸變分框架，然后用準牛頓法對其進行加速，并使用此創(chuàng)新設(shè)計了通過變分框架的快速方法.用aCGH數(shù)據(jù)作為原始輸入，經(jīng)過以上處理，便可以從低秩分量中識別出復(fù)發(fā)性CNV，從稀疏分量中識別出個體CNV.

FLLat使用潛在特征模型對aCGH數(shù)據(jù)進行建模，其中每個樣本均通過固定數(shù)量的特征的加權(quán)組合來建模.這些特征代表了樣本組CNV的關(guān)鍵區(qū)域，并與權(quán)重相結(jié)合，描述了每個單獨樣本中的CNV區(qū)域.FLLat在特征的估計中使用了融合最小絕對值收斂和選擇算子，這在估計中既保證了數(shù)據(jù)的平滑度，也保證了數(shù)據(jù)的稀疏性.

以上這些方法雖然能較好的從多樣本數(shù)據(jù)中檢測出CNV，但是都不能對個體CNV進行針對性的檢測，因此本文提出了可以對個體CNV進行針對性檢測的算法IndivCNV.

2 方法

IndivCNV算法的基本框架如圖1所示，其輸入數(shù)據(jù)格式為大小為L×S的矩陣X，其中L代表探針數(shù)，S代表一組數(shù)據(jù)中包含的樣本個數(shù).該算法通過以下5個主要步驟實現(xiàn)對個體CNV的檢測：(1) 基于全變分正則化的信號層次化分解，(2) 應(yīng)用融合最小絕對值收斂和選擇算子，(3) 計算約束權(quán)重與特征數(shù)量J，(4) 模型參數(shù)估計，(5) 用分層矩陣能量譜識別個體CNV，下面將會針對每一個步驟的相關(guān)理論和實現(xiàn)過程進行詳細闡述.

2.1 基于全變分正則化的信號層次化分解

本文使用潛在特征模型來模擬多樣本數(shù)據(jù)，并且提出逐層分解信號的策略，通過將CNV的原始數(shù)據(jù)重建為不同特征模式的組合來發(fā)現(xiàn)原始數(shù)據(jù)中的CNV模式.將兩個秩為j的矩陣相乘的形式用j個秩為1的列向量與行向量相乘的加和來等價表示，以此來表示原始矩陣的分層分解，即

(1)

該模型說明了樣本組的CNV的重要特征是由J個特征共同總結(jié)的.具體來說，每個特征代表CNV的特定模式.然后，給定樣本的權(quán)重確定每個特征對該樣本的貢獻程度.換句話說，通過這些特征的權(quán)重可以知道不同特征的發(fā)生頻率，以此來推斷復(fù)發(fā) CNV和個體 CNV.

2.2 應(yīng)用融合最小絕對值收斂和選擇算子

CNV區(qū)域傾向于在整個染色體的連續(xù)區(qū)域中發(fā)生，區(qū)域具有相同的拷貝數(shù).對于未顯示CNV的染色體的其余部分，預(yù)期的信號強度應(yīng)為零.因此，如果我們將生物芯片數(shù)據(jù)視為沿著染色體的1維信號，則信號的大部分都為零，非零區(qū)域出現(xiàn)在平滑區(qū)域中.通過這種1維信號的稀缺性和平滑性的組合可以自然地想到融合最小絕對值收斂和選擇算子信號近似器(fused lasso signal approximator，F(xiàn)LSA[13]).FLSA可以解決優(yōu)化問題

(2)

其中u=(u1,…,up)T是估計所述有序結(jié)果的參數(shù)的向量.第一個懲罰項負責懲罰每個參數(shù)大小，這可以促使解決方案稀疏，第二個懲罰項負責懲罰相鄰參數(shù)之間的絕對差異，這可以促使解決方案平滑.有2個相應(yīng)的調(diào)諧參數(shù)，λ1和λ2，分別控制稀疏性和平滑性.

(3)

2.3 約束權(quán)重與特征數(shù)量 J

(4)

約束(4)對每行V的大小設(shè)置了限制，即對應(yīng)于給定特征的權(quán)重.在此認為這是限制權(quán)重大小的最合適方式.首先，它使估計的特征之間的直接比較更有意義；其次，它可以防止大部分權(quán)重僅分布在少數(shù)幾個特征上.

模型(1)中需要對特征J的數(shù)量做出選擇.從理論上講，J可以取{1，2,…，S}中的任何值，其中S是樣本數(shù).J的最好的選擇對于任何給定的數(shù)據(jù)集都是難以確定的，并可能取決于許多因素，例如，噪聲的水平，調(diào)諧參數(shù)λ1和λ2的值，以及S的值.因此，J的值通常留給用戶指定，默認設(shè)置為min {15，S/2}.本方法也提供選擇J的半自動過程，這是基于解釋的變化百分比(PVE).對于給定的J值，PVE被定義為

(5)

2.4 模型參數(shù)估計

(6)

(7)

2.4.2 選擇融合最小絕對值收斂和選擇算子調(diào)諧參數(shù)λ1和λ2.通常，給定模型的最佳調(diào)諧參數(shù)的選擇都是一個困難的任務(wù)，并且隨著調(diào)諧參數(shù)數(shù)量的增加會更加復(fù)雜.為了簡化對最佳調(diào)諧參數(shù)的搜索，本方法通過引入λ0和α∈(0,1)來重新定義參數(shù)λ1和λ2，使得λ1=αλ0，λ2=(1-α)λ0.在此可以認為λ0是整體調(diào)諧參數(shù)，它和α一起確定對稀疏度與平滑度的重視程度.通過固定α可能采取的值，可以有效地將對兩個參數(shù)λ1和λ2的搜索簡化為僅對一個參數(shù)λ0的搜索.

(8)

2.5 分層矩陣能量譜

(9)

(10)

其中T代表設(shè)定的占比閾值，I是大小為L×S的矩陣，代表個體 CNV.得到最終的個體CNV矩陣I以后，需要按照樣本將數(shù)據(jù)區(qū)分為S個大小為L×1的矩陣，每個矩陣代表每個樣本的結(jié)果.此時，需要再選定一個閾值H，若某探針處的絕對值大于H，則認為該處有個體CNV，反之則認為是正常.因為個體CNV在樣本間有很大的差異，所以需要按上述對每個樣本的結(jié)果數(shù)據(jù)都分別判斷.

3 實驗結(jié)果

3.1 模擬數(shù)據(jù)

3.1.1 模擬數(shù)據(jù)介紹. 為了評估IndivCNV算法對個體CNV的檢測性能，本節(jié)將采用模擬數(shù)據(jù)進行實驗，并與三種現(xiàn)有方法(PLA、FLLat、fastRPCA)進行比較.在文獻[15]里，詳細地定義了六種不同的復(fù)發(fā) CNV場景.在本文的研究里，將采用這六種場景來生成模擬數(shù)據(jù).在每一種場景下生成50組數(shù)據(jù)，每組數(shù)據(jù)是50×5000的矩陣，其中50代表50個樣本，5000代表每個樣本上的5000個探針.在生成每組數(shù)據(jù)時，無CNV區(qū)域的信號值設(shè)為0；復(fù)發(fā) CNV區(qū)域位于探針1876到3125之間，其模式參考圖2，將缺失變異區(qū)域的信號值設(shè)為-1，擴增變異區(qū)域設(shè)為1.每個樣本還需要在不與復(fù)發(fā) CNV區(qū)域重合的部分，隨機選取一個位置，添加一個長度為500探針的個體 CNV，個體 CNV的信號值從{-2，-1,1,2}中隨機選取，最后再向整個數(shù)據(jù)加入高斯噪聲.

6種不同場景生成模擬數(shù)據(jù)的過程展示在圖3，圖中黃色區(qū)域代表擴增，藍色區(qū)域代表缺失.其中第一行是根據(jù)文獻[15]中對不同場景的描述生成的只有復(fù)發(fā) CNV的數(shù)據(jù)，第二行是在復(fù)發(fā) CNV的基礎(chǔ)上隨機添加個體 CNV的數(shù)據(jù)，第三行是添加了噪聲水平為1的高斯噪聲的最終模擬數(shù)據(jù).每組數(shù)據(jù)的縱向代表樣本，橫向代表探針.從圖上可以看出，這六種場景可以分為兩類，場景1、3、5為一類，它們只有一個復(fù)發(fā)CNV區(qū)域；場景2、4、6為一類，它們含有多個復(fù)發(fā)CNV區(qū)域.本文的研究任務(wù)是從這些最終的模擬數(shù)據(jù)里準確恢復(fù)出個體 CNV.

3.1.2 檢測結(jié)果熱圖展示. 在圖4中展示了在6種場景下不同方法對個體 CNV的檢測結(jié)果.從圖中可以看出來IndivCNV檢測出了絕大部分的個體 CNV，并且能很好地把個體 CNV與復(fù)發(fā) CNV區(qū)分開來，沒有將復(fù)發(fā) CNV誤判為個體 CNV.FastRPCA可以分辨出一部分個體 CNV，但是沒有將噪聲很好地剔除，因此難以識別檢測出的個體 CNV的模式；而PLA則傾向于將一個完整的個體 CNV切割成多個小段，有明顯的缺失； FLLat的特點是它做檢測時不對復(fù)發(fā) CNV與個體 CNV進行區(qū)分，導(dǎo)致兩種類型的CNV都存在于結(jié)果數(shù)據(jù)中.由以上分析可知IndivCNV在檢測個體CNV時確實更加有優(yōu)勢，但是從圖中可以看出它還是存在一定的缺陷，因為它更趨向于檢測出發(fā)生個體CNV頻率較高的位置的變異，而對于發(fā)生頻率小的個體CNV，則很難檢測出.

3.1.3 檢測結(jié)果ROC曲線. 為了可量化地評估這些方法，本研究進一步通過ROC曲線評估各方法在六種場景下的個體 CNV識別性能.ROC(receiver operating characteristic curve)是一種顯示分類模型在所有分類閾值下的效果的圖表，其橫軸是假陽性率 (False Negative Rate，F(xiàn)PR)，縱軸是真陽性率 (True Negative Rate，TPR).FPR指的是所有非個體 CNV區(qū)域中被誤判為個體 CNV的比率，該值越小越好，TPR指的是在所有檢測出來為個體 CNV的區(qū)域里，確實是個體 CNV的比率，該值越大越好.ROC曲線的作用在于，在很多分類器分析中，得到的預(yù)測值通常不是0或1，而是一個0-1之間的概率值，此時就需要人為設(shè)定一個閾值，比如設(shè)定大于0.6則為1，反之則為零.但是不同的閾值所帶來的預(yù)測結(jié)果一定有差異，此時就可以用ROC曲線來刻畫不同閾值給分類器帶來的影響.通過上文對FPR和TPR含義的介紹可知，ROC曲線越靠近左邊沿和上邊沿，說明模型越好，因為此時TPR足夠大，F(xiàn)PR足夠小，說明分類的正確率很高.而ROC曲線上不同的點對應(yīng)著模型對不同閾值的預(yù)測水平，簡單來講，閾值越大，點越靠近左下，反之越靠近右上.

圖5展示了各方法在6種場景下的ROC曲線.這些ROC曲線是通過對各方法檢測出來的結(jié)果數(shù)據(jù)設(shè)定不同的閾值生成的.從圖上可以看出， IndivCNV檢測個體 CNV的性能優(yōu)于其他三種方法.例如在場景1的ROC曲線中，當FPR=0.1時，IndivCNV的TPR就已達到0.8，而FLLat的TPR只有0.45，PLA和fastRPCA的TPR僅有0.3；在場景2中，雖然當FPR值大于0.3時，F(xiàn)LLat和IndivCNV的曲線基本重合，但是IndivCNV在FPR=0.05時TPR就已經(jīng)達到了0.7，這說明IndivCNV在低FPR水平就可以實現(xiàn)較高水平的TPR；在場景3、5、6中，呈現(xiàn)出同樣的趨勢：當FPR較高時，F(xiàn)LLat與IndivCNV的曲線十分接近，但是始終都低于IndivCNV，只有在場景4中曲線的后半段FLLat超過了IndivCNV，盡管如此，其前半段依舊遠低于IndivCNV的ROC.

綜上所述，與fastRPCA和PLA相比，IndivCNV和FLLat算法對個體CNV的識別結(jié)果具有更高的TPR.然而，F(xiàn)LLat的性能與IndivCNV雖然較為接近，但僅表現(xiàn)在FPR較高的情況，當FPR較低時，其ROC曲線依舊遠低于IndivCNV.因此，在對個體CNV的檢測中， IndivCNV算法具有更明顯的優(yōu)勢.

3.2 真實數(shù)據(jù)

為了證明IndivCNV在真實數(shù)據(jù)上的可用性，本實驗引入異質(zhì)性乳腺癌CNA真實數(shù)據(jù)集對算法進行驗證.這個數(shù)據(jù)集中包含了112個乳腺癌樣本的SNP array數(shù)據(jù)，每個樣本都有23條染色體上的不同數(shù)據(jù)，每條染色體的探針各不相同，由Illumina 109 K SNP array平臺采集.在進行實驗時，首先將每個樣本不同染色體上的數(shù)據(jù)分割開來，然后將處理所得的CNV分段在基因組區(qū)域?qū)R，成為一個大小為112×pi的變異強度矩陣，其中112代表樣本數(shù)，pi代表在第i條染色體上的探針數(shù)，即分割完成后有22個變異強度矩陣(因為乳腺癌是常染色體上的疾病，所以僅對前22條常染色體進行實驗)，并分別對這22個信號矩陣進行實驗分析.在實驗過程中，使用IndivCNV對數(shù)據(jù)進行分析，閾值T設(shè)為0.1.為了消除每個樣本中的波譜偏差，需通過局部中值減去信號數(shù)據(jù)，中值計算的窗口大小是染色體長度的四分之一.

對于IndivCNV算法在乳腺癌數(shù)據(jù)中所發(fā)現(xiàn)的個體CNV區(qū)域，本研究通過乳腺癌相關(guān)文獻報道的CNV區(qū)域?qū)λ惴ńY(jié)果進行驗證.對于IndivCNV算法所發(fā)現(xiàn)的個體CNV區(qū)域，其中許多區(qū)域被現(xiàn)有文獻報道為乳腺癌CAN驅(qū)動區(qū)域.例如，IndivCNV算法成功識別出17號染色體上的ERBB2基因[16]，該基因曾被多項研究報道為乳腺癌CAN驅(qū)動變異.同時，IndivCNV在14號染色體發(fā)現(xiàn)AKT1基因[17]，而該基因則被報道與乳腺癌的發(fā)生發(fā)展密切相關(guān).表1匯總了IndivCNV所發(fā)現(xiàn)的個體CNV與現(xiàn)有文獻報道發(fā)現(xiàn)與乳腺癌有密切關(guān)系的基因重合的結(jié)果.上述結(jié)果表明，IndivCNV算法所發(fā)現(xiàn)的個體CNV區(qū)域與已報道CNV驅(qū)動變異區(qū)域具有較高的一致性.

表1 IndivCNV檢測出與現(xiàn)有文獻報道發(fā)現(xiàn)與乳腺癌有密切關(guān)系的基因重合的結(jié)果

表2 IndivCNV的復(fù)發(fā)CNV模式所匹配到的乳腺癌相關(guān)基因

正如第3節(jié)所說，IndivCNV在做個體CNV模式檢測的過程中，會將復(fù)發(fā)CNV的模式剔除.在此，本實驗在用該真實數(shù)據(jù)檢測時，將剔除的復(fù)發(fā)CNV數(shù)據(jù)也另行保存，并對復(fù)發(fā)CNV模式進行驅(qū)動基因匹配.表2中匯總了在IndivCNV的復(fù)發(fā)CNV模式中發(fā)現(xiàn)的乳腺癌驅(qū)動基因，表3匯總了IndivCNV在真實數(shù)據(jù)檢測出的個體CNV模式在復(fù)發(fā)CNV模式之外發(fā)現(xiàn)的驅(qū)動基因.由表2、3可以看出，個體CNV的檢測可以很大程度上彌補復(fù)發(fā)CNV對驅(qū)動基因發(fā)現(xiàn)的不足，例如，在表2復(fù)發(fā)CNV的檢測結(jié)果里，未發(fā)現(xiàn)1號染色體和13號染色體上有與乳腺癌相關(guān)的基因，而在個體CNV模式里則發(fā)現(xiàn)了1號染色體上的ARID1A基因，13號染色體上的BRCA2基因和RB1基因，這幾個基因都是乳腺癌相關(guān)基因，并被權(quán)威癌癥基因數(shù)據(jù)庫Cancer Gene Census所收錄[18-20].上述結(jié)果表明，IndivCNV算法的個體CNV發(fā)現(xiàn)結(jié)果可有效彌補現(xiàn)有方法發(fā)現(xiàn)結(jié)果的不足，同時也證明了個體CNV檢測對于癌癥研究的重要性.

表3 IndivCNV的個體CNV模式在其復(fù)發(fā)模式之外檢測到的乳腺癌相關(guān)基因

4 結(jié)論

CNV是導(dǎo)致癌癥發(fā)生發(fā)展的重要因素之一，由于現(xiàn)有研究更側(cè)重于對復(fù)發(fā)CNV的研究，對個體CNV的關(guān)注程度不夠，忽略了個體CNV的研究價值，因此本文通過分析個體CNV的模式，提出了一種新的適用于發(fā)現(xiàn)個體CNV的算法IndivCNV.IndivCNV首先需要使原始信號趨于平滑，因此采用了全變分正則化的方式達到此目的；接著將原始數(shù)據(jù)的每個樣本建模為固定數(shù)量的特征的加權(quán)和，這一步使用了潛變量模型和融合最小絕對值收斂和選擇算子懲罰；然后使用信號層次化分解，將不同模式的CNV用不同層的矩陣表示；最后利用分層矩陣能量譜，根據(jù)復(fù)發(fā)CNV模式能量占比大，個體CNV模式的能量占比小的原理，將復(fù)發(fā)CNV與個體CNV區(qū)分開來，最終達到檢測個體CNV的目的.

在本文的實驗中，首先將IndivCNV應(yīng)用到六種不同場景的模擬數(shù)據(jù)上，同時將fastRPCA、PLA、FLLat這三種算法也應(yīng)用到該模擬數(shù)據(jù)上，以ROC曲線為性能判斷標準，根據(jù)檢測結(jié)果選定不同閾值繪制ROC，以此進行性能對比，實驗結(jié)果表明，IndivCNV檢測個體CNV的性能顯著高于已有的三種方法的性能.然后又使用IndivCNV檢測異質(zhì)性乳腺癌CNA真實數(shù)據(jù)集中的個體CNV，檢測個體CNV結(jié)果中包含許多現(xiàn)有文獻已報道過與乳腺癌相關(guān)的基因，還發(fā)現(xiàn)了復(fù)發(fā)CNV模式?jīng)]有發(fā)現(xiàn)的與乳腺癌相關(guān)的基因，因此IndivCNV的性能在實際數(shù)據(jù)上也得到了驗證.綜上所述，IndivCNV在個體CNV方面的檢測性能確實有了大幅提升.