李慧月， 王沐榛， 喻 葉， 田雪梅

(華南師范大學(xué)生命科學(xué)學(xué)院， 廣州 510631)

胰腺癌是人類最具侵襲性和致命性的惡性腫瘤之一，胰腺癌的早期診斷較為困難，大多數(shù)患者發(fā)現(xiàn)時已經(jīng)是晚期. 手術(shù)切除、化療在胰腺癌治療中的效果十分有限[1]. 免疫治療是近年來腫瘤治療研究的熱點，在非小細胞肺癌和黑色素瘤等實體瘤中展現(xiàn)了不錯的臨床療效，但由于胰腺癌獨特的腫瘤微環(huán)境和較低的腫瘤免疫原性，限制了免疫治療對胰腺癌的療效[2-4]. 腫瘤浸潤性免疫細胞是腫瘤微環(huán)境的主要細胞成分，其組成和功能可隨著宿主免疫狀態(tài)的改變而改變，并與腫瘤細胞之間存在復(fù)雜的相互作用，影響腫瘤的進展和預(yù)后[5-7]. 全面揭示胰腺癌腫瘤浸潤性免疫細胞的特征，了解胰腺癌的免疫表型對胰腺癌患者的早期診斷和治療具有重要意義.

加權(quán)聯(lián)合表達網(wǎng)絡(luò)分析(Weighted Gene Co-expression Network Analysis, WGCNA)可用于探索疾病臨床特征與基因簇之間的相關(guān)性，以確定癌癥的相關(guān)模塊和中心基因[8]. 該方法已被廣泛用于在轉(zhuǎn)錄組水平尋找生物標志物[9]. 本研究使用反卷積CIBERSORT法計算樣本的腫瘤浸潤性免疫細胞的相對比例，結(jié)果顯示B細胞所在的模塊相關(guān)性最高，接著通過WGCNA識別出與B細胞相關(guān)的重要模塊和樞紐基因，隨后利用公共數(shù)據(jù)庫進行一系列驗證，最終分析篩選出胰腺癌中與B細胞免疫浸潤相關(guān)的分子標志物.

1 材料和方法

1.1 數(shù)據(jù)來源和預(yù)處理

人胰腺癌樣本的mRNA測序數(shù)據(jù)從TCGA數(shù)據(jù)庫( https:∥portal．gdc．cancer．gov/)的“PAAD”隊列下載得到，共獲得182個樣本數(shù)據(jù)，包括178例胰腺癌組織樣本和4例正常胰腺樣本. 同時下載臨床資料數(shù)據(jù)(Clinical). 在R語言環(huán)境下運行“l(fā)imma”包對基因表達數(shù)據(jù)進行歸一化處理后用于后續(xù)分析.

1.2 腫瘤浸潤性免疫細胞的評估

CIBERSORT反卷積法是基于已知參考集LM22，通過分析 mRNA 表達數(shù)據(jù)來定量評估每個樣本中 22 種免疫細胞亞型的相對比例的分析方法[10]，目前已廣泛運用于腫瘤相關(guān)免疫細胞的研究中. 在本研究中，使用R包“CIBERSORT”，基于參考基因表達特征集“LM22”，將TCGA中胰腺癌的mRNA表達矩陣導(dǎo)入算法進行反卷積分析和比對，得到胰腺癌樣本中22種免疫細胞的相對比例.

1.3 構(gòu)建加權(quán)基因共表達網(wǎng)絡(luò)和模塊分析

加權(quán)基因共表達網(wǎng)絡(luò)分析WGCNA是用來分析不同樣本之間基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法，可以用來鑒定高度協(xié)同變化的基因表達矩陣，并根據(jù)表達矩陣的內(nèi)連性以及基因與表型之間的關(guān)聯(lián)度鑒定候選核心基因. 本研究使用R包“WGCNA”鑒定相關(guān)基因模塊，篩選可能與胰腺癌免疫相關(guān)的基因集. 通過分析不同加權(quán)系數(shù)下模塊的尺度獨立性和平均連通性來確定模塊分析的軟閾值β. 軟閾值確定后，構(gòu)建無尺度拓撲分布網(wǎng)絡(luò)，根據(jù)基因間的Pearson相關(guān)系數(shù)將相關(guān)矩陣轉(zhuǎn)換為鄰接矩陣，進一步轉(zhuǎn)換為拓撲重疊矩陣(Topological Overlap Matrix，TOM)，得到基因間的相異度(1-TOM). 利用層次聚類函數(shù)，將表達譜相似的基因歸類于同一個基因模塊.

