基于TCGA篩選結腸癌預后相關lncRNA及建立預后風險模型

何天 曹天生

廣州市花都區人民醫院胃腸外科，廣州 510800

結腸癌是常見的消化道腫瘤，據全球癌癥統計數據，其在惡性腫瘤中發病率居第5，隨著飲食結構的改變，結腸癌發病率逐漸升高［1-2］。由于結腸癌發病隱匿，疾病初期無明顯臨床癥狀，導致多數患者就診時已處于腫瘤的中晚期，因此大多失去標準治療的最佳時機。目前，手術切除仍是治愈的主要方式，而對有轉移的晚期患者，其5 年生存率不足20%［3］。結腸癌患者預后判斷主要依賴TNM 分期，血清學檢查中癌胚抗原、CA19-9等生物標記物對預后評價具有重要的參考作用，但敏感性和特異性欠佳［4］。因此，臨床上亟需新的生物學標志分子和有效的預后評估模型以提高對結腸癌患者診治的準確性。

長鏈非編碼RNA（long non-coding RNA，lncRNA）指大于200 個核苷酸的非編碼RNA，一開始被認定是轉錄“垃圾”，不具有調節機體過程等生物學功能［5］。隨著生物信息學技術與芯片技術的發展，發現lncRNA 在多種遺傳生物學過程中，像轉錄前、后的調控和表觀遺傳調控等，都發揮著重要的作用［6］。研究證實當lncRNA 轉錄發生改變時可促進癌癥的發生和轉移［7］。關于lncRNA 如何影響腫瘤進程，哈佛醫學院Pandolfi 課題組提出的競爭性內源RNA（competing endogenous RNA，ceRNA）假說認為，lncRNA 可通過結合微小RNA（microRNA，miRNA）來阻隔miRNA 與其他信使RNA（mRNA）作用，從而降低了miRNA 對mRNA 的抑制調控，促進相應mRNA 的表達，進而影響疾病的進程［8］。已有多種異常表達lncRNA 被證實作為ceRNA 與結腸癌進展相關，如Peng 等［9］發現lncRNA HOTAIR 在結腸癌中的表達明顯升高，作為ceRNA 調節miR-34a 促進結腸癌的進展，與結腸癌的遠處轉移和不良預后有關。為構建1 個由lncRNA 組成的結腸癌預后風險評估模型，本研究基于癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數據庫，獲取預后相關lncRNA，構建可用于評估結腸癌患者預后的多基因風險模型，并初步探索lncRNA 參與結腸癌進程的機制。

材料與方法

1、原始數據下載及數據預處理

數據提取時間：建庫至2022年3月1日。從TCGA 上下載結腸癌樣本（癌和癌旁組織）轉錄組測序數據及臨床資料。使用Rstudio 軟件對獲取的數據進行初步處理及篩選，根據以下標準排除部分患者：沒有生存狀態和生存時間信息的患者；缺乏完整lncRNA 表達數據的患者；對重復樣本的患者隨機去重。基因表達數據中，去除表達量較低的基因，通過從全球數據中心（Global Data Center，GDC）數據庫下載人類基因注釋文件（genecode.v22）進行注釋。Genecode 數據庫中有關于lncRNA 的種類說明，依此篩選過濾得出lncRNA的表達矩陣。

2、差異表達基因分析

依據癌旁樣本挑選匹配的癌癥樣本構建配對樣本lncRNA 表達矩陣，利用“edgeR”R 包篩選，閾值設置為：|log 2 fold change|>1.5 和P-value <0.05，獲得差異表達lncRNA（differentially expressed lncRNA，DElncRNA）。根據分組繪制lncRNA 主成分分析（principal component analysis，PCA）圖觀察配對樣本間lncRNA 的表達差異，使用“pheatmap”包繪制表達熱圖、“ggplot2”包繪制火山圖可視化上調及下調lncRNA。

3、挑選預后相關lncRNA

