基于部分DNA甲基化預(yù)測腫瘤浸潤免疫細胞的比例

宋春暉，秦玉芳，陳明

1.上海海洋大學信息學院，上海 201306；2.農(nóng)業(yè)農(nóng)村部漁業(yè)信息重點實驗室，上海 201306

前言

DNA 甲基化狀態(tài)是基因組研究中一個突出的表觀遺傳學標記。隨著全基因組DNA甲基化數(shù)據(jù)的獲得，大量研究為DNA 甲基化在細胞過程和疾病如腎癌［1］、結(jié)直腸癌［2］中發(fā)揮的作用提供了證據(jù)，DNA 甲基化狀態(tài)為疾病的診斷和治療提供了很好的幫助。

腫瘤微環(huán)境中存在的浸潤性免疫細胞在癌癥進展、患者存活和對癥治療中起重要作用［3］。組織混合物細胞成分測定的實驗方法有流式細胞術(shù)和尖端單細胞技術(shù)，如Drop-seq［4］、10X 基因組學和sci-RNAseq［5］，但是實驗價格昂貴、勞動強度大，需要新鮮的組織，并且對細胞分離過程中的技術(shù)變化敏感。因此除了實驗方法，采用計算方法預(yù)測腫瘤組織中細胞成分的含量十分必要。

基于轉(zhuǎn)錄組的基因表達數(shù)據(jù)對復(fù)雜混合物（包括實體瘤）進行分析進而推斷細胞類型比例的方法目前有了一些研究［6-9］。Newman 等［6］采用支持向量回歸的方法CIBERSORT，利用微陣列數(shù)據(jù)對未知的復(fù)雜混合物進行細胞類型比例估計。基于細胞類型特異性mRNA 含量的重組，EPIC 考慮特征差異顯著的細胞類型，從大量的腫瘤轉(zhuǎn)錄組表達數(shù)據(jù)中估計腫瘤和免疫細胞類型的比例［7］。然而由于臨床上組織樣品的化學固定［10］，通過轉(zhuǎn)錄組分析方法測量的RNA 分子更容易降解，而DNA 甲基化是更穩(wěn)定的分子且具有高度的細胞類型特異性［11-12］，因此基于DNA 甲基化的方法成為細胞反卷積的一種更有效的替代方法。目前，對混合物中細胞含量的估計主要分為兩類，一類是基于參考的方法，Houseman等［13］將混合物樣本作為組成細胞類型的DNA甲基化與比例的加權(quán)組合［14］，利用約束投影/二次規(guī)劃推斷混合物樣本中的細胞類型比例；Teschendorff 等［15］利用來自NIH 表觀基因組的細胞類型特異性高靈敏位點信息構(gòu)建甲基化參考數(shù)據(jù)來推斷全血樣本中細胞類型的比例。另一類是基于無參考的方法，最初的兩種方法FaST-LMM-EWASher 和RefFreeEWAS 分別由James 等［16］和Houseman 等［17］提出，隨后Pavlo 等［18］在此基礎(chǔ)上，基于約束性非負矩陣分解并結(jié)合一種新的生物相關(guān)正則化函數(shù)開發(fā)了MeDeCom，并用于預(yù)測混合物中細胞類型的比例；Rahmani等［19］基于主成分分析開發(fā)了ReFACTor，該方法不需要細胞計數(shù)的先驗知識，提供了細胞類型組成的改進估計。雖然原則上基于無參考的方法可以應(yīng)用于任何組織，但這種方法的預(yù)測準確率較低；基于參考的方法能獲得較高的細胞比例預(yù)測準確率，然而在實際臨床上往往只知道一部分細胞類型的甲基化，即部分參考的情況，因此可以利用容易獲得的腫瘤混合物中部分細胞類型的表觀基因組信息去推斷出所有組成細胞類型的比例?；诓糠謪⒖嫉姆椒ㄒ殉晒?yīng)用于基于轉(zhuǎn)錄組數(shù)據(jù)的腫瘤浸潤免疫細胞的分解［20］。

