基于因果推斷方法檢驗mRNA在年齡和低級別膠質瘤中的中介效應*

王碧玨 郭志飛 楊海濤 李 治 王菊平 曹紅艷，6△ 周立業△

【提 要】 目的 構建年齡、mRNA表達水平和低級別膠質瘤(LGG)患者結局之間的因果中介推斷模型，識別介導年齡和LGG結局的中介因子，為LGG的治療和預后提供方法指導。方法 以TCGA數據庫中LGG數據為例，采用CIT和VanderWeele方法建立因果中介模型，識別年齡和LGG結局之間的介導基因。結果 識別出11個在年齡和LGG結局之間起中介作用的基因，中介間接效應的RR值均大于1，表明這11個基因的表達異常可能是LGG不良結局的危險因素。結論 基于因果推斷的方法能夠識別介導年齡和低級別膠質瘤結局的重要基因，為LGG的預后治療提供了潛在的預后標志物和治療靶點。

人腦膠質瘤多起源于神經膠質細胞，屬于神經上皮源性腫瘤，是最常見的原發性顱內惡性腫瘤，發病率(3～6.4)/10萬，分為低級別膠質瘤(low-grade gliomas，LGG)(WHOⅠ、Ⅱ級)和高級別膠質瘤(high-grade gliomas，HGG)(WHOⅢ、Ⅳ級)[1-2]。LGG約占膠質瘤的15%～25%，具有侵襲性生長、放化療抵抗等惡性生物學特征，并且有向高級別膠質瘤進展的傾向[3-4]。LGG的發生發展是基因、環境共同調控的結果，年齡≥40歲是LGG患者不良預后的重要因素，不可逆及持續增長的年齡如何影響疾病一直是研究的熱點，但并未得到一致的結論[5]。基因表達在特定蛋白的合成中調控細胞增殖與分化，參與LGG癌變的發生發展，是LGG預后和靶向治療的重要候選生物標志物[6]。年齡和基因表達均與LGG預后有關聯，故假設年齡、基因表達和LGG預后之間存在因果中介關系，即年齡可能調控基因表達影響LGG的預后。

因果推斷方法根據推斷理論的不同分為檢驗法和模型法，其中常用的檢驗法和模型法分別為因果推斷檢驗(causal inference test，CIT)[7]方法和VanderWeel[8]中介模型。本文以TCGA(The Cancer Genome Atlas)數據為例，聯合使用CIT因果中介方法和VanderWeele模型構建以年齡(以40歲為二分類的分類標準)為暴露、基因表達為中介、LGG患者預后為結局的因果推斷模型，探討年齡、mRNA表達水平和LGG患者結局之間的關系，識別中介因子。所識別的中介基因將為LGG患者的靶向治療提供理論依據，有利于改善患者預后，提高生命質量。

材料與方法

1.數據資料來源與整理

使用TCGA Assembler軟件在TCGA網站下載LGG數據[9-10]，從中選取包含臨床和mRNA表達數據的患者共515例，刪除零表達值所占比例≥30%的基因，最終mRNA有17236個變量，結局為截止至隨訪時間是否死亡。

2.方法

(1)CIT

Millstein(2009)提出了判斷因果關系的CIT方法，將因果推斷問題視為一條數學條件的鏈，所謂“因果關系”，即中介的差異在一定比例上決定了結局的某些差異，即使這個比例很小。CIT方法的適用條件為：①連續型或離散型暴露變量，且不限制為單變量，也可是一組變量；②連續型潛在中間變量；③連續型或二分類結局變量[11-12]。CIT方法的基本思想是，引入暴露L，潛在中間變量G和結局變量T，三者需要滿足以下四個條件：①L與T相關，②G與T|L相關，③L與G|T相關，④L與T|G相互獨立。以年齡為L，基因表達為G，生存結局為T，四個條件對應的模型如下：

(1)

(2)

(3)

(4)

H0：γ1=0，H1：γ1≠0，