構建癌組織樣本DElncRNA 表達矩陣，與患者生存資料合并，使用統計學方法進一步篩選。先采用COX 回歸模型單變量分析評估DElncRNA 表達水平與患者總生存期間的關系，選取明顯影響生存期的lncRNA（P<0.05）。隨后利用“glmnet”R 包進行Lasso 回歸分析，過濾基因間的共線性因素，根據參數中的Lambda 值，選擇Lambda.min 作為臨界點，進行篩選。再以lncRNA 表達值的中位數作為分界值，將患者分位高低表達組，進行K-M 生存分析，挑選出表達量與生存率明顯相關的lncRNA（P<0.05）。最后使用“survival”R 包的逐步回歸法構建多元COX 回歸模型，衡量各個lncRNA 表達情況對生存的影響程度，用風險比（hazard ratio，HR）體現，以P<0.05 作為標準，得出預后相關lncRNA［10］。

4、構建預后風險模型

根據多因素COX 回歸分析中計算得出預后相關lncRNA 的回歸系數（β），按照如下公式構建結腸癌預后風險評估模型［11］。β 表示相對危險度，是對模型中各個因素相對危險程度的量化。當β>0時，對應因素的取值越大，患者死亡的風險越高；β<0 時，對應因素取值越大，患者死亡的風險越低。

公式中Risk score（RS）為風險評分，n 為關鍵lnRNA 基因數量，Exp為lncRNA表達量，β為回歸系數。

5、評估預后風險模型

先計算模型的C 指數值，其可用來評價模型的預測能力，主要用于計算COX 回歸模型中預測值與真實之間的區分度［12］。一般認為，如C指數介于0.50～0.70之間則為低準確度；在0.71～0.90 之間則視為中等準確度；高于0.90 則可看做高準確度。然后對此模型構建時間依賴的受試者工作特征曲線（receiver operating characteristic curve，ROC），并計算ROC 曲線下的面積（area under the curve，AUC），以評估模型預測患者3 年、5 年生存期的能力。AUC 和C 指數意義相似，一般認為0.85

6、構 建ceRNA 網絡及基因本體論（Gene Ontology，GO）、京都基因與基因組大百科全書數據庫（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析

進一步構建lncRNA-miRNA-mRNA 相互作用ceRNA網絡。先利用LncBOOK 在線數據庫，根據相互作用得分，選擇前30個靶miRNA。然后通過數據庫miRDB、TargeScan和miRTarBase 預測miRNA 下游的靶mRNA，取三者交集。在TCGA 下載結腸癌miRNA-seq 表達數據，與RNA-seq 進行整合，獲得一個含有lncRNA、miRNA、mRNA 的表達數據。計算miRNA 與靶mRNA 間的皮爾遜相關系數（Pearson correlation coefficient，PCC），得出與miRNA 表達明顯相關（|PCC|>0.4，P<0.05）的靶mRNA。去掉部分沒有靶mRNA的miRNA，構成ceRNA 網絡，使用Cytoscape 軟件構圖以可視化。GO 分析從生物過程（biological process，BP）、細胞組分（cellular component，CC）、分子功能（molecular function，MF）3 方面注釋一系列生物學過程［14］。KEGG 富集分析可以將基因組信息與高階的功能信息聯系起來［15］。依據匹配到的mRNA，使用“ClusterProfiler”R包，進行GO、KEGG富集分析探索預后相關lncRNA 影響結腸癌進展可能的機制。整個結腸癌TCGA數據提取及分析流程見圖1。

圖1 結腸癌TCGA數據提取分析流程圖

結果

1、TCGA中結腸癌患者臨床病理信息

從TCGA 共下載514 個結腸癌樣本轉錄組數據和452 個患者臨床病理信息。依據排除標準得到364 個樣本（癌旁樣本28 個）和336 例患者信息，對其臨床病理信息進行統計（表1）。其中，男性患者占53.6%，女性患者占46.4%；患者年齡31～90 歲，平均65.9 歲；生存或隨訪時間0.03～137.40個月，中位生存時間為5.62個月。