本研究基于腫瘤組織中已知細胞類型的甲基化，利用非負矩陣分解的框架去估計所有細胞類型的比例，簡記為MethyPR。對模擬數(shù)據(jù)的評估表明本研究的方法較現(xiàn)有方法在識別細胞類型比例上有了明顯的提高；其次，在體外制備的混合物上驗證了本研究的方法能很好地還原出所有細胞類型的比例；最后，將MethyPR 應(yīng)用于癌癥基因組圖譜（TCGA）甲基化數(shù)據(jù)，能很好地識別出癌癥特異性腫瘤浸潤性免疫細胞，為靶向免疫治療的設(shè)計提供依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)集獲取

1.1.1 模擬混合物數(shù)據(jù)本研究從基因表達綜合數(shù)據(jù)庫（GEO）下載了6 種純化免疫細胞［CD4+T 細胞、CD8+T細胞、自然殺傷細胞、B細胞、單核細胞（Mon）和粒細胞（Gra）（GSE35069）］和一種乳腺癌細胞（MCF-7）（GSE44837）的甲基化譜，并按照一定比例產(chǎn)生免疫細胞和腫瘤細胞的混合物。

1.1.2 實驗獲得的真實樣本數(shù)據(jù)真實樣本的甲基化數(shù)據(jù)以及組成樣本的細胞類型甲基化數(shù)據(jù)來源于Onuchic等［21］的研究，具體來說，該數(shù)據(jù)集包含兩個樣本集和構(gòu)成樣本的細胞類型甲基化數(shù)據(jù)。第一個樣本集由6個樣本組成，分別由3個成對細胞系的組合組成（MCF-7/HMEC、MCF-7/CD8+T 細胞和MCF-7/CAF），分別按75%：25%和95%:5%的比例混合組成。第二個樣本集是由亞硫酸氫鹽測序生成的29個乳腺腫瘤樣本組成，每個樣本由不同的乳腺癌細胞系、一種正常乳腺細胞系（HMEC）、一種成纖維細胞系（CAF）和一種免疫細胞（CD8+T）組成，同時使用H&E染色，估計每個樣本的癌性、正常、成纖維和免疫細胞的比例。該數(shù)據(jù)集還包含了構(gòu)成樣本的細胞系的甲基化數(shù)據(jù)，分別是6種不同的乳腺癌細胞系（MCF-7、T47D、MDA-MB-231、MDA-MB-361、HCC1954、HCC1569）、HMEC、CAF和CD8+T細胞。

1.1.3 TCGA癌癥樣本數(shù)據(jù)本研究使用GDC客戶端工具從GDC 數(shù)據(jù)門戶（https://gdc.nci.nih.gov）下載胸腺瘤、結(jié)直腸癌、急性髓系腫瘤和彌漫性大B 細胞淋巴瘤樣本的Infinium HumanMethylation 450 芯片的三級甲基化數(shù)據(jù)。另外，從Arneson 等［22］的研究得到免疫細胞的甲基化數(shù)據(jù)，該數(shù)據(jù)集包含11 種細胞類型（單核細胞、樹突狀細胞、巨噬細胞、中性粒細胞、嗜酸性粒細胞、調(diào)控性T 細胞、幼稚T 細胞、記憶T 細胞、CD8+T 細胞、自然殺傷細胞和B 細胞）的甲基化數(shù)據(jù)（GSE35069、GSE59250、GSE71837），并且按照文獻［22］的方法，采用單核細胞作為橋接細胞類型糾正批次效應(yīng)，最終得到419個甲基化位點處的腫瘤浸潤性免疫細胞甲基化數(shù)據(jù)。

1.2 數(shù)據(jù)預(yù)處理