H0：βi=0，H1：βi≠0，

H0：γ3=0，H1：γ3≠0，

H0：γ4≠0，H1：γ4=0.

CIT方法為基于似然的假設檢驗方法，若以上四個條件均滿足，則推斷潛在中間變量是結局的原因，即對每一個條件計算P值，將所有P值中的最大值作為整體檢驗的結果。

(2)VanderWeele模型

(5)

(6)

將總效應分解為限制性直接效應(controlled direct effects，CDE)、自然直接效應(natural direct effects，NDE)和自然間接效應(natural indirect effects，NIE)。通過計算直接效應和間接效應的RR值，進一步估計處理因素對結局的直接效應和間接效應。間接效應的計算見公式(7)，利用bootstrap法重復1000次計算RR值95%置信區間。

RRNIE=exp(θm×βx+θinteraction×βx)

(7)

(3)統計方法實現

應用R軟件進行統計分析，CIT第四步在cit包實現，bootstrap在boot包實現，顯著水平α=0.05。

結 果

1.CIT方法分析結果

(1)第一步分析結果

第一步校正性別協變量后，研究年齡與LGG的相關性，P值為0.000316，RR值為0.124，95%CI為(0.060，0.192)，表明年齡與LGG有關聯，即年齡對LGG患者的結局有直接效應。

(2)第二步分析結果

第二步在校正性別協變量及控制年齡基礎上，研究LGG與每個基因表達水平的相關性。以每個基因的表達為自變量、LGG為應變量擬合logistic回歸模型(2)，對檢驗水準α進行Bonferroni校正后，有557個(3.23%)基因與LGG結局有關聯。圖1為第二步驟的Mahattan圖，橫坐標表示基因在23條染色體上的位置，縱坐標為-log(P)尺度。所得有意義的基因將進入到第三步的分析。

(3)第三步分析結果

第三步在校正性別協變量及控制LGG結局的基礎上，研究每個基因的表達水平與年齡之間的相關性。具體過程是以每個基因的mRNA表達水平為應變量、年齡為自變量來擬合線性回歸模型(3)，用所得P值表示第j個基因的表達水平與年齡的關聯性，對檢驗水準α進行Bonferroni校正后，有25個(4.49%)基因與年齡有關聯。圖2為第三步驟的Mahattan圖。所得有意義的基因將進入到第四步的分析中。

(4)第四步分析結果

第四步在校正性別協變量及控制基因表達水平基礎上，研究LGG與年齡之間是否相互獨立。具體過程是以LGG為應變量、年齡為自變量來擬合線性回歸模型(4)。該步驟利用R軟件的“cit”包進行分析，可得到P值和基于置換的FDR值(q值)，q值小于0.05表示控制基因的表達水平后LGG與年齡之間相互獨立。有11個(44%)基因的表達水平有統計學意義，分別為ARL3、DPY19L1、E2F7、FAM171A1、IGFBP2、LAMA2、MRPS16、PER3、PLAT、RANBP17和RASL10A。因此年齡可能通過調節這11個基因的表達水平來影響LGG結局。11個基因的位置信息和CIT方法四個條件的P值見表1。

圖1 17236個基因的Mahattan圖

圖2 557個基因的Mahattan圖

表1 CIT方法分析結果及基因的位置信息

2.VanderWeele分析結果

CIT方法識別出11個在年齡和LGG間起中介作用的基因，采用VanderWeele模型對中介的間接效應(NIE)進行估計，并利用bootstrap方法計算RR值的95%置信區間。圖3是11個基因的表達水平分別在年齡和LGG結局中的間接效應的RR值。

結果表明11個基因在年齡和LGG之間的效應均有意義，且均大于1，說明這11個基因的表達異常可能是LGG不良結局的危險因素。

圖3 11個基因的mRNA表達水平分別在年齡和LGG之間的間接效應大小

討 論

本研究基于因果推斷方法檢驗mRNA在年齡和LGG之間的中介效應，年齡通過調控11個基因的表達，進而影響LGG的預后，為LGG患者的治療提供了潛在的預后標志物和治療靶點。