表1 336例結腸癌患者的臨床病理信息

2、配對樣本差異分析結果

構建28 對癌與癌旁組織配對的lncRNA 表達矩陣，以|log2 FC|>1.5 和P<0.05 為標準，共獲得868 個DElncRNA，繪制PCA 圖（圖2）、表達熱圖（圖3）及火山圖（圖4）。其中PCA 圖顯示DElncRNA 在癌與癌旁組織中的表達具有明顯差異；熱圖顯示出在癌組織樣本表達量增加的上調lncRNA及表達量下降的下調lncRNA；火山圖顯示上調lncRNA 548個、下調lncRNA 320個。

圖2 配對樣本長鏈非編碼RNA表達主成分分析圖

圖3 配對樣本長鏈非編碼RNA表達熱圖

圖4 長鏈非編碼RNA火山圖

3、結腸癌預后關鍵lncRNA的篩選

構建癌組織DElncRNA 表達矩陣，先進行單因素COX回歸分析，篩選出40 個lncRNA。隨后進行Lasso 回歸模型分析，得到34 個候選lncRNA。以lncRNA 表達值的中位數作為分界值，進行高低分組生存分析，發現其中的14 個lncRNA 表達譜與生存時間明顯相關（P<0.05），繪制K-M 生存曲線（圖5）。

圖5 14 個長鏈非編碼RNA 高低表達組間生存分析：ELFN1-AS1（A）、RP5-884M 6.1（B）、LINC00461（C）、RP11-161I 6.2（D）、LINC01132（E）、RP11-108K 3.1（F）、RP1-170O19.17（G）、RP11-108K 3.2（H）、AC114730.3（I）、RP1-79C4.4（J）、RP4-816N1.7（K）、RP3-380B8.4（L）、RP5-823G15.5（M）、WFDC21P（N）

最后進行多元COX 回歸模型分析，得到7 個預后相關lncRNA（ELFN1-AS1、RP5-884M6.1、LINC00461、LINC01132、RP1-79C4.4、RP4-816N1.7、RP3-380B8.4）。計算得出每個lncRNA 的HR 值及95%可信區間，繪制風險森林圖（圖6）。其中1 個lncRNA（LINC01132）的HR<1，提示其表達量增加可降低結腸癌不良預后的風險，是患者的保護因素；6 個 lncRNA（ELFN1-AS1、RP5-884M6.1、LINC00461、RP1-79C4.4、RP4-816N1.7、RP3-380B8.4）的HR>l，提示它們表達量增加可增加結腸癌不良預后的風險，是患者的危險因素。

4、預后風險模型的構建

根據多元COX 回歸模型，提取每個lncRNA 的回歸系數，依據公式構建由7 個lncRNA 組成的結腸癌預后風險模型。

Risk score=0.256 8×（ELFN1-AS1 表達量）+0.289 1×（RP5-884M6.1 表達量）+0.244 6×（LINC00461 表達量）-0.424 4×（LINC01132表達量）+0.235 5×（RP1-79C4.4表達量）+0.298 1×（RP4-816N1.7表達量）+0.341 4×（RP3-380B8.4表達量）

5、預后風險模型的評估

通過計算模型C 指數值為0.82（圖6），介于0.71～0.90之間，為中等準確度。構建時間依賴ROC，并計算AUC值。預測3 年、5 年結腸癌患者生存率的AUC 值分別為0.79 和0.84（圖7），0.7

圖7 風險模型時間依賴受試者工作特征曲線

根據模型計算出每個患者的RS 值，以中位數作為截斷值分組，其中高、低風險組各有168 例，繪制風險評分曲線圖、生存狀態散點圖和7 個lncRNA 的表達熱圖（圖8）。從生存狀態散點圖中我們可以發現，高風險組中有更多患者處于死亡狀態，且生存時間較短，并與RS呈正相關，證明模型預測的準確性。表達熱圖中，高風險組的ELFN1-AS1、RP5-884M6.1、LINC00461、RP1-79C4.4、RP4-816N1.7、RP3-380B8.4 表達量較高，而LINC01132 基因的表達量較低，與以上分析結果一致。