對于下載得到的高維甲基化數(shù)據(jù)，首先排除甲基化位點缺失超過10%的樣本以及樣本缺失超過10%的甲基化位點。然后，用R 包“ChAMP”包進行一系列處理：使用ChAMP.impute 方法對缺失值進行填充，使用ChAMP.filter 方法過濾常染色體上與性相關(guān)的和位于X、Y 染色體上的位點，以避免任何與性別相關(guān)的信號，同時過濾所有與單核苷酸多態(tài)性（Single Nucleotide Polymorphism, SNP）重疊的和報道為交叉反應(yīng)的位點［23］。例如，對于模擬實驗中6種免疫細胞的高維甲基化數(shù)據(jù)，排除6個在每種細胞中缺失值超過10%的位點，接下來過濾10 028個與性相關(guān)和位于X、Y 染色體上的位點，使用一般的450K SNP列表過濾了59 901個SNP重疊的位點，并移除了11個交叉反應(yīng)的位點，最終得到6種免疫細胞對應(yīng)的412 481個位點處的甲基化水平。

為了獲得對每種細胞類型差異顯著的位點，對上述數(shù)據(jù)進行進一步處理，使用ChAMP.DMP方法提取細胞類型中差異顯著的位點［24］。接下來，遵循文獻［25］對上述得到的數(shù)據(jù)集中的每種細胞類型的甲基化水平與其他細胞類型甲基化水平進行t 檢驗比較，得到差異顯著（P＜0.000 1）的甲基化位點，按照P值進行降序排列，選擇出1 000個甲基化位點，并選擇最終的甲基化位點為每次t檢驗中選擇出的位點的交集。例如：對于模擬實驗中的6 種免疫細胞，通過“ChAMP”步驟的處理，得到25 820 個對于每種細胞類型差異顯著的位點，經(jīng)過t檢驗的處理后，最終得到包含412個顯著差異的位點處的甲基化水平。

1.3 模型

1.3.1 模型建立MethyPR 的模型如圖1 所示。M∈Rm×n是來自k個細胞類型組成的一個行為m個位點，列對應(yīng)n個樣本的DNA 甲基化水平矩陣。其中Mij∈M，表示第i個位點在第j個樣本的甲基化探針占探針總數(shù)的比例，所以0 ≤Mij≤1。本研究用W∈Rm×k表示每個位點的細胞類型甲基化水平，即特征矩陣，其中Wih∈W；此外，用H∈Rk×n表示k個細胞類型在樣本中的比例矩陣，Hhj∈H。根據(jù)文獻［26］，常用的DNA甲基化混合物模型是

其中，誤差項εij服從正態(tài)分布。式（3）和（4）中的約束要求細胞比例是正的，并在每個樣本中的總和為1；式（5）中的約束要求細胞類型的甲基化水平在［0,1］范圍內(nèi)。

圖1 MethyPR模型的框架Figure 1 Framework of MethyPR model

在大多數(shù)情況下，矩陣W的數(shù)據(jù)并不是完全已知，即研究中討論的細胞類型部分已知的情況。此時模型可以修改為：

矩陣W分為W1和W2，即式（6）中的Wih(1 ≤h≤k1)和Wih(k1+ 1 ≤h≤k)，W=(W1,W2)，樣本細胞比例矩陣H分為H1和H2，即式（6）中的Hhj(1 ≤h≤k1) 和Hhj(k1+ 1 ≤h≤k)，，其中W1和H1是已知的細胞類型所對應(yīng)的甲基化水平和所占的比例，W2和H2是未知細胞類型所對應(yīng)的甲基化水平和所占的比例。

給定腫瘤樣本的甲基化矩陣M和已知細胞甲基化矩陣W1，目標是得到所有細胞類型對應(yīng)的比例矩陣H。這個問題可以轉(zhuǎn)化為尋找合適的H1、H2和W2，使全局誤差平方和最小：

1.3.2 模型求解采用迭代的非負矩陣分解方法來求解模型，具體流程如下：