分析各模塊的組成部分. 將模塊與臨床特征進行關(guān)聯(lián)分析，本研究將22種免疫細胞作為臨床特征納入共表達網(wǎng)絡(luò)，采用Pearson檢驗來計算各模塊與免疫細胞浸潤程度之間的相關(guān)性. 當P<0.05時，選取相關(guān)系數(shù)最高的模塊定義為中心模塊.

為了確定基因在中心模塊中的功能，使用網(wǎng)頁工具“Metascape”(http:∥metascape.org)用于通路富集分析[11]，設(shè)置P<0.05為截止點.

1.4 樞紐基因的識別

通過WGCNA分析模塊內(nèi)連通性篩選前30個候選基因，顯著性模塊定義為P<0.05. 選取中心模塊中的所有基因，使用檢索工具String(https:∥string-db.org/)檢索相互作用的基因[12]，并構(gòu)建PPI網(wǎng)絡(luò)，設(shè)置置信度評分大于0.4為截止標準，以TSV格式導(dǎo)出. 然后將PPI網(wǎng)絡(luò)導(dǎo)入到Cytoscape工具中進行可視化操作[13]，利用CytoHubba插件的Degree評分來識別前30個中心節(jié)點. 使用網(wǎng)頁工具Venny 2.1.0(https:∥bioinfogp.cnb.csic.es/tools/venny/)對WGCNA分析得到的30個候選基因和PPI分析得到的30個中心節(jié)點繪制韋恩圖，取交集內(nèi)所包含的基因作為本研究的樞紐基因.

TIMER(https:∥cistrome.shinyapps.io/timer/)是一個系統(tǒng)分析不同癌癥類型的免疫浸潤的綜合資源[14]，它通過反卷積的統(tǒng)計方法，從基因表達譜推斷腫瘤浸潤性免疫細胞的相對比例[15]. 使用TIMER來計算目標免疫細胞的浸潤水平與樞紐基因表達之間的Spearman相關(guān)性.

1.5 差異表達基因和生存分析

原始數(shù)據(jù)在R語言環(huán)境(version3.5.3, https:∥www.r-project.org/)下使用“affy”包進行預(yù)處理和標準化[16]. 使用“l(fā)imma”包篩選胰腺癌患者和對照組表達數(shù)據(jù)中的差異表達基因(DEGs)，設(shè)定閾值：|log2FC|>2，P<0.05；采用Kaplan-Meier分析和logrank檢驗評估樞紐基因表達對胰腺癌患者總生存率的影響[17].

為了驗證樞紐基因的表達方式，使用了來自癌癥微陣列數(shù)據(jù)庫Oncomine(https:∥www.oncomine.org)的4個獨立的微陣列數(shù)據(jù)集進行meta分析[18]. Oncomine將大量公布的癌癥微陣列轉(zhuǎn)錄組數(shù)據(jù)與先進的分析工具結(jié)合在一起，能夠分析人類癌癥中原癌和抑癌基因的通路、調(diào)控網(wǎng)絡(luò)和功能網(wǎng)絡(luò).

1.6 統(tǒng)計分析

根據(jù)數(shù)據(jù)分布特征，采用非參數(shù)檢驗或t檢驗分析實驗組和對照組之間差異的統(tǒng)計學(xué)意義. 所有R語言分析均用R3.5.3軟件進行，P<0.05被認為有統(tǒng)計學(xué)意義.

2 結(jié)果

2.1 mRNA表達數(shù)據(jù)