圖8 風險評分3圖聯動。A：評分曲線圖；B：生存狀態散點圖；C：lncRNA表達熱圖

對高、低風險組進行生存分析，繪制總生存期的K-M生存曲線（圖9），顯示高風險組患者的總體生存率更低，兩組比較差異有統計學意義（P<0.000 1），且高風險組某一節點存活人數更少、單位時間內的死亡人數更多。以上說明此模型可較好評估結腸癌患者預后。

圖9 風險評分高低組間生存分析圖。A：生存曲線圖；B：隨訪某一時間節點兩組間存活例數；C：隨訪某一時間段死亡例數

6、ceRNA網絡構建及GO、KEGG富集分析

對模型中的5 個lncRNA（ELFN1-AS1、RP5-884M6.1、LINC00461、LINC01132、RP1-79C4.4）構建ceRNA 網絡。通過數據庫lncBOOK 預測miRNA，通過miRDB、TargeScan 和miRTarBase 數據庫預測mRNA，并根據相關性篩選構成lncRNA-miRNA-mRN 的ceRNA 網絡，使用軟件Cytoscape繪圖。

對下調lncRNA（LINC01132）匹配到的154個mRNA及上調lncRNA（ELFN1-AS1、RP5-884M6.1、LINC00461、RP1-79C4.4）匹配到的的351 個mRNA 進行GO 富集分析。下調lncRNA 相關mRNA 涉及肌肉系統進程（muscle system process）、肌肉器官發展（muscle organ development）等生物過程；富集在肌纖維膜（sarcolemma）、小凹（Caveola）等細胞組件；參與整合蛋白綁定（integrin binding）、離子通道綁定（ion channel binding）等分子功能（圖10）。上調lncRNA 相關的mRNA 涉及細胞外基質組織（extracellular matrix organization）、細胞外結構組織（extracellular structure organization）等生物過程；富集在內含膠原蛋白細胞外基質（collagen-containing extracellular matrix）、細胞重要邊緣（cell leading edge）等細胞組件；參與整合蛋白綁定（integrin binding）、細胞外基質結構成分（extracellular matrix structural constituent）等分子功能（圖11）。

圖10 下調長鏈非編碼RNA基因本體論（GO）富集分析圖

圖11 上調長鏈非編碼RNA基因本體論（GO）富集分析圖

再對lncRNA 相關mRNA 進行KEGG 富集分析，分組顯示前15 條通路，發現下調lncRNA 相關基因顯著富集在肌動蛋白細胞骨架的調控（regulation of actin cytoskeleton）、癌癥中的蛋白聚糖（proteoglycans in cancer）、PI3K-Akt 信號通路（PI3K-Akt signaling pathway）等通路（圖12）；上調lncRNA相關基因顯著富集在細胞粘附分子（Cell adhesion molecules）、局灶性粘連（Focal adhesion）、吞噬體（phagosome）等通路（圖13）。

圖12 下調長鏈非編碼RNA京都基因與基因組大百科全書數據庫（KEGG）富集分析圖

圖13 上調長鏈非編碼RNA京都基因與基因組大百科全書數據庫（KEGG）富集分析圖

討論

lncRNA 在細胞增殖分化、血管再生和細胞活力等方面發揮了重要的作用，當其表達發生改變時可影響腫瘤的發生和轉移［16］。隨著生物學研究方法不斷進步，越來越多lncRNA 被深入研究，通過對其進行整合分析，構建了多種癌癥不同的預后風險模型，對腫瘤預后判斷提供了重要幫助［17-18］。大量研究表明，lncRNA 參與了結直腸癌的進展，并可能成為診斷、治療及預后判斷的潛在生物標志物［19］。