（1）隨機給定W2矩陣的初始值；

（4）不斷重復(fù)（2）和（3）的步驟，直至收斂或者達到預(yù)定的最大迭代次數(shù)。

其中，t為迭代次數(shù)。在第（1）步中，使用RPMM 算法［27］來初始化未知細胞類型的甲基化水平。在運用MethyPR時要注意，該方法適用于下載的定量基因組的原始數(shù)據(jù)，如原始reads 數(shù)或CpG 位點甲基化水平，不適合使用logit 變換的數(shù)據(jù)，因為變換會破壞數(shù)據(jù)中存在的線性關(guān)系。

1.3.3 細胞類型數(shù)量確定和模型質(zhì)量評估采用赤池信息準則（Akaike Information Criterion,AIC）度量來確定腫瘤樣本中細胞類型的數(shù)量。AIC 是衡量統(tǒng)計模型擬合優(yōu)良性的一種方法，它不僅考慮了模型的擬合優(yōu)度，也考慮了當模型成分數(shù)量增加時可能發(fā)生的過度擬合。由于本實驗中的樣本量小，因此采用適用于樣本量小的改進版本AICc。AICc 的公式為：

其中，n表示樣本量；k表示模型參數(shù)的數(shù)量（細胞類型的數(shù)量）；SSR 為殘差平方和。不同的k值在模型中對應(yīng)不同的AICc 值，在所有結(jié)果中，AICc 最小時所對應(yīng)的k就是最優(yōu)的細胞類型數(shù)。

采用兩個評價指標來評估MethyPR 的性能，分別為均方根誤差（Root Mean Squared Error, RMSE）和Pearson 相關(guān)系數(shù)。RMSE 是指對預(yù)測值與真實值差平方的平均值求平方根，這是回歸問題常用的性能指標；Pearson 相關(guān)系數(shù)用于衡量真實值和預(yù)測值之間的相關(guān)程度，其值在［-1,1］之間，絕對值越接近1，相關(guān)性越強。本研究用模型預(yù)測的細胞類型比例和真實細胞比例的Pearson相關(guān)系數(shù)來衡量評估方法的精度，另一方面，使用真實細胞類型的比例和預(yù)測的細胞比例之間的RMSE作為評價指標。

2 結(jié)果與分析

2.1 模擬混合物驗證

首先，在模擬混合物上對本研究提出的方法進行基準測試，模擬混合物包含6 種純化的免疫細胞（CD4+T 細胞、CD8+T 細胞、巨噬細胞、B 細胞、單核細胞和粒細胞）和一種乳腺癌癌癥細胞（MCF-7）。下載這7種細胞類型的甲基化水平之后，將腫瘤細胞以90%的含量添加到混合物中，剩余部分由6 種免疫細胞類型按隨機比例構(gòu)成。

對于上述基準數(shù)據(jù)集，將現(xiàn)有的3 種方法（QP［13］、CIBERSORT［6］、EDEC［21］）與本研究提出的方法進行對比。QP 是基于參考的方法，直接使用組成細胞類型的參考甲基化譜，將其與混合物樣本甲基化譜作為輸入，利用二次規(guī)劃來推斷混合物樣本中的細胞類型比例。CIBERSORT 內(nèi)核為支持向量機的一個實例（ν - SVR），將其與細胞類型差異顯著的特征矩陣相結(jié)合，從而預(yù)測細胞類型的比例，其中參數(shù)ν給出了訓練誤差的上界和支持向量的界。EDEC假設(shè)混合物樣本的對應(yīng)位點甲基化譜為細胞類型特異性甲基化譜和細胞類型比例的線性組合，隨機初始化細胞類型甲基化譜，通過二次規(guī)劃求解約束最小二乘問題，得到估計的細胞類型比例。