ARL3在膠質瘤中表達下調[13]，其通過參與免疫微環境和血管生成，進而影響膠質瘤的生長、侵襲和轉移。DPY19L1在膠質瘤中表達下調，有可能通過調節大腦皮層中谷氨酸能神經元的徑向遷移而作用于膠質瘤[14-15]。LAMA2基因編碼層粘連蛋白α2鏈，是一種重要的細胞外基質蛋白，有研究發現層粘連蛋白可能通過促進腦室下區少突膠質細胞祖細胞的存活來調節產后少突膠質細胞的產生[16]。因此，這些基因可能是膠質瘤的預后標志物和治療靶點。E2F7、PER3和RANBP17基因均為與細胞增殖、細胞周期有關的因子，其中，E2F7在膠質瘤中高表達，可能通過與CDK1的協同作用參與膠質瘤的發生發展[17]。E2F7的過表達可能通過上調Beclin-1誘導自噬而導致LGG的不良預后和結局[18]，因此提示E2F7可能是神經膠質瘤的不良預后指標。PER3是人類Period基因家族中最重要的一員，可通過調控生物節律等重要生物通路而參與腫瘤的發生[19]。目前已發現PER3在乳腺癌[20]、頭頸部鱗狀細胞癌[21]等癌癥中表達異常，且與腫瘤浸潤程度、復發率和預后有較大關系。但其在LGG中的作用仍需要進一步的驗證。RANBP17與細胞外基質和細胞黏附功能相關，是膠質瘤獨立的預后因素[22]，且RANBP17的表達水平與年齡呈負相關[23]，因此可以推斷RANBP17可能是在年齡和膠質瘤之間發揮作用的中介因子。IGFBP2和RASL10A為抑癌基因。IGFBP2在膠質瘤中通常為過表達[24]，且IGFBP2的瞬時表達對膠質瘤進展和生存有顯著影響，抑制IGFBP2也能抑制腫瘤細胞的擴散[25]。LGG中的IDH突變促進腫瘤抑制因子的信號傳導，同時抑制IGFBP2來改善患者的生存[26]，提示IGFBP2是膠質瘤進展的驅動因素和預后因素，是治療膠質瘤的可行治療靶點。RASL10A又稱RRP22，其具有腫瘤抑制作用[27]，可以抑制腫瘤細胞生長，促進非依賴性細胞死亡，降低侵襲性，抑制軟瓊脂中膠質瘤細胞生長，從而起到抑癌功能[28]。

在11個基因中，FAM171A1、PLAT和MRPS16與膠質瘤的關系尚不清楚，但有研究發現他們與其他癌癥有關。FAM171A1是一種調節人體細胞形狀和侵襲性生長的物質，其表達上調可能會促進癌細胞的侵襲性增長[29]。PLAT基因在乳腺癌和卵巢癌中發揮作用[30]。MRPS16是已知的與胚胎或未產期致死的基因[31]，截斷型MRPS16變體可能是胎兒死亡的潛在原因。本研究通過嚴格的因果推斷方法檢驗了其在年齡和LGG結局間的影響，可為LGG今后的研究提供新的生物學假設。

本研究僅控制了性別協變量，CIT方法在處理混雜因素上具有穩健性，可彌補潛在混雜因素干擾的不足。VanderWeele中介模型在暴露和中介變量存在交互作用時，其效應估計也較為穩健。本文未研究年齡和基因表達間的交互作用，今后將進一步探討。LGG預后隨訪數據的結局變量包括生存時間和是否死亡，但目前未見針對生存分析的CIT方法，今后將進一步研究生存分析的CIT方法，提高模型的識別能力。

綜上所述，本研究基于因果推斷的方法識別出介導年齡和LGG患者結局的重要基因，年齡通過調控其基因表達水平進而對LGG結局產生影響，為膠質瘤患者的治療提供了重要的方法參考。