本研究獲得了來自癌癥基因組圖譜(The cancer genome atlas,TCGA)數(shù)據(jù)庫的182個樣本RNA seq數(shù)據(jù)，包括178例胰腺癌組織樣本和4例正常胰腺樣本. 將mRNA表達數(shù)據(jù)合并成一個行為基因名、列為樣本名的標準化矩陣，得到一個包含19 669個基因、182個樣本的表達矩陣. 為剔除樣本中與其他基因聯(lián)系不夠密切的基因，簡化運算，基于R環(huán)境下使用“apply”函數(shù)選取了變化量較大的(方差處于前四分之一)的4 916個基因進行進一步分析.

2.2 腫瘤浸潤性免疫細胞的評估

CIBERSORT計算胰腺癌中的22種腫瘤浸潤性免疫細胞成分，選擇7種免疫細胞：B cells memory ，NK cells activated，Dendritic cells activated，Macrophages M1，Macrophages M2，Neutrophils，Monocytes作為WGCNA的臨床特征數(shù)據(jù)，保存為文本格式用于后續(xù)分析.

2.3 加權(quán)基因共同表達網(wǎng)絡(luò)的構(gòu)建和模塊分析

利用R包“WGCNA”構(gòu)建4 916個基因表達值的共表達網(wǎng)絡(luò)，計算平均連接值和Pearson相關(guān)值對樣本進行聚類. 當以0.9作為相關(guān)系數(shù)時，軟閾值β=12. 通過WGCNA分析，構(gòu)建了13個共表達模塊，每個模塊至少包含30個基因. 以模塊特征基因和免疫細胞的Pearson相關(guān)系數(shù)表示二者的相關(guān)性，相關(guān)性排名前3的模塊分別為：記憶B細胞所在的cyan模塊、樹突狀細胞所在的royalblue模塊以及NK細胞所在的greenyellow模塊，相關(guān)系數(shù)分別為0.52、0.45和0.37(圖1)，選擇其中相關(guān)系數(shù)最高的記憶B細胞所在的cyan模塊內(nèi)的基因進行后續(xù)分析. 使用Metascape工具對該模塊中的基因進行通路分析和過程富集分析(圖2)，“B細胞增殖”是第3個高度富集的條目，進一步說明了cyan模塊和B細胞的高度相關(guān)性.

圖2 cyan模塊內(nèi)基因富集的生物學(xué)過程條目

2.4 樞紐基因的識別

通過以下步驟篩選胰腺癌中與B細胞相關(guān)的樞紐基因. 首先，基于WGCNA分析在cyan模塊中進行基因連通性排序，篩選出前30個候選基因；其次，在PPI網(wǎng)絡(luò)中，運用Cytoscape的插件CytoHubba對網(wǎng)絡(luò)中的基因進行連接度Degree值排序，得到Degree評分前30個節(jié)點作為中心節(jié)點，置信度評分大于0.4 (圖3)；最后，進行韋恩分析，得到候選基因與中心節(jié)點的交集，共篩選出9個樞紐基因：CD79B、MYC、BANK1、TIMELESS、CD19、ATF3、ITGAL、IKZF3、RRAGB.

為了研究這些樞紐基因與B細胞的關(guān)系，分析了這些基因在TIMER數(shù)據(jù)庫中的表達數(shù)據(jù)，結(jié)果顯示：9個基因的表達值與B細胞浸潤水平呈正相關(guān)(除ATF3外，所有基因的P值均小于0.05)(圖4).

圖4 樞紐基因與B細胞浸潤程度的相關(guān)性

2.5 差異表達基因和生存分析

利用R語言的“l(fā)imma”包對TCGA中的胰腺癌組織和正常組織的基因表達數(shù)據(jù)進行處理，得到了9個基因的表達水平，發(fā)現(xiàn)ITGAL在腫瘤組織中的表達水平明顯高于正常對照組(|log2FC|>2，P<0.05)(圖5). 采用Kaplan-Meier分析方法對9個樞紐基因進行分析，結(jié)果顯示MYC、BANK1、ITGAL、RRAGB具有統(tǒng)計學(xué)意義(P<0.05)(圖6). 對于MYC和BANK1，高表達患者生存預(yù)后較差，而ITGAL和RRAGB高表達的患者預(yù)后較好.結(jié)合差異表達基因和生存分析，選擇ITGAL作為潛在的生物標志物進行進一步分析.

