引用本文： 倪耀軍, 徐忠能, 陳勝, 王俊. 非小細胞肺癌細胞 A549 順鉑耐藥潛在機制的全轉錄組測序分析. 中國胸心血管外科臨床雜志, 2021, 28(1): 35-42. doi: 10.7507/1007-4848.202003132 復制
肺癌是最常見的惡性腫瘤之一，發病率與死亡率位居惡性腫瘤之首，并呈逐年上升趨勢[1-2]。其中，非小細胞肺癌（non-small cell lung cancer，NSCLC）約占肺癌總數的 85%。NSCLC 易出現血行播散和區域淋巴結轉移，患者預后差，5 年生存率不足 15%。近年來隨著分子靶向藥物、抗血管藥物以及免疫靶向藥物的臨床應用，NSCLC 的治療取得了很大進展，尤其是分子靶向藥物的應用顯著提高了 NSCLC 患者的 5 年生存率[3-4]。但研究[5-7]發現，我國 NSCLC 患者表皮生長因子受體（EGFR）突變率約為 20%～30%，存在棘皮動物微管相關類蛋白4/間變性淋巴瘤激酶（EML4-ALK）融合基因的僅占 4%，因此能夠接受分子靶向治療的 NSCLC 患者很少。對于大多數晚期或術后復發的 NSCLC 患者，基于鉑類藥物的化療仍是一線治療方案，其中以順鉑的治療效果最佳[8-9]。然而，固有性耐藥及獲得性耐藥嚴重限制了以順鉑為基礎的聯合化療療效。研究[10]顯示 NSCLC 患者順鉑耐藥率高達 63%，耐藥的發生嚴重影響順鉑的療效，儼然已經成為 NSCLC 患者治療及管理所面臨的主要臨床挑戰。研究[11-12]證實，腫瘤通過多種方法調控自身對順鉑的敏感性，包括順鉑在細胞中的累積及排出、DNA 損傷修復、藥物的失活作用和凋亡抑制等，雖然許多基礎實驗已經闡明腫瘤順鉑耐藥的部分機制，但是目前并沒有很好的辦法替代順鉑的抗腫瘤地位或改變耐藥腫瘤對順鉑的敏感性。因此，積極尋找順鉑耐藥的關鍵分子機制尤為重要。
全轉錄組是指特定物種或細胞在某一狀態下轉錄產生的所有 RNA 總和，包括 mRNA 和非編碼 RNA。非編碼 RNA 包括 microRNA（miRNA），長鏈非編碼 RNA（long non-coding RNA，lncRNA）和環狀 RNA（circle RNA，circRNA）。本研究進行全轉錄組測序和分析，比較野生型人肺腺癌細胞株 A549 和 A549/DDP 耐藥細胞株中 lncRNA、miRNA 和 mRNA 的表達譜。采用整合分析方法以及 ceRNA 調控網絡分析，識別與 DDP 耐藥相關的信號通路和關鍵基因。
1 材料與方法
1.1 實驗材料
人肺腺癌細胞株 A549 購于中國科學院上海細胞庫；DDP 購自 Sigma 公司；1640 細胞培養基和胎牛血清（FBS）購自于 Gibco 公司；總 RNA 提取試劑盒 Trizol 購于 Invitrogen 公司。肺腺癌細胞株 A549 在含有 10% FBS 的 1640 培養基中，置于 37℃、CO2 體積分數為 5%、相對濕度為 90%的培養箱中進行常規培養。通過間歇誘導的方法建立 DDP 耐藥細胞株（A549/DDP）[13]，即在 A549 細胞培養基中遞增 DDP 濃度直至耐藥，建成耐藥細胞系后，在無 DDP 培養基中繼續培養 20 周后，該細胞株的耐藥性仍穩定，本研究中的維持耐藥濃度為 1 000 ng/mL。
1.2 方法
1.2.1 全轉錄組測序
產生 DDP 耐藥性的 A549/DDP 細胞（實驗組，3 個樣本）和 DPP 敏感的 A549 細胞（對照組，3 個樣本）用 Trizol 法提取 RNA，全轉錄組測序（miRNA-seq、circRNA-seq 和 lncRNA-seq）由上海生工公司完成，測序平臺為 Illumina HiseqTM。
1.2.2 測序數據質量控制
用 Trimmomatic（v 3.6）[14]對測序數據進行質量控制。去掉序列（reads）兩端連續質量低于 10 的堿基；去掉少于 80% 的堿基質量大于 20 的 reads；過濾掉長度小于 35 bp 的 reads，對于 miRNA-seq 數據，過濾掉長度小于 16 bp 的 reads。從而得到過濾后的可用序列（clean reads）。
1.2.3 lncRNA-seq 數據差異表達基因以及富集分析
首先將 clean reads 定位到人類的參考基因組（Hg38，ENSEMBL）上[15]，進行基因組比對。對于已知基因，根據基因組數注釋數據，將其分為 lncRNA 和蛋白質編碼基因（protein coding gene），統稱為 mRNA。logFC>1 并且 q<0.01（矯正過的 P 值）作為篩選實驗組與對照組之間差異基因（differently expressed genes，DEG）的閾值。用 clusterProfiler 對 DEG 進行 GO 功能富集分析和 KEGG 通路富集分析[16]。超幾何檢驗顯著性閾值 P<0.05 視為顯著富集結果。GO 功能分為三個大類：分子功能（molecular function，MF），細胞組分（cellular component，CC），生物過程（biological process，BP）。
1.2.4 circRNA-seq 數據表達分析以及功能富集
將 clean reads 通過比對到人類參考基因組（ENSEMBL）上以后，通過 CIRCexplorer2[17]來鑒定 circRNA，并對篩選獲得的 circRNA candidate 進行注釋。使用反向剪接序列（back-spliced reads）的數目來估計 circRNA 的表達量，基于 circRNA 的表達量通過 edgeR[18]軟件包篩選兩組樣本中實驗組與對照組之間的顯著差異表達 circRNA，設定閾值為|logFC|>1，并且 P<0.05。另外，circRNA 在至少 3 個樣本中表達值大于 0，其余 circRNA 認為低表達豐度 circRNA 而不進入下游差異分析。對差異表達 circRNA 的主基因（host gene）使用 clusterProfiler 分析其富集的 GO 功能和 KEGG 通路，超幾何檢驗顯著性閾值 P<0.05 視為顯著富集結果。
1.2.5 miRNA-seq 數據分析
利用 bowtie 軟件（v0.12.9）[19]將 clean reads map 到人類參考基因組（ENSEMBL），統計 reads 數目和百分比，過濾不能比對到參考基因組的 reads。利用 HTSeq 工具（V 0.9.1）獲得每個樣本中成熟 miRNA 比對上的 reads 數并通過 CPM（count per million）的方法對 miRNA 表達豐度進行定量。利用 edgeR 包 quasi-likelihood F-tests 方法篩選實驗組與對照組之間的顯著差異表達 miRNA，設定閾值為|logFC|>1，FDR<0.01。
1.2.6 全轉錄組數據聯合分析
1.2.6.1 miRNA 靶點預測
通過 miRanda[17]工具預測差異表達的 miRNA 與差異表達 circRNA，lncRNA，mRNA 的 3’UTR 區域的結合靶點，使用參數-sc 140-en-20-strict。
1.2.6.2 miRNA 功能預測
通過預測 miRNA 與 mRNA 的調控關系及 miRNA 和 mRNA 的表達變化情況，本研究選擇 mRNA 與 miRNA 有結合靶點并且表達變化與 miRNA 相反的 mRNA 作為 miRNA 最終的靶基因，得到 miRNA-mRNA 的靶向關系對。然后根據 mRNA 間接預測 miRNA 的功能，利用 R 的 clusterProfiler 包分別對表達上調 miRNA 靶基因和表達下調 miRNA 靶基因以及所有差異 miRNA 的靶基因（mRNA）的 GO 功能富集［biological process（BP）］和 KEGG 通路分析，篩選 P<0.05 的結果作為顯著富集結果。
1.2.6.3 ceRNA 網絡構建
首先得到 miRNA-mRNA 的靶向關系對，分別結合 miRNA-lncRNA 和 miRNA-circRNA 靶向關系對，以及 miRNA、lncRNA、mRNA、circRNA 在樣本中表達譜，通過 miRsponge 預測與 mRNA 有關的 lncRNA 和 circRNA，組成 ceRNA 關系對[20]。預測得到的 ceRNA 關系對及結合 miRNA 網絡均通過 cytoscape[21]工具構建關系網絡。
2 結果
2.1 lncRNA-seq 數據分析結果
根據設定的閾值|log2FC|>1 & q<0.01，總共篩選到了 4 517 個差異表達轉錄本（mRNA）。取實驗組和對照組樣本的差異 mRNA 的交集，與對照組相比，其中 882 個 mRNA 是下調的，554 個 mRNA 是上調的；219 個 lncRNA 是下調的，273 個 lncRNA 是上調的。差異 mRNA 和差異 lncRNA 的火山圖（圖1 a、1 b）和聚類分析熱圖（圖1 c、1 d）顯示差異 mRNA 和差異 lncRNA 能夠把細胞樣本明顯分為實驗組和對照組，說明我們篩選到的差異 mRNA 具有明顯的差異特性。
圖1
A549/DDP 細胞和 A549 細胞差異 mRNA 和差異 lncRNA 的火山圖（a、b）和聚類分析熱圖（c、d）
紅色代表上調，藍色代表下調；c 和 d 中橫坐標代表樣本，縱坐標代表差異 mRNA
對差異變化顯著的 mRNA 進行了 KEGG 通路和 GO 功能富集分析，差異 mRNA 主要富集在 ko05202（transcriptional misregulation in cancers）、ko05222（small cell lung cancer）、ko04668［腫瘤壞死因子（TNF）信號通路］和 ko04210（apoptosis）等信號通路（圖2）。
圖2
A549/DDP 細胞和 A549 細胞差異 mRNA KEGG 富集（a）和 GO（b）富集氣泡圖
縱軸表示通路或 GO 功能注釋信息，橫軸表示通路對應的富集因子；
2.2 circRNA-seq 數據分析結果
6 個樣本中共鑒定 9 324 個 circRNA，對應3 681 個主基因（host gene）。共得到 123 個差異表達 circRNA，其中表達上調 17 個，表達下調 106 個。差異 circRNA 的分布熱圖（圖3a ）顯示差異 circRNA 顯著把 A549/DDP 細胞和 A549 細胞區分開。對差異變化顯著的 circRNA-hosting gene 進行了 GO 功能富集分析和 KEGG 通路富集分析，hsa04115（p53 信號通路）為最顯著的信號通路（P=0.006），多種生物過程的陽性調節（positive regulation of multi-organism progress），溶酶體膜和核內小分子 GTP 結合蛋白（lysosomal membrane and Ran GTPase binding）等為比較顯著的 GO 功能條目（圖3b、3c）。
圖3
A549/DDP 細胞和 A549 細胞差異 circRNA 的聚類圖（a）和功能富集圖（b～e）
a：差異 circRNA 二維聚類分布熱點圖，顏色有綠色-黑色-紅色表達 circRNA 表達豐度逐漸增高；b：差異 circRNA 富集的主要 KEGG 信號通路；c～e：差異 circRNA 富集的主要 GO（BP，CC 和 MF）功能條目；b～e：縱坐標代表信號通路信息和 GO 功能描述信息，橫坐標富集到該 Term 的差異基因數，
2.3 miRNA-seq 測序數據分析
對兩組樣本中 miRNA 進行差異表達分析，共得到 145 個差異表達的 miRNA，其中包括 29 個上調的 miRNA 和 126 個下調的 miRNA。聚類分析結果顯示這些差異 miRNA 可以把樣本分為具有顯著特點的耐藥組和對照組（圖4a）。通過靶點預測以及相反的表達變化，預測 miRNA（up）-mRNA（down）和 miRNA（down）-mRNA（up）的調控關系（結果未展示）。接著 miRNA 的靶基因進行功能富集分析，結果顯示差異 miRNA 的靶基因主要富集在心肌病以及癌癥相關的信號通路上（圖4b、4c），如：肥厚性心肌病（hypertrophic cardiomyopathy，HCM），擴張型心肌病（dilated cardiomyopathy，DCM），心律失常性右心室心肌病（arrhythmogenic right ventricular cardiomyopathy，ARVC），胰腺癌（pancreatic cancer）和 p53 信號通路等，其中 p53 信號通路也富集了差異 circRNA-hosting gene。另外，差異 miRNA 靶基因主要富集在組織細胞外基質（extracellular matrix organization），細胞外結構組織（extracellular structure organization）和上皮細胞增殖（epithelial cell proliferation）等功能上。
圖4
A549/DDP 細胞和 A549 細胞差異 miRNA 的聚類圖（a）和功能富集圖（b、c）
a：差異 miRNA 聚類圖，紅色代表上調，綠色代表下調；b：差異 miRNA 顯著富集的 KEGG 信號通路；c：差異 miRNA 功能富集的顯著 GO 功能條目；b、c：縱坐標代表信號通路信息和 GO 功能描述信息，橫坐標富集到該 Term 的差異基因數，
2.4 ceRNA 網絡分析
首先對 miRNA 與 lncRNA，mRNA，circRNA 靶點預測。根據 ceRNA 機制，結合已得到的 miRNA-lncRNA/circRNA 及 miRNA（up）-mRNA（down）和 miRNA（down）-mRNA（up），分別預測 lncRNA-mRNA 和 circRNA-mRNA 的 ceRNA 關系對，并構建 miRNA-circRNA-mRNA 和 miRNA-lncRNA-mRNA 等 ceRNA 網絡（圖5 ）。根據拓撲學性質分析，hsa-novel-32-mature（up）和 hsa-miR-365a-5p（up）的 degree 分別為 36 和 25（圖5b），hsa-miR-665（down）和 hsa-miR-370-3p（down）的 degree 分別為 99 和 79（圖5c ），高于本網絡中其它成員的 degree，推測這些 miRNA 在 A549 順鉑耐藥性產生機制中作用顯著。
圖5
ceRNA 三元網絡圖
a：miRNA（down）-circRNA-mRNA（up），三元網絡；b：miRNA（up）-lncRNA-mRNA（down）三元網絡；c：miRNA（down）-lncRNA-mRNA（up）三元網絡；紅色節點表示表達上調，綠色節點表示表達下調，三角形節點代表 miRNA，V 形節點代表 circRNA，圓形節點代表 mRNA，菱形節點代表 lncRNA
對得到 miRNA-lncRNA-mRNA 和 miRNA-circRNA-mRNA 網絡進行整合，僅保留兩者都存在 mRNA 的關系對，構建 miRNA-circRNA-lncRNA-mRNA 四元網絡（圖6）。此網絡包括了 12 個 miRNA，4 個 circRNA，23 個 lncRNA 和 9 個 mRNA 節點。對miRNA-circRNA-lncRNA-mRNA 四元網絡進行拓撲學分析，得到 degree 排名在前 10 名的分子（表1 ）。Degree 越大表示該分子在此網絡中和其它成員相互作用越多，作用越關鍵。
圖6
miRNA-circRNA-lncRNA-mRNA 四元網絡
紅色節點表示表達上調，綠色節點表示表達下調，三角形節點代表 miRNA，V 形節點為 circRNA，菱形節點代表 lncRNA，圓形節點代表 mRNA；藍色連接線為 miRNA-lncRNA 調控關系，橙色連接線為 miRNA-mRNA 關系對，黑色連接線為 lncRNA-mRNA 關系對，粉色連接線為 miRNA-circRNA 調控關系，青色連接線為 circRNA-mRNA 關系對
表1
miRNA-circRNA-lncRNA-mRNA 四元網絡節點信息（degree top10）
3 討論
順鉑是一種重要的肺癌化療藥物，但是固有性及獲得性耐藥嚴重限制其化療效果。據報道，耐藥機制涉及許多細胞系多種基因過表達和/或低表達，從而引起一系列病理變化[22-24]。目前對基因調控在 NSCLC 順鉑耐藥中的認識仍然有限。在本研究中，我們評估了 lncRNA、miRNA 和 cirRNA 在順鉑耐藥細胞系 A549/CDDP 中的表達譜。與野生型 A549 細胞系相比，微陣列分析顯示了一組差異表達的非編碼 RNA，其中 4517 個 lncRNA 和 123 個 cirRNA 以及 145 個 miRNA 在 A549/CDDP 細胞中有差異表達。
基于差異表達的 mRNA，我們進行了 KEGG 和 GO 功能富集分析，以闡明哪些特定的信號通路可能參與了控制 A549/DDP 和 A549 的細胞特征。KEGG 結果表明，差異 mRNA 顯著富集在和癌癥相關的通路上，如：ko05202（transcriptional misregulation in cancers），ko05222（small cell lung cancer），ko04668（TNF 信號通路）和 hsa04115（p53 信號通路）等。TNF 信號通路為 TNF 信號通路，TNF 是一種可以直接殺傷腫瘤細胞，而對正常細胞沒有明顯細胞毒性的細胞因子，此信號通路參與細胞的生存、死亡和分化。肺癌晚期多會發現 TNF 受體表達缺失，人 A549 肺腺癌細胞對 TNF-α 和順鉑的細胞毒性作用具有抗性[22-24]。本研究表明 TNF 信號通路富集了差異 mRNA，我們推測 A549/DDP 耐藥機制可能與 TNF 信號通路異常調控有關。
p53 基因是一種抑癌基因，編碼的蛋白質是一種轉錄因子，其控制著細胞周期的啟動，在絕大多數腫瘤細胞會出現該基因的突變[25]。眾多報道顯示 p53 信號通路與多腫瘤產生的的順鉑耐藥相關。Wang 等[26]根據分子成像研究表明順鉑通過阻斷細胞周期發揮細胞毒作用，可能部分受 p53 信號通路調控。Mg2+依賴性蛋白磷酸酶 1δ（PPM1D）通過減弱細胞周期檢測點激酶 1（Chk1）和 p53 的激活，使卵巢癌細胞對順鉑產生耐藥性[27]。在維生素 C（VC）輔助治療腸癌的研究中發現 VC 可以通過上調 p53 來增強順鉑的敏感性和凋亡[28]。通過分析 A549 與 A549/CDDP 細胞系的全轉錄組測序結果發現差異 mRNA 顯著富集在 p53 信號通路上，提示我們 p53 信號通路可能參與了 A549/CDDP 耐藥機制。
對 lncRNA，miRNA 和 circRNA 表達譜進行綜合分析得到 ceRNA 調控網絡，其中 hsa-miR-125a-5p 和 hsa-miR-125b-5p 在 A549/DPP 抗性細胞中表達下調，是 miRNA-circRNA-lncRNA-mRNA 四元網絡圖中的關鍵 miRNA。成熟的 hsa-miR-125a-5p 來自于 pre-miR-125a 的 5'端，在 NSCLC 癌組織中 hsa-miR-125a-5p 的表達低于鄰近正常肺組織，Spearman 相關檢驗結果表明 hsa-miR-125a-5p 的表達與病理分期或者淋巴結轉移呈負相關[29]。Liu 等[30]研究發現喉癌組織和 Hep-2 喉癌干細胞（Hep-2-cscs）中，miRNA-125a 表達降低，在體內外，miRNA-125a 的功能增強顯著提高了 Hep-2-CSCs 對順鉑的敏感性。Jiang 等[31]發現不同肺癌細胞系中 hsa-miR-125a-5p 的表達均低于人支氣管上皮細胞，并且抑制了 A549 細胞的增殖以及誘導了細胞凋亡，通過分別阻斷野生型的 p53 和 hsa-miR-125a-5p，發現 hsa-miR-125a-5p 通過 p53 依賴通路誘導人肺癌細胞凋亡。另外 Zhang 等[32]也證明 miRNA 125a 可以通過靶向 p53 的 3'‐UTR 翻譯和轉錄抑制 p53 基因。hsa-miR-125b-5p 與 hsa-miR-125a-5p 高度同源，被報道在多種癌細胞中作用顯著，比如肺癌與乳腺癌細胞[33-34]。與卵巢癌細胞系與 OV2008 細胞相比，miR-125b 在耐順鉑細胞株中（C13*）的表達水平升高，上調 microRNA-125b 可顯著抑制順鉑誘導的細胞毒性和凋亡，并可增加 OV2008 和 C13*細胞對順鉑的耐藥性[35]。而我們的研究表明 A549/DPP 抗性細胞中 micRNA-125b 表達是下調的，與上述報道有一定差異，可能是由于不同的細胞種類或者不同的處理方式引起的表達差異。基于前期研究和我們的分析結果推測 hsa-miR-125a-5p 和 hsa-miR-125b-5p 參與了 A549/DPP 抗性機制，可能是逆轉順鉑抗性的潛在靶點。
本研究通過全轉錄組測序比較野生型 A549 細胞和 A549/DPP 順鉑耐藥細胞表達譜的差別，初步篩選了與 NSCLC 順鉑耐藥相關基因，獲得具有潛在意義的耐藥基因，并通過生物信息學分析闡釋其可能機制。但是，我們在研究中獲得的差異表達及機制層面的分析結果尚需進行基礎實驗及臨床大樣本驗證，這也是我們目前正在進行的工作，其深入機制仍有待進一步探索。
利益沖突：無。
作者貢獻：王俊負責實驗設計及部分數據分析；倪耀軍負責論文撰寫及實驗的具體實施；徐忠能和陳勝負責生物信息分析解讀及部分數據分析。
肺癌是最常見的惡性腫瘤之一，發病率與死亡率位居惡性腫瘤之首，并呈逐年上升趨勢[1-2]。其中，非小細胞肺癌（non-small cell lung cancer，NSCLC）約占肺癌總數的 85%。NSCLC 易出現血行播散和區域淋巴結轉移，患者預后差，5 年生存率不足 15%。近年來隨著分子靶向藥物、抗血管藥物以及免疫靶向藥物的臨床應用，NSCLC 的治療取得了很大進展，尤其是分子靶向藥物的應用顯著提高了 NSCLC 患者的 5 年生存率[3-4]。但研究[5-7]發現，我國 NSCLC 患者表皮生長因子受體（EGFR）突變率約為 20%～30%，存在棘皮動物微管相關類蛋白4/間變性淋巴瘤激酶（EML4-ALK）融合基因的僅占 4%，因此能夠接受分子靶向治療的 NSCLC 患者很少。對于大多數晚期或術后復發的 NSCLC 患者，基于鉑類藥物的化療仍是一線治療方案，其中以順鉑的治療效果最佳[8-9]。然而，固有性耐藥及獲得性耐藥嚴重限制了以順鉑為基礎的聯合化療療效。研究[10]顯示 NSCLC 患者順鉑耐藥率高達 63%，耐藥的發生嚴重影響順鉑的療效，儼然已經成為 NSCLC 患者治療及管理所面臨的主要臨床挑戰。研究[11-12]證實，腫瘤通過多種方法調控自身對順鉑的敏感性，包括順鉑在細胞中的累積及排出、DNA 損傷修復、藥物的失活作用和凋亡抑制等，雖然許多基礎實驗已經闡明腫瘤順鉑耐藥的部分機制，但是目前并沒有很好的辦法替代順鉑的抗腫瘤地位或改變耐藥腫瘤對順鉑的敏感性。因此，積極尋找順鉑耐藥的關鍵分子機制尤為重要。
全轉錄組是指特定物種或細胞在某一狀態下轉錄產生的所有 RNA 總和，包括 mRNA 和非編碼 RNA。非編碼 RNA 包括 microRNA（miRNA），長鏈非編碼 RNA（long non-coding RNA，lncRNA）和環狀 RNA（circle RNA，circRNA）。本研究進行全轉錄組測序和分析，比較野生型人肺腺癌細胞株 A549 和 A549/DDP 耐藥細胞株中 lncRNA、miRNA 和 mRNA 的表達譜。采用整合分析方法以及 ceRNA 調控網絡分析，識別與 DDP 耐藥相關的信號通路和關鍵基因。
1 材料與方法
1.1 實驗材料
人肺腺癌細胞株 A549 購于中國科學院上海細胞庫；DDP 購自 Sigma 公司；1640 細胞培養基和胎牛血清（FBS）購自于 Gibco 公司；總 RNA 提取試劑盒 Trizol 購于 Invitrogen 公司。肺腺癌細胞株 A549 在含有 10% FBS 的 1640 培養基中，置于 37℃、CO2 體積分數為 5%、相對濕度為 90%的培養箱中進行常規培養。通過間歇誘導的方法建立 DDP 耐藥細胞株（A549/DDP）[13]，即在 A549 細胞培養基中遞增 DDP 濃度直至耐藥，建成耐藥細胞系后，在無 DDP 培養基中繼續培養 20 周后，該細胞株的耐藥性仍穩定，本研究中的維持耐藥濃度為 1 000 ng/mL。
1.2 方法
1.2.1 全轉錄組測序
產生 DDP 耐藥性的 A549/DDP 細胞（實驗組，3 個樣本）和 DPP 敏感的 A549 細胞（對照組，3 個樣本）用 Trizol 法提取 RNA，全轉錄組測序（miRNA-seq、circRNA-seq 和 lncRNA-seq）由上海生工公司完成，測序平臺為 Illumina HiseqTM。
1.2.2 測序數據質量控制
用 Trimmomatic（v 3.6）[14]對測序數據進行質量控制。去掉序列（reads）兩端連續質量低于 10 的堿基；去掉少于 80% 的堿基質量大于 20 的 reads；過濾掉長度小于 35 bp 的 reads，對于 miRNA-seq 數據，過濾掉長度小于 16 bp 的 reads。從而得到過濾后的可用序列（clean reads）。
1.2.3 lncRNA-seq 數據差異表達基因以及富集分析
首先將 clean reads 定位到人類的參考基因組（Hg38，ENSEMBL）上[15]，進行基因組比對。對于已知基因，根據基因組數注釋數據，將其分為 lncRNA 和蛋白質編碼基因（protein coding gene），統稱為 mRNA。logFC>1 并且 q<0.01（矯正過的 P 值）作為篩選實驗組與對照組之間差異基因（differently expressed genes，DEG）的閾值。用 clusterProfiler 對 DEG 進行 GO 功能富集分析和 KEGG 通路富集分析[16]。超幾何檢驗顯著性閾值 P<0.05 視為顯著富集結果。GO 功能分為三個大類：分子功能（molecular function，MF），細胞組分（cellular component，CC），生物過程（biological process，BP）。
1.2.4 circRNA-seq 數據表達分析以及功能富集
將 clean reads 通過比對到人類參考基因組（ENSEMBL）上以后，通過 CIRCexplorer2[17]來鑒定 circRNA，并對篩選獲得的 circRNA candidate 進行注釋。使用反向剪接序列（back-spliced reads）的數目來估計 circRNA 的表達量，基于 circRNA 的表達量通過 edgeR[18]軟件包篩選兩組樣本中實驗組與對照組之間的顯著差異表達 circRNA，設定閾值為|logFC|>1，并且 P<0.05。另外，circRNA 在至少 3 個樣本中表達值大于 0，其余 circRNA 認為低表達豐度 circRNA 而不進入下游差異分析。對差異表達 circRNA 的主基因（host gene）使用 clusterProfiler 分析其富集的 GO 功能和 KEGG 通路，超幾何檢驗顯著性閾值 P<0.05 視為顯著富集結果。
1.2.5 miRNA-seq 數據分析
利用 bowtie 軟件（v0.12.9）[19]將 clean reads map 到人類參考基因組（ENSEMBL），統計 reads 數目和百分比，過濾不能比對到參考基因組的 reads。利用 HTSeq 工具（V 0.9.1）獲得每個樣本中成熟 miRNA 比對上的 reads 數并通過 CPM（count per million）的方法對 miRNA 表達豐度進行定量。利用 edgeR 包 quasi-likelihood F-tests 方法篩選實驗組與對照組之間的顯著差異表達 miRNA，設定閾值為|logFC|>1，FDR<0.01。
1.2.6 全轉錄組數據聯合分析
1.2.6.1 miRNA 靶點預測
通過 miRanda[17]工具預測差異表達的 miRNA 與差異表達 circRNA，lncRNA，mRNA 的 3’UTR 區域的結合靶點，使用參數-sc 140-en-20-strict。
1.2.6.2 miRNA 功能預測
通過預測 miRNA 與 mRNA 的調控關系及 miRNA 和 mRNA 的表達變化情況，本研究選擇 mRNA 與 miRNA 有結合靶點并且表達變化與 miRNA 相反的 mRNA 作為 miRNA 最終的靶基因，得到 miRNA-mRNA 的靶向關系對。然后根據 mRNA 間接預測 miRNA 的功能，利用 R 的 clusterProfiler 包分別對表達上調 miRNA 靶基因和表達下調 miRNA 靶基因以及所有差異 miRNA 的靶基因（mRNA）的 GO 功能富集［biological process（BP）］和 KEGG 通路分析，篩選 P<0.05 的結果作為顯著富集結果。
1.2.6.3 ceRNA 網絡構建
首先得到 miRNA-mRNA 的靶向關系對，分別結合 miRNA-lncRNA 和 miRNA-circRNA 靶向關系對，以及 miRNA、lncRNA、mRNA、circRNA 在樣本中表達譜，通過 miRsponge 預測與 mRNA 有關的 lncRNA 和 circRNA，組成 ceRNA 關系對[20]。預測得到的 ceRNA 關系對及結合 miRNA 網絡均通過 cytoscape[21]工具構建關系網絡。
2 結果
2.1 lncRNA-seq 數據分析結果
根據設定的閾值|log2FC|>1 & q<0.01，總共篩選到了 4 517 個差異表達轉錄本（mRNA）。取實驗組和對照組樣本的差異 mRNA 的交集，與對照組相比，其中 882 個 mRNA 是下調的，554 個 mRNA 是上調的；219 個 lncRNA 是下調的，273 個 lncRNA 是上調的。差異 mRNA 和差異 lncRNA 的火山圖（圖1 a、1 b）和聚類分析熱圖（圖1 c、1 d）顯示差異 mRNA 和差異 lncRNA 能夠把細胞樣本明顯分為實驗組和對照組，說明我們篩選到的差異 mRNA 具有明顯的差異特性。
圖1
A549/DDP 細胞和 A549 細胞差異 mRNA 和差異 lncRNA 的火山圖（a、b）和聚類分析熱圖（c、d）
紅色代表上調，藍色代表下調；c 和 d 中橫坐標代表樣本，縱坐標代表差異 mRNA
對差異變化顯著的 mRNA 進行了 KEGG 通路和 GO 功能富集分析，差異 mRNA 主要富集在 ko05202（transcriptional misregulation in cancers）、ko05222（small cell lung cancer）、ko04668［腫瘤壞死因子（TNF）信號通路］和 ko04210（apoptosis）等信號通路（圖2）。
圖2
A549/DDP 細胞和 A549 細胞差異 mRNA KEGG 富集（a）和 GO（b）富集氣泡圖
縱軸表示通路或 GO 功能注釋信息，橫軸表示通路對應的富集因子；
2.2 circRNA-seq 數據分析結果
6 個樣本中共鑒定 9 324 個 circRNA，對應3 681 個主基因（host gene）。共得到 123 個差異表達 circRNA，其中表達上調 17 個，表達下調 106 個。差異 circRNA 的分布熱圖（圖3a ）顯示差異 circRNA 顯著把 A549/DDP 細胞和 A549 細胞區分開。對差異變化顯著的 circRNA-hosting gene 進行了 GO 功能富集分析和 KEGG 通路富集分析，hsa04115（p53 信號通路）為最顯著的信號通路（P=0.006），多種生物過程的陽性調節（positive regulation of multi-organism progress），溶酶體膜和核內小分子 GTP 結合蛋白（lysosomal membrane and Ran GTPase binding）等為比較顯著的 GO 功能條目（圖3b、3c）。
圖3
A549/DDP 細胞和 A549 細胞差異 circRNA 的聚類圖（a）和功能富集圖（b～e）
a：差異 circRNA 二維聚類分布熱點圖，顏色有綠色-黑色-紅色表達 circRNA 表達豐度逐漸增高；b：差異 circRNA 富集的主要 KEGG 信號通路；c～e：差異 circRNA 富集的主要 GO（BP，CC 和 MF）功能條目；b～e：縱坐標代表信號通路信息和 GO 功能描述信息，橫坐標富集到該 Term 的差異基因數，
2.3 miRNA-seq 測序數據分析
對兩組樣本中 miRNA 進行差異表達分析，共得到 145 個差異表達的 miRNA，其中包括 29 個上調的 miRNA 和 126 個下調的 miRNA。聚類分析結果顯示這些差異 miRNA 可以把樣本分為具有顯著特點的耐藥組和對照組（圖4a）。通過靶點預測以及相反的表達變化，預測 miRNA（up）-mRNA（down）和 miRNA（down）-mRNA（up）的調控關系（結果未展示）。接著 miRNA 的靶基因進行功能富集分析，結果顯示差異 miRNA 的靶基因主要富集在心肌病以及癌癥相關的信號通路上（圖4b、4c），如：肥厚性心肌病（hypertrophic cardiomyopathy，HCM），擴張型心肌病（dilated cardiomyopathy，DCM），心律失常性右心室心肌病（arrhythmogenic right ventricular cardiomyopathy，ARVC），胰腺癌（pancreatic cancer）和 p53 信號通路等，其中 p53 信號通路也富集了差異 circRNA-hosting gene。另外，差異 miRNA 靶基因主要富集在組織細胞外基質（extracellular matrix organization），細胞外結構組織（extracellular structure organization）和上皮細胞增殖（epithelial cell proliferation）等功能上。
圖4
A549/DDP 細胞和 A549 細胞差異 miRNA 的聚類圖（a）和功能富集圖（b、c）
a：差異 miRNA 聚類圖，紅色代表上調，綠色代表下調；b：差異 miRNA 顯著富集的 KEGG 信號通路；c：差異 miRNA 功能富集的顯著 GO 功能條目；b、c：縱坐標代表信號通路信息和 GO 功能描述信息，橫坐標富集到該 Term 的差異基因數，
2.4 ceRNA 網絡分析
首先對 miRNA 與 lncRNA，mRNA，circRNA 靶點預測。根據 ceRNA 機制，結合已得到的 miRNA-lncRNA/circRNA 及 miRNA（up）-mRNA（down）和 miRNA（down）-mRNA（up），分別預測 lncRNA-mRNA 和 circRNA-mRNA 的 ceRNA 關系對，并構建 miRNA-circRNA-mRNA 和 miRNA-lncRNA-mRNA 等 ceRNA 網絡（圖5 ）。根據拓撲學性質分析，hsa-novel-32-mature（up）和 hsa-miR-365a-5p（up）的 degree 分別為 36 和 25（圖5b），hsa-miR-665（down）和 hsa-miR-370-3p（down）的 degree 分別為 99 和 79（圖5c ），高于本網絡中其它成員的 degree，推測這些 miRNA 在 A549 順鉑耐藥性產生機制中作用顯著。
圖5
ceRNA 三元網絡圖
a：miRNA（down）-circRNA-mRNA（up），三元網絡；b：miRNA（up）-lncRNA-mRNA（down）三元網絡；c：miRNA（down）-lncRNA-mRNA（up）三元網絡；紅色節點表示表達上調，綠色節點表示表達下調，三角形節點代表 miRNA，V 形節點代表 circRNA，圓形節點代表 mRNA，菱形節點代表 lncRNA
對得到 miRNA-lncRNA-mRNA 和 miRNA-circRNA-mRNA 網絡進行整合，僅保留兩者都存在 mRNA 的關系對，構建 miRNA-circRNA-lncRNA-mRNA 四元網絡（圖6）。此網絡包括了 12 個 miRNA，4 個 circRNA，23 個 lncRNA 和 9 個 mRNA 節點。對miRNA-circRNA-lncRNA-mRNA 四元網絡進行拓撲學分析，得到 degree 排名在前 10 名的分子（表1 ）。Degree 越大表示該分子在此網絡中和其它成員相互作用越多，作用越關鍵。
圖6
miRNA-circRNA-lncRNA-mRNA 四元網絡
紅色節點表示表達上調，綠色節點表示表達下調，三角形節點代表 miRNA，V 形節點為 circRNA，菱形節點代表 lncRNA，圓形節點代表 mRNA；藍色連接線為 miRNA-lncRNA 調控關系，橙色連接線為 miRNA-mRNA 關系對，黑色連接線為 lncRNA-mRNA 關系對，粉色連接線為 miRNA-circRNA 調控關系，青色連接線為 circRNA-mRNA 關系對
表1
miRNA-circRNA-lncRNA-mRNA 四元網絡節點信息（degree top10）
3 討論
順鉑是一種重要的肺癌化療藥物，但是固有性及獲得性耐藥嚴重限制其化療效果。據報道，耐藥機制涉及許多細胞系多種基因過表達和/或低表達，從而引起一系列病理變化[22-24]。目前對基因調控在 NSCLC 順鉑耐藥中的認識仍然有限。在本研究中，我們評估了 lncRNA、miRNA 和 cirRNA 在順鉑耐藥細胞系 A549/CDDP 中的表達譜。與野生型 A549 細胞系相比，微陣列分析顯示了一組差異表達的非編碼 RNA，其中 4517 個 lncRNA 和 123 個 cirRNA 以及 145 個 miRNA 在 A549/CDDP 細胞中有差異表達。
基于差異表達的 mRNA，我們進行了 KEGG 和 GO 功能富集分析，以闡明哪些特定的信號通路可能參與了控制 A549/DDP 和 A549 的細胞特征。KEGG 結果表明，差異 mRNA 顯著富集在和癌癥相關的通路上，如：ko05202（transcriptional misregulation in cancers），ko05222（small cell lung cancer），ko04668（TNF 信號通路）和 hsa04115（p53 信號通路）等。TNF 信號通路為 TNF 信號通路，TNF 是一種可以直接殺傷腫瘤細胞，而對正常細胞沒有明顯細胞毒性的細胞因子，此信號通路參與細胞的生存、死亡和分化。肺癌晚期多會發現 TNF 受體表達缺失，人 A549 肺腺癌細胞對 TNF-α 和順鉑的細胞毒性作用具有抗性[22-24]。本研究表明 TNF 信號通路富集了差異 mRNA，我們推測 A549/DDP 耐藥機制可能與 TNF 信號通路異常調控有關。
p53 基因是一種抑癌基因，編碼的蛋白質是一種轉錄因子，其控制著細胞周期的啟動，在絕大多數腫瘤細胞會出現該基因的突變[25]。眾多報道顯示 p53 信號通路與多腫瘤產生的的順鉑耐藥相關。Wang 等[26]根據分子成像研究表明順鉑通過阻斷細胞周期發揮細胞毒作用，可能部分受 p53 信號通路調控。Mg2+依賴性蛋白磷酸酶 1δ（PPM1D）通過減弱細胞周期檢測點激酶 1（Chk1）和 p53 的激活，使卵巢癌細胞對順鉑產生耐藥性[27]。在維生素 C（VC）輔助治療腸癌的研究中發現 VC 可以通過上調 p53 來增強順鉑的敏感性和凋亡[28]。通過分析 A549 與 A549/CDDP 細胞系的全轉錄組測序結果發現差異 mRNA 顯著富集在 p53 信號通路上，提示我們 p53 信號通路可能參與了 A549/CDDP 耐藥機制。
對 lncRNA，miRNA 和 circRNA 表達譜進行綜合分析得到 ceRNA 調控網絡，其中 hsa-miR-125a-5p 和 hsa-miR-125b-5p 在 A549/DPP 抗性細胞中表達下調，是 miRNA-circRNA-lncRNA-mRNA 四元網絡圖中的關鍵 miRNA。成熟的 hsa-miR-125a-5p 來自于 pre-miR-125a 的 5'端，在 NSCLC 癌組織中 hsa-miR-125a-5p 的表達低于鄰近正常肺組織，Spearman 相關檢驗結果表明 hsa-miR-125a-5p 的表達與病理分期或者淋巴結轉移呈負相關[29]。Liu 等[30]研究發現喉癌組織和 Hep-2 喉癌干細胞（Hep-2-cscs）中，miRNA-125a 表達降低，在體內外，miRNA-125a 的功能增強顯著提高了 Hep-2-CSCs 對順鉑的敏感性。Jiang 等[31]發現不同肺癌細胞系中 hsa-miR-125a-5p 的表達均低于人支氣管上皮細胞，并且抑制了 A549 細胞的增殖以及誘導了細胞凋亡，通過分別阻斷野生型的 p53 和 hsa-miR-125a-5p，發現 hsa-miR-125a-5p 通過 p53 依賴通路誘導人肺癌細胞凋亡。另外 Zhang 等[32]也證明 miRNA 125a 可以通過靶向 p53 的 3'‐UTR 翻譯和轉錄抑制 p53 基因。hsa-miR-125b-5p 與 hsa-miR-125a-5p 高度同源，被報道在多種癌細胞中作用顯著，比如肺癌與乳腺癌細胞[33-34]。與卵巢癌細胞系與 OV2008 細胞相比，miR-125b 在耐順鉑細胞株中（C13*）的表達水平升高，上調 microRNA-125b 可顯著抑制順鉑誘導的細胞毒性和凋亡，并可增加 OV2008 和 C13*細胞對順鉑的耐藥性[35]。而我們的研究表明 A549/DPP 抗性細胞中 micRNA-125b 表達是下調的，與上述報道有一定差異，可能是由于不同的細胞種類或者不同的處理方式引起的表達差異。基于前期研究和我們的分析結果推測 hsa-miR-125a-5p 和 hsa-miR-125b-5p 參與了 A549/DPP 抗性機制，可能是逆轉順鉑抗性的潛在靶點。
本研究通過全轉錄組測序比較野生型 A549 細胞和 A549/DPP 順鉑耐藥細胞表達譜的差別，初步篩選了與 NSCLC 順鉑耐藥相關基因，獲得具有潛在意義的耐藥基因，并通過生物信息學分析闡釋其可能機制。但是，我們在研究中獲得的差異表達及機制層面的分析結果尚需進行基礎實驗及臨床大樣本驗證，這也是我們目前正在進行的工作，其深入機制仍有待進一步探索。
利益沖突：無。
作者貢獻：王俊負責實驗設計及部分數據分析；倪耀軍負責論文撰寫及實驗的具體實施；徐忠能和陳勝負責生物信息分析解讀及部分數據分析。
