引用本文： 倪耀军, 徐忠能, 陈胜, 王俊. 非小细胞肺癌细胞 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 越大表示该分子在此网络中和其它成员相互作用越多，作用越关键。
表1
miRNA-circRNA-lncRNA-mRNA 四元网络节点信息（degree top10）
图6
miRNA-circRNA-lncRNA-mRNA 四元网络
红色节点表示表达上调，绿色节点表示表达下调，三角形节点代表 miRNA，V 形节点为 circRNA，菱形节点代表 lncRNA，圆形节点代表 mRNA；蓝色连接线为 miRNA-lncRNA 调控关系，橙色连接线为 miRNA-mRNA 关系对，黑色连接线为 lncRNA-mRNA 关系对，粉色连接线为 miRNA-circRNA 调控关系，青色连接线为 circRNA-mRNA 关系对
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 越大表示该分子在此网络中和其它成员相互作用越多，作用越关键。
表1
miRNA-circRNA-lncRNA-mRNA 四元网络节点信息（degree top10）
图6
miRNA-circRNA-lncRNA-mRNA 四元网络
红色节点表示表达上调，绿色节点表示表达下调，三角形节点代表 miRNA，V 形节点为 circRNA，菱形节点代表 lncRNA，圆形节点代表 mRNA；蓝色连接线为 miRNA-lncRNA 调控关系，橙色连接线为 miRNA-mRNA 关系对，黑色连接线为 lncRNA-mRNA 关系对，粉色连接线为 miRNA-circRNA 调控关系，青色连接线为 circRNA-mRNA 关系对
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 顺铂耐药相关基因，获得具有潜在意义的耐药基因，并通过生物信息学分析阐释其可能机制。但是，我们在研究中获得的差异表达及机制层面的分析结果尚需进行基础实验及临床大样本验证，这也是我们目前正在进行的工作，其深入机制仍有待进一步探索。
利益冲突：无。
作者贡献：王俊负责实验设计及部分数据分析；倪耀军负责论文撰写及实验的具体实施；徐忠能和陈胜负责生物信息分析解读及部分数据分析。