圖5 差異基因表達的火山圖

圖6 樞紐基因的Kaplan-Meier曲線

為了驗證ITGAL在腫瘤組織和正常組織中的差異表達，對Oncomine數(shù)據(jù)庫中3個數(shù)據(jù)集的4個分析進行了meta分析. 結(jié)果顯示：ITGAL在腫瘤組織中顯著高表達，這與TCGA數(shù)據(jù)集的分析一致(Median Rank Values=1 809,P<0.05)(圖7).

圖7 基于Oncomine數(shù)據(jù)集的基因表達的meta分析

3 討論

胰腺癌是全球最致命的惡性腫瘤之一，常規(guī)的治療手段如手術(shù)、化療和放療效果不佳. 近年來，免疫治療在一些惡性腫瘤中取得了很大的成功，不斷有新的靶點被發(fā)現(xiàn)，免疫治療也是人們不斷探索的治療胰腺癌的方法[19]. 盡管在癌癥研究領(lǐng)域使用WCGNA已經(jīng)發(fā)現(xiàn)了的一些關(guān)鍵生物標志物[20-21]，但大多數(shù)研究局限使用差異表達基因來構(gòu)建加權(quán)基因共表達網(wǎng)絡(luò)，這可能導(dǎo)致一些在癌癥過程中起重要作用的基因的缺失. 因此，本研究對TCGA中胰腺癌數(shù)據(jù)進行嚴格的質(zhì)量控制和評估后，將表達數(shù)據(jù)中的4 916個基因全部納入WCGNA分析，共分為13個模塊. 接下來，本研究發(fā)現(xiàn)了與B細胞相關(guān)性最高的模塊. 此外，近期的研究表明，B細胞在免疫治療中發(fā)揮著非常重要的作用. JONSSON等[22]發(fā)現(xiàn)三級淋巴結(jié)構(gòu)中的B細胞與T淋巴細胞協(xié)同作用，最終可能靶向腫瘤細胞. 經(jīng)過免疫檢查點抑制劑治療后，有應(yīng)答者的B細胞受體的多樣性高于無應(yīng)答者，提示B細胞可能具有更好的識別腫瘤抗原的能力[23]. 因此，本研究將B細胞浸潤水平作為臨床信息. 結(jié)合基因表達數(shù)據(jù)和臨床信息，鑒定出9個與B細胞浸潤水平相關(guān)的關(guān)鍵模塊的樞紐基因，提示了這些基因參與胰腺癌進展的潛在機制. 在TIMER數(shù)據(jù)庫中分析這9個基因與B細胞的關(guān)系，均顯示出顯著的正相關(guān). 隨后進行差異表達基因分析和生存分析，結(jié)果顯示ITGAL表達最為顯著.ITGAL在腫瘤組織中高表達，當它在胰腺癌組織中高表達時，有著更好的預(yù)后.

ITGAL編碼淋巴細胞功能相關(guān)抗原-1 (LFA-1)的α亞基，表達于淋巴細胞表面[24]. 一方面，LFA-1主要與內(nèi)皮細胞上的胞內(nèi)黏附分子-1 (ICAM-1)相互作用，介導(dǎo)B細胞活化和產(chǎn)生免疫球蛋白，通過促進B細胞黏附來降低B細胞活化的閾值[25-26]. ICAM-1在腫瘤組織中的過表達，與腫瘤侵襲轉(zhuǎn)移、免疫功能調(diào)節(jié)有關(guān). 胰腺癌組織中ICAM-1的表達明顯高于正常組織，提示ICAM-1與胰腺癌的進展的相關(guān)性，ICAM-1可能作為胰腺癌早期診斷的生物標志[27-28]. 因此，ITGAL作為與ICAM-1相互作用的分子有可能是胰腺癌免疫治療的潛在靶點.