已知有6 種免疫細胞，本研究假設(shè)模擬混合物的細胞類型數(shù)量為7～15，將細胞數(shù)量、模擬混合物數(shù)量以及預(yù)測結(jié)果與模擬結(jié)果的誤差作為AICc指標的輸入，進而推斷混合物中細胞類型的總數(shù)量。由圖2發(fā)現(xiàn)，當細胞類型數(shù)量為7 時AICc 值最小，最小AICc值識別的細胞類型數(shù)量與我們構(gòu)建模擬混合物時的數(shù)量一致。接下來，對于腫瘤含量為90%的混合物，將混合物樣本和6 種純化免疫細胞經(jīng)過預(yù)處理后對應(yīng)的位點處的甲基化水平作為方法的輸入，而乳腺癌細胞的甲基化水平未知，用4種方法分別估計了模擬混合物中6種免疫細胞類型的比例。

如圖3 所示，分別計算了4 種方法預(yù)測的免疫細胞實際比例與預(yù)測比例的RMSE 和Pearson 相關(guān)系數(shù)。從圖3a 可以看出，本研究方法在組成混合物樣本的所有免疫細胞中取得最小的RMSE 分別為0.005 4、0.003 3、0.006 4、0.006 0、0.008 4、0.005 6）；圖3b顯示了4 種方法在所有免疫細胞類型上的RMSE 平均值，可以觀察到本研究方法明顯優(yōu)于其他方法。圖3c 和圖3d 為4 種方法在Pearson 相關(guān)系數(shù)指標下的比較，MethyPR 的性能顯著高于其他方法，平均Pearson相關(guān)系數(shù)指標達到0.97。

圖2 識別模擬混合物中細胞類型數(shù)量的AICc線圖Figure 2 AICc line graph for identifying the number of cell types in the simulated mixture

為了測試本研究方法對不同噪聲水平和生成腫瘤混合物時不同細胞類型含量下的穩(wěn)健性，進一步進行了以下實驗。固定混合物中腫瘤含量為90%，在模擬混合物中加入不同的噪聲水平，從0.1 逐漸增加到0.5，比較MethyPR 和其他3 種方法在不同噪聲數(shù)據(jù)上的預(yù)測性能，其中每次模擬重復(fù)20次。另外，固定噪聲水平為0.1，將腫瘤細胞的含量從90%降低到50%（每次減少10%）添加到混合物中，混合物樣本中的剩余部分由6種免疫細胞類型按隨機比例構(gòu)成，每次模擬同樣重復(fù)20 次。在這兩個參數(shù)的不同設(shè)置下，估計6種免疫細胞類型的比例。

圖3 MethyPR與其他算法在細胞水平上的比較Figure 3 Comparison of MethyPR and other algorithms at the cellular level

在腫瘤含量分別為50%、60%、70%、80%、90%的情況下，計算AICc 值。由表1 可知，混合物的細胞類型數(shù)量均在7時取得最小，這與構(gòu)建混合物的數(shù)量相一致，說明采用AICc 值預(yù)測腫瘤細胞類型數(shù)量得到了很好的結(jié)果。由表2 可知，在不同的噪聲水平下，也得到了類似的結(jié)果。

比較4 種方法在不同噪聲水平和不同腫瘤細胞含量下的預(yù)測性能（圖4）。在不同的噪聲水平下，得到細胞比例的RMSE 和Pearson 相關(guān)系數(shù)如圖4a 和圖4b 所示。MethyPR 在不同噪聲水平下均優(yōu)于其他方法，對于低噪聲水平的性能是最佳的，但是在高噪聲下，QP 的性能逐漸趕上MethyPR。從圖4c 可以看到，MethyPR在不同腫瘤細胞含量下取得穩(wěn)定的較小RMSE，而其他方法隨著未知腫瘤含量的增加，RMSE越來越大；圖4d中，隨著腫瘤細胞含量的增加，4種方法的預(yù)測性能逐漸下降，但在所有情況下，MethyPR具有最高的Pearson相關(guān)系數(shù)。

