基于mRNA生物信息學分析的肺腺癌免疫基因預后風險模型建立與評估*

葛美玲，胡月，丁杰，劉艷紅，高玒

(南京大學醫學院附屬鼓樓醫院，a.生物樣本庫，b.病理科，南京 210008)

在過去的數十年里，肺癌發病率和死亡率逐年升高，目前已經成為癌癥相關死亡的主要原因之一。據WHO最新數據顯示，肺癌的發病率(11.6%)和死亡率(18.4%)在所有惡性腫瘤中位居第一[1]。在我國，按發病人數順位排序，肺癌也位于惡性腫瘤發病首位[2]。肺腺癌(LUAD)作為非小細胞肺癌(NSCLC)最常見的組織學亞型，占肺癌發病率的40%。越來越多的證據表明，LUAD發生和發展與腫瘤免疫逃逸密切相關。腫瘤細胞不僅能夠通過表達抑制性配體與細胞毒性T細胞表面受體結合抑制其功能，還能通過分泌免疫抑制因子，調節輔助淋巴細胞降低機體的抗腫瘤反應[3]。因此，尋找特異性的免疫分子標志物，開發有效的腫瘤免疫監控手段，探索腫瘤免疫治療的新方案，對于LUAD患者的治療及預后具有重要的意義[4]。本研究從分子水平出發，篩選有效的免疫基因標志物，構建預后風險模型，探索免疫基因在LUAD診斷、治療及預后中的重要作用，并為尋找新的檢測手段和治療靶標提供依據。

1 資料與方法

1.1患者及臨床樣本

1.1.1患者入組 隨機選取2019年5月至9月在南京鼓樓醫院進行手術的38例早期LUAD患者。患者手術前均未接受消融、化療或放療。患者人口學資料：男性23例(60.53%)，女性15例(39.47%)，年齡(61.29±10.06)歲。收集上述患者的腫瘤及正常石蠟組織樣本76份。

1.1.2病理學評估 將1.1中的石蠟組織制作病理切片，由2位病理科醫師檢查所有HE染色切片，按照2017年IASLC(國際肺癌研究協會)第八版TNM分期標準及2015年WHO肺腺癌分型標準(貼壁型、腺泡型、乳頭型、微乳頭型、實體型)對樣本進行分期及分型。鑒于肺腺癌的高度異質性，以腫瘤中百分比最高的成分判定主要組織學分型。患者臨床病理資料如下：原位肺腺癌3例(7.89%)，浸潤性肺腺癌35例(92.11%)；浸潤性肺癌病理分型：貼壁型20例(57.14%)，腺泡型10例(28.57%)，乳頭型3例(8.57%)，實體型2例(5.71%)。腫瘤分期：0期3例(7.89%)，Ⅰ期35例(92.11%)。腫瘤分化程度：低分化2例(5.26%)，中分化17例(44.74%)，高分化19例(50.00%)。

1.1.3臨床樣本全轉錄組mRNA高通量測序 對1.1.1中的腫瘤組織及其對應的正常肺組織的石蠟標本進行全轉錄組mRNA高通量測序。樣本符合以下要求：(1)病理分期為0期或Ⅰ期；(2)FFPE樣本的腫瘤陽性細胞比例>80%。采用測序數據的FPKM值衡量基因的表達量并進行后續分析。

1.2訓練集及測試集數據