研究中，我們使用TCGA 數據庫結腸癌患者的轉錄組數據進行分析并構建預后風險模型。分析方法上，我們在COX回歸分析的基礎上進行Lasso回歸分析，能有效過濾其中的共線性因素，構建了一個更為簡練而準確的模型。并對每個lncRNA 進行單獨的生存分析，保證了模型中每個lncRNA 均對結腸癌患者生存率具有顯著影響。在模型準確性評估方面，通過多種指標進行評估，包括C 指數值、ROC及AUC值、RS值可視化和高、低風險組生存分析，均表明模型具有較好的結腸癌預后預測能力。本研究篩選出7 個lncRNA 構成的結腸癌預后風險模型，與單一的生物標志物相比，多個綜合的生物標志物系統可以提高容錯性和預測的準確性。目前此模型在國內外的研究中未有報道，提示其可為結腸癌的預后判斷提供新的參考，及為結腸癌的分子機制研究提供新的方向。

模型中的7 個lncRNA 均被Ensembl 數據庫所注解，關于它們在腫瘤中的作用已有不同程度的認識。關于下調lncRNA（LINC01132），研究發現其在膠質母細胞瘤維持缺氧誘導因子的高活性和腫瘤細胞遷移率［20］；有研究通過生物分析發現其在卵巢癌中低表達，與患者的生存率明顯相關，可用于評估卵巢癌患者預后，與本研究結果相符［21］。對于ELFN1-AS1則有較多的報道，如Lei等［22］發 現ELFN1-AS1 通過調節MIR-4644/TRIM44 軸加速了結直腸癌的增殖和遷移。LINC00461 是一個明星lncRNA，研究發現LINC00461 在多種癌癥中高表達，通過不同的途徑影響著腫瘤的發生、進展，其中包括結直腸癌［23-24］。Fu等［25］通過TCGA 數據分析發現RP4-816N1.7 與結腸腺癌預后密切相關。模型中的4個lncRNA（LINC01132、ELFN1-AS1、LINC00461、RP4-816N1.7）均有相關的研究報道，大多與結腸癌有關，從側面驗證本研究預后模型的準確性。然而根據PubMed數據庫的搜索，沒有發現模型中其余3個lncRNA（RP5-884M6.1、RP3-380B8.4、RP1-79C4.4）研究報道，其在腫瘤，特別是結腸癌中的作用機制有待發掘。

為了探究模型中的lncRNA 潛在的生物學功能及如何影響結腸癌進程，本研究構建ceRNA 網絡，通過GO 分析發現，lncRNA 相關mRNA 基因主要富集在細胞膜、細胞質等細胞組件，參與整合蛋白綁定、離子通道綁定、細胞外基質結構成分等分子功能，涉及肌肉系統進程、肌肉器官發展、細胞外基質組織、細胞外結構組織等生物過程。KEGG 富集分析提示下調lncRNA 可能是通過肌動蛋白細胞骨架的調控、癌癥中蛋白聚糖、PI3K-Akt 信號通路等抑制結腸癌進展，上調lncRNA 可能是通過細胞粘附分子、局灶性粘連、吞噬體等通路促進結腸癌進展。

雖然我們構建了一個有效的結腸癌預后風險模型，但也存在著一些不足。首先，從TCGA 中提取的結腸癌樣本量不夠大；其次，模型的預后價值評估不是十分令人滿意；再者，lncRNA、miRNA 和mRNA 間的相關程度及lncRNA 如何影響結腸癌進程主要通過在線數據庫進行預測，需要進一步的實驗研究。

綜上所述，本研究對TCGA 中結腸癌轉錄組數據進行整理、分析，篩選出7 個預后相關lncRNA，構建預后風險模型，通過評估證明其具有較好預測患者生存預后的準確性，同時每個lncRNA 是潛在單獨的預后生物標志物，對結腸癌患者臨床預后評估具有一定的參考價值。通過GO、KEGG富集分析對模型中lncRNA 影響結腸癌機制進行初步探索。在臨床上，結合結腸癌患者的lncRNA 分子水平的風險評分，可篩選出高風險患者，制定個體化的治療方案，有利于患者病情的預后評估和綜合診治。