表1 生成混合物中不同腫瘤細胞含量下的AICc值Table1 AICc values of the generated mixture containing different tumor cell contents

表2 生成混合物在不同噪音水平下的AICc值Table 2 AICc values of the generated mixture at different noise levels

圖4 不同方法在不同噪聲水平和不同腫瘤細胞含量下的性能評估Figure 4 Performance evaluations of different methods at different noise levels and different tumor cell contents

2.2 體外制備的細胞系混合物驗證

在Onuchic 等［21］體外制備的細胞混合物上驗證本研究提出的方法。從GEO 下載MCF-7、HMEC、CAF和CD8+T細胞的甲基化譜，由3對組合（MCF-7/HMEC、MCF-7/CD8+T 細胞和MCF-7/CAF）各按兩種比例（75%:25%和95%:5%）生成6 個樣本，每個樣本通過靶向亞硫酸氫鹽測序進行分析。采用數(shù)據(jù)預(yù)處理方法得到149 個在不同乳腺癌細胞類型具有顯著差異的位點并用于本研究方法。在癌癥細胞MCF-7 未知的情況下（圖5a），MethyPR 估計的免疫細胞比例和真實值之間有很強的一致性（R=0.996），同時組成混合物的MCF-7 細胞比例與真實值也有很強的一致性（R=0.993）。

臨床病理學家根據(jù)H&E 染色對29 個乳腺腫瘤樣本進行細胞類型組成評估，估計癌性、正常、基質(zhì)和免疫細胞的比例。本研究假定癌性細胞未知，將其余3 種細胞類型差異顯著的特征矩陣和29 個樣本的混合甲基化矩陣作為輸入，用MethyPR 估計29 個樣本中細胞的比例，觀察到估計的CD8+T 免疫細胞與真實免疫細胞的比例具有較高的一致性（R=0.81），而癌癥細胞的相關(guān)性與之相比較低（R=0.71），見圖5b。從圖5 可以看出，MethyPR 可以準確地預(yù)測免疫細胞和剩余其他細胞類型的比例。

圖5 MethyPR在體外制備的混合物上的驗證Figure 5 Validation of MethyPR on the mixture prepared in vitro

2.3 腫瘤免疫浸潤細胞成分估計

從TCGA的4種腫瘤樣本中分別選取了100個樣本，然后將MethyPR 應(yīng)用到這些數(shù)據(jù)上，估計樣本中浸潤免疫細胞的比例。如表3所示，不同類型的腫瘤樣本表現(xiàn)出不同的免疫細胞浸潤模式。幼稚T 細胞在胸腺瘤樣本中占有最高的比例（0.031），這與之前的實驗研究一致［28］。在結(jié)直腸癌樣本中，以T細胞為主［29］，而急性髓系腫瘤樣本具有高比例的單核細胞，彌漫性大B細胞淋巴瘤樣本中的B細胞比例最高，這分別與文獻［21］和文獻［30］結(jié)果相一致。

表3 腫瘤浸潤免疫細胞在TCGA樣本中的比例Table 3 Proportions of tumor-infiltrating immune cells in TCGA samples

3 結(jié)論

本研究提出了一種使用甲基化數(shù)據(jù)對腫瘤微環(huán)境進行穩(wěn)健反卷積的方法MethyPR，該方法基于非負矩陣分解方法，利用容易獲得的免疫細胞類型的表觀基因組信息，從DNA 甲基化估計細胞組成比例。作為解決基于參考與無參考方法局限性的新方法，MethyPR 能夠基于部分參考數(shù)據(jù)推斷混合物中的細胞類型比例。在模擬混合物、實驗混合物和真實混合物上，本研究方法都表現(xiàn)出良好的性能，其預(yù)測精度高，可以幫助減少估計腫瘤組成的時間和金錢成本，快速獲得腫瘤混合物的免疫細胞比例，為表觀基因組研究提供了新的思路。