1.2.1全轉錄組mRNA數據下載及預處理 從TCGA (https://portal.gdc.cancer.gov/)中下載肺腺癌(LUAD)的mRNA測序數據(HTSeq-FPKM)及臨床數據作為訓練集進行后續分析。在GEO平臺中下載數據集GSE101929、GSE50081、GSE41271、GSE42127 mRNA測序數據及臨床數據作為測試集，使用svaR軟件包矯正各數據集測序數據的多批次效應，將得到的結果用于后續分析。

1.2.2患者臨床數據篩選 選取截至2020年6月收錄于TCGA數據庫以及GEO數據庫中GSE101929、GSE50081、GSE41271、GSE42127數據集中的LUAD患者。納入標準：(1)病理診斷為LUAD，且無其他惡性腫瘤；(2)mRNA測序數據及隨訪資料完整，隨訪天數大于30 d。(3)同例患者存在不同類型腫瘤樣本時，僅選擇1個腫瘤樣本的數據。不同類型樣本選擇順序：冰凍組織>石蠟組織。本研究中訓練集共納入患者439例，測試集共納入患者476例。

1.3肺腺癌預后模型的建立及驗證

1.3.1訓練集中腫瘤及正常樣本差異表達基因分析 對訓練集mRNA測序數據進行預處理，利用limma R[5]軟件包從訓練集中鑒定了439例腫瘤組織樣本及52例正常肺組織樣本之間的差異表達基因(differently expressed genes，DEGs)。以log2-fold change (|FC|)>1和錯誤發現率(false discovery rate，FDR)<0.05作為閾值。

1.3.2預后模型的建立及驗證 將腫瘤與正常樣本的DEGs與ImmPort數據庫(https://www.immport.org/home)提供的免疫相關基因取交集篩選出候選基因。采用單變量COX回歸對候選基因與患者總生存時間(OS)進行分析，篩選P<0.05的基因。然后采用Lasso回歸篩選上述基因并構建模型[6-7]，計算調整參數lambda值(tuning parameter，λ)。通過識別與交叉檢驗誤差均值最小值對應的λ值作為依據，篩選變量計算回歸系數，構建預后風險模型[8]。由回歸模型得到的回歸系數乘以其mRNA水平產生了每個樣本的預后風險評分(risk score)。選擇訓練集的風險評分中位值作為閾值進行分組，將風險分值大于閾值的患者歸為高危組(high risk group)，分值小于或等于閾值的患者歸為低危組(low risk group)。采用相同的公式和閾值計算測試集中患者的風險評分并進行分組。分別繪制2個隊列患者的生存曲線并計算高、低危組患者OS的差異。分別繪制兩隊列患者的總生存時間的ROC曲線并計算曲線下面積(area under curve，AUCROC)，評估風險模型的預測性能[9]。

1.4統計學分析 所有統計分析均使用R軟件(4.0.2)進行。使用Kaplan-Meier法進行生存分析并繪制曲線。采用Log-rank檢驗(Log-rank)分析生存差異的統計學意義。采用LASSO回歸構建預后風險模型。通過ROC曲線評價預后模型的效能。采用卡方檢驗(Chi-square test)或Fisher精確檢驗(Fisher′s exact test)分析高、低危組間的臨床病理資料。所有統計檢驗均為雙側檢驗，以P<0.05為差異有統計學意義。

2 結果

2.1基于訓練集的肺腺癌預后模型構建 對訓練集數據中腫瘤與正常組織mRNA測序數據進行比較，篩選出12 946個DEGs進行下一步分析。其中，與正常組織相比，10 207個DEGs在腫瘤組織中被上調，2 739個DEGs被下調。將這些DEGs與免疫基因進行交集篩選出1 091個基因。采用單因素COX 回歸分析這些基因與OS的相關性，篩選出86個候選基因(P<0.05)。將上述候選基因納入LASSO回歸分析，并計算最優λ值。根據最優λ值重新擬合數據，篩選出12個基因包括：纖維母細胞生長因子2(fibroblast growth factor 2,FGF2)、軸突導向因子7A(semaphorin 7A，SEMA7A)、抑制素亞基α (inhibin subunit alpha,INHA)、葡萄糖-6-磷酸異構酶 (glucose-6-phosphate isomerase,GPI)、血管生成素樣蛋白4(angiopoietin like 4,ANGPTL4)、腫瘤壞死因子受體超家族成員11 A (TNF receptor superfamily member 11a ,TNFRSF11A)、細胞色素b-245 β鏈 (cytochrome b-245 beta chain,CYBB)、精氨酸酶2 (arginase 2,ARG2)、尾加壓素2 (urotensin 2,UTS2)、白血病抑制因子受體亞基α (LIF receptor subunit alpha,LIFR)、SHC適配器蛋白3(SHC adaptor protein 3,SHC3)、細胞毒性T淋巴球相關蛋白4(cytotoxic T-lymphocyte associated protein 4, CTLA4)，然后計算LASSO回歸系數構建公式(圖1A、B)。公式如下：風險評分= (-0.001 77×CYBBmRNA水平)+(0.170 269×FGF2mRNA水平)+(-0.045 58×ARG2mRNA水平)+(0.0051 21×SEMA7AmRNA水平)+(0.005 625×GPImRNA水平)+(0.003 913×INHAmRNA水平)+(-0.133 98×UTS2mRNA水平)+(0.043 29×ANGPTL4mRNA水平)+(-0.016 18×LIFRmRNA水平)+(0.144 495×TNFRSF11AmRNA水平)+(-0.083 3×SHC3mRNA水平)+(-0.043 53×CTLA4mRNA水平)。

注：A，Lasso回歸分析中的調節參數(λ)與模型誤差圖。橫坐標為λ的對數值，縱坐標為模型誤差，圖形最上方為模型篩選出來的對應λ的變量個數。以最小二項偏差篩選出12個變量； B，Lasso 回歸系數路徑圖。圖中的每條線代表了1個變量，縱軸為回歸系數，橫坐標為λ的對數值。10倍交叉驗證法對log(λ)進行測試篩選出12個非零的變量。

LASSO回歸中FGF2、SEMA7A、INHA、GPI、ANGPTL4、TNFRSF11A的正系數提示其高表達與較差的OS有關，而CYBB、ARG2、UTS2、LIFR、SHC3、CTLA4的負系數提示其高表達與較好的OS有關。根據中位風險評分(0.334)，將訓練集中所有患者分為高危組(221例)和低危組(218例)。生存分析顯示低危組患者3年、5年生存率分別為80.2%和62.5%，高危組患者3年、5年生存率分別為44.2%和23.7%。高危組的總生存率明顯少于低危組患者(P<0.001,圖2A～C)。將所有患者的12個基因表達譜可視化為熱圖(圖2D)。訓練集的預后模型AUCROC為0.759 (圖2E)。表明該模型在生存監測中具有良好的效力。基于自變量分析建立了預測1年、2年、3年生存率的列線圖，年齡、病理分期(stage)、腫瘤大小分期(T)、淋巴結轉移(N)、遠處轉移(M)及風險評分(risk score)均作為參數納入列線圖(圖2F)。

注：A，高危組及低危組患者OS比較； B，高危組及低危組患者的風險曲線； C，高危組及低危組患者的生存狀況； D，高危組及低危組患者12個模型基因表達熱圖； E，以訓練集中患者年齡、性別、病理分期、腫瘤大小、淋巴結轉移、遠處轉移、風險評分分別作為變量的生存預測模型的ROC曲線； F，以訓練集中患者性別、病理分期、腫瘤大小、淋巴結轉移、遠處轉移、風險評分分別作為變量預測患者1年、2年、3年生存率的列線圖。

2.2基于測試集的肺腺癌預后模型驗證 在測試集中評估預后風險模型的預后預測性能，使用上述風險評分公式計算測試集患者的風險評分，使用與訓練集患者風險評分相同的閾值按同樣的方式將測試集劃分為高危組(122例)和低危組(354例)。在兩組人群的生存分析中，低危組患者3年、5年總生存率分別為72.8%和59.3%，高危組患者3年、5年總生存率分別為53.7%和42.5%，高危組的總生存率明顯低于低危組(P<0.001,圖3A～C)。將測試集中所有患者的12個模型基因表達譜可視化為熱圖(圖3D)。模型 AUCROC為 0.707(圖3E)。結果表明，預后模型對預測測試集患者生存預后的表現較好。基于自變量分析建立了預測患者1年、2年、3年生存率的列線圖，年齡、病理分期、風險評分均作為參數被納入列線圖 (圖3F)。

注：A，高危組及低危組患者OS比較； B，高危組及低危組患者的風險曲線；C，高危組及低危組患者的生存狀況；D，高危組及低危組患者中12個基因表達熱圖； E，以測試集中患者年齡、性別、病理分期、風險評分分別作為變量的預測模型的ROC曲線； F，以測試集中患者年齡、性別、病理分期、風險評分分別作為變量預測患者1年、2年、3年生存率的列線圖。

2.3預后風險模型對早期肺腺癌患者預后的預測潛能 對訓練集及測試集中腫瘤早期(Ⅰ期)患者的總生存時間與風險評分進行分析，訓練集的AUCROC為 0.709(圖4A)，測試集的 AUCROC為 0.639(圖4B)。表明該模型對早期LUAD的預后也有較好的預測潛能。對訓練集患者的風險評分與各項臨床病理數據進行分析，該評分隨著患者Stage(圖5A)和T分期(圖5B)進展逐漸升高，Ⅰ期患者腫瘤樣本的風險評分明顯高于正常組織(P<0.05)，淋巴結轉移患者的得分明顯高于未轉移患者(P<0.05,圖5C)，但遠處轉移患者與未轉移患者之間的風險評分差異無統計學意義(P>0.05,圖5D)。

注：A，訓練集中患者年齡、性別、病理分期、腫瘤大小、淋巴結轉移、遠處轉移及風險評分分別作為變量的生存預測模型的ROC曲線；B，以測試集中早期LUAD(Ⅰ期)患者年齡、性別、病理分期、風險評分分別作為變量的生存預測模型ROC曲線。

注：A，正常樣本及腫瘤Ⅰ、Ⅱ、Ⅲ、Ⅳ期樣本的風險評分比較； B，不同腫瘤大小分期(T1、T2、T3、T4)患者的風險評分比較； C，未發生淋巴結轉移及淋巴結轉移患者的風險評分比較； D，未發生遠處轉移及遠處轉移患者的風險評分比較。

2.4預后模型與臨床早期肺腺癌患者病理參數的相關性 本研究對臨床早期肺腺癌患者的轉錄組測序數據進行了分析，對比38例患者的腫瘤與正常組織mRNA表達情況(圖6)，CYBB、FGF2、LIFR、SHC3基因在正常組織中的表達量高于腫瘤組織，差異具有統計學意義(P<0.05)；GPI、UTS2、TNFRSF11A、CTLA4基因在正常組織中的表達量低于腫瘤組織，差異有統計學意義(P<0.05)。依據上述構建的預后風險模型閾值(0.334)，將早期肺腺癌患者分為高危組和低危組。高危組患者12例，低危組患者26例。比較不同預后模型分組患者的臨床病理參數差異(表1)。兩組患者在性別、年齡、吸煙史、分化程度、TNM分期及腫瘤大小分期方面差異無統計學意義(P>0.05)，但高危組患者的腫瘤最大徑大于低危組患者，該組乳頭型或實體型肺腺癌樣本比例也高于低危組，差異均有統計學意義(P=0.029，P=0.044)。

表1 高、低危組患者各項臨床病理參數比較

注：縱坐標為基因表達量，取測序數據FPKM值的對數值(Log2)；橫坐標為12個基因的基因名稱。紅色圖標代表腫瘤組織基因表達量，藍色圖標代表正常組織基因表達量。

3 討論

大量的研究表明，免疫微環境在肺腺癌的發生、發展中起著重要的作用。探索異常表達的免疫基因對于確定腫瘤發生、發展機制，發現潛在治療靶點具有重要意義。對于LUAD這種異質性很高的腫瘤而言，單一功能的基因模型很難完整評估所有患者的免疫狀態。本研究基于不同來源的數據集，構建了包含不同功能的免疫差異基因的預后模型，試圖開發一種方便可靠的方法來監測LUAD患者的免疫狀態并評估預后。

本研究建立的模型篩選了12個免疫基因。這些基因涉及不同的蛋白質、通路及功能：CTLA4主要表達在活化的 T 細胞及Treg細胞表面，是第1個被證實的靶向的免疫檢查點分子[10]，并在抑制機體免疫反應及促進腫瘤免疫逃逸中起著重要作用[11]。FGF2是一種多肽類生長因子，可以激活Ras-Raf-MapK、PI3K/Akt/mTOR等信號通路，促進腫瘤細胞生長、增殖、分化及轉移等行為[12]。UTS2是一種生長抑素樣環狀多肽，能夠促進腫瘤血管形成，該基因促新生血管形成的作用與FGF2相似[13]。LIFR是白血病抑制因子的受體，能夠介導由白細胞介素6相關因子引起的信號參與免疫調節，并參與PI3K/Akt通路從而影響LUAD發生、發展[14]；TNFRSF11A作為腫瘤壞死因子超家族細胞因子，能夠抑制腫瘤細胞的能動性和遷移[15]。GPI是細胞內糖酵解和糖異生過程中的關鍵酶，以自分泌的形式與細胞自身或鄰近細胞表面受體相結合[16]，影響肺腫瘤細胞的凋亡、遷移及浸潤等行為[17]。ARG2可能通過阻礙聚乙二醇化的ARG1與配體結合來抑制后者的抗腫瘤作用[18]。ANGPTL4能夠促進腫瘤細胞肺轉移[19]；SHC3被認為是T細胞受體和B細胞受體信號傳導的負調節因子[20]；CYBB作為NOX2的亞基，其表達缺失能夠減少腫瘤細胞的肺轉移[21]。INHA與受體TGFβRⅢ結合能夠拮抗激活素，具備抑制多種細胞腫瘤化的效應[22]。信號素蛋白SEMA7A通過調控mTOR信號通路，并被EGFR通路誘導表達，促進肺腺癌發生，可以作為EGFR突變的肺腺癌預測生物標志物及新治療靶點[23]。

為了進一步探索預后模型在早期肺腺癌中的應用前景，筆者分析了訓練集中不同病理分組患者的風險評分，結果發現該評分在早期腫瘤組織中顯著高于正常組織，具有成為早期肺腺癌診斷標志物的潛力。生存分析也揭示預后模型對于早期腫瘤患者的生存時間具有一定的預測能力。同時，筆者依照模型基因及參數結合測序數據對38例臨床早期肺腺癌患者樣本進行了驗證，發現與正常肺組織相比，CYBB、FGF2、LIFR、SHC3基因在早癌組織中表達降低，GPI、UTS2、TNFRSF11A、CTLA4基因在早癌組織中表達升高。隨后，將患者進行了高、低危分組，并評估了兩組患者的臨床病理參數的差異，結果發現腫瘤最大徑及病理分型在高、低危組患者之間存在差異。高風險評分與更大的腫瘤直徑以及特異的病理分型(乳頭或實體型)有關。已有研究表明[24]，相對于貼壁型肺腺癌患者，實體型為主的肺腺癌患者5年生存率低，預后較差。因此，筆者認為具有高風險評分的早期肺腺癌患者可能預后更差。

綜上所述，本研究構建了1個免疫差異基因的預后風險模型，該模型對于肺腺癌患者預后具有良好的預測能力，相關免疫基因可能成為早期肺腺癌診斷的生物標志物。然而，本研究尚存在許多局限性。本研究構建的模型是建立在生物信息及統計學分析基礎上的，篩選出的模型基因及其對下游信號通路可能產生的影響還需要生物學實驗的進一步驗證。雖然筆者對該模型進行了臨床驗證，發現其對于早期肺腺癌預后有一定的預測能力，但本研究選擇的樣本量較少而且患者的生存及預后數據尚不完善，這種預測能力的可靠性還需進一步探索。我們希望本研究結果能為肺腺癌患者早期診斷、風險分層、預后評估及免疫治療提供新的研究方向。