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

何天 曹天生

廣州市花都區(qū)人民醫(yī)院胃腸外科，廣州 510800

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

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

材料與方法

1、原始數(shù)據(jù)下載及數(shù)據(jù)預處理

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

2、差異表達基因分析

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

3、挑選預后相關lncRNA

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

4、構建預后風險模型

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

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

5、評估預后風險模型

先計算模型的C 指數(shù)值，其可用來評價模型的預測能力，主要用于計算COX 回歸模型中預測值與真實之間的區(qū)分度［12］。一般認為，如C指數(shù)介于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 指數(shù)意義相似，一般認為0.85

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

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

圖1 結(jié)腸癌TCGA數(shù)據(jù)提取分析流程圖

結(jié)果

1、TCGA中結(jié)腸癌患者臨床病理信息

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

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

2、配對樣本差異分析結(jié)果

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

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

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

圖4 長鏈非編碼RNA火山圖

3、結(jié)腸癌預后關鍵lncRNA的篩選

構建癌組織DElncRNA 表達矩陣，先進行單因素COX回歸分析，篩選出40 個lncRNA。隨后進行Lasso 回歸模型分析，得到34 個候選lncRNA。以lncRNA 表達值的中位數(shù)作為分界值，進行高低分組生存分析，發(fā)現(xiàn)其中的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%可信區(qū)間，繪制風險森林圖（圖6）。其中1 個lncRNA（LINC01132）的HR<1，提示其表達量增加可降低結(jié)腸癌不良預后的風險，是患者的保護因素；6 個 lncRNA（ELFN1-AS1、RP5-884M6.1、LINC00461、RP1-79C4.4、RP4-816N1.7、RP3-380B8.4）的HR>l，提示它們表達量增加可增加結(jié)腸癌不良預后的風險，是患者的危險因素。

4、預后風險模型的構建

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

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 指數(shù)值為0.82（圖6），介于0.71～0.90之間，為中等準確度。構建時間依賴ROC，并計算AUC值。預測3 年、5 年結(jié)腸癌患者生存率的AUC 值分別為0.79 和0.84（圖7），0.7

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

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

圖8 風險評分3圖聯(lián)動。A：評分曲線圖；B：生存狀態(tài)散點圖；C：lncRNA表達熱圖

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

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

6、ceRNA網(wǎng)絡構建及GO、KEGG富集分析

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

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

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

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

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

圖12 下調(diào)長鏈非編碼RNA京都基因與基因組大百科全書數(shù)據(jù)庫（KEGG）富集分析圖

圖13 上調(diào)長鏈非編碼RNA京都基因與基因組大百科全書數(shù)據(jù)庫（KEGG）富集分析圖

討論

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

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

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

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

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

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