人诱导多能干细胞分化为心肌细胞关键基因的生物信息学分析

涂思梅,曲鑫建,2*

(1.大连理工大学生命科学与药学学院,中国辽宁 盘锦 124221;2.广西中医药大学海洋药物研究院，中国广西 南宁 530200)

诱导多能干细胞(induced pluripotent stem cells,iPSCs)是通过体细胞重编程技术获得的具有类似胚胎干细胞(embryonic stem cells,ESCs)自我更新和多向分化潜能的细胞[1～2]。人类iPSCs在细胞治疗、疾病建模和药物开发等领域已经扮演了重要的角色。在当今社会,心血管疾病已经严重地威胁人类的生命健康,而治疗心血管疾病的关键工具和重要因素是获得正常生理功能的心肌细胞。由iPSCs分化而来的心肌细胞,在产生符合治疗和研究心血管疾病的细胞模型方面具有巨大的优势。因此,应用该细胞进行心血管疾病的补偿替代性治疗或开发新型药物具有良好应用前景和重要意义。

相关研究在探索人类多能干细胞(human pluripotent stem cells,hPSCs,包括 iPSCs和 ESCs)向心肌细胞诱导分化的分子机制中发现,Wnt信号、骨形态发生蛋白质(bone morphogenetic proteins,BMPs)和Activin/Nodal等在hPSCs诱导中胚层细胞的产生中具有重要作用,其中,BMP4和Activin A是体外诱导中胚层产生最常用的因子[3];BMP2和成纤维细胞生长因子8(fibroblast growth factor 8,FGF8)可能会协同促进中胚层细胞向心肌细胞分化[4]。另有研究发现,趋化因子及其受体也可以在多能干细胞中发挥重要的作用,其参与多能干细胞的增殖、分化以及多能性维持[5～7]。但是,目前尚没有对诱导iPSCs向心肌细胞分化的重要基因进行系统性筛选与分析的报道。

高通量测序技术的应用,产生了大量的iPSCs向心肌细胞分化的基因表达芯片数据。在本研究中,我们从GEO(Gene Expression Omnibus)数据库中下载了iPSCs向心肌细胞分化的高通量测序数据集,首先,应用R包分析细胞分化过程中的差异表达基因,并对其进行基因本体论(Gene Ontology,GO)与京都基因和基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)信号通路分析;其次,构建差异表达基因的蛋白质-蛋白质互作(protein-protein interaction,PPI)网络;随后,利用R包加权基因共表达网络分析(weighted gene co-expression network analysis,WGCNA)对共表达的基因进行聚类,并将其划分为不同的基因模块;最后,通过Cytoscape对重要基因模块进行可视化和富集分析,以期预测并筛选出iPSCs向心肌细胞分化过程中的关键基因。

1 材料与方法

1.1 基因表达数据的准备

从GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)筛选并下载iPSCs向心肌细胞分化的基因表达谱(GSE137920),获得iPSCs向心肌细胞分化的4个阶段的数据集,分别是未分化的iPSCs(D0)样品3例、中胚层细胞(D2)样品3例、早期心肌细胞(D7)样品3例、心肌细胞(D14)样品3例[8]。

1.2 基因差异表达分析

对基因表达谱GSE137920进行ensembl数据整理,把ensembl ID转换为基因的名称。把Bioconductor 3.12中的 limma 3.46.0、edgeR 3.32.0和pheatmap 1.0.12安装在R 4.0.3中,以识别D0和D2或D7或D14样品之间的差异表达基因,然后进行热图的绘制。差异表达基因获取的标准:错误发现率(false discovery rate,FDR)＜0.05 和|log2(FC)|≥2(FC:fold change),FDR表示矫正后的P值,log2(FC)表示实验组与对照组样品表达量的比值。

1.3 差异表达基因的功能及信号通路分析

GO是从分子功能(molecular function,MF)、生物过程(biological process,BP)和细胞组分(cellular component,CC)3个方面注释基因及其产物的功能。利用GO分析可以找到富集差异基因的GO分类条目,从而找出不同样品中差异基因的生物学功能。KEGG是一个集成了化学、基因组和系统功能信息的数据库,能够把基因及表达信息作为一个整体的网络进行研究,通常被用于细胞信号转导通路的分析。本文首先在R中利用org.Hs.eg.db 3.12.0把基因的名称转换为基因ID,再利用clusterProfiler v3.0.4实现 D0和 D2样品、D0和D7样品以及D0和D14样品差异表达基因的GO功能和KEGG信号通路富集分析,最后利用enrichplot 1.10.1和ggplot 2 3.3.2把输出结果进行可视化。

1.4 共同差异表达基因的获取

利用在线工具Draw Venn Diagram(http://bioinformatics.psb.ugent.be/webtools/Venn/)绘制 D2、D7和D14三个样品中差异表达基因的Venn图,分析得到3个样品中共同差异表达基因。

1.5 共同差异表达基因的PPI网络构建

将共同差异表达的基因导入STRING数据库(http://string-db.org/)[9],获得其所编码蛋白质之间的相互作用关系,PPI网络构建标准:相互作用分数(interaction score)≥0.9。利用Cytoscape 3.7.2中的cytoHubba插件将PPI网络可视化,筛选关键基因。

1.6 加权基因共表达网络构建

WGCNA[10]是根据基因表达模式的不同,挖掘出相似表达模式的基因,通常称为模块(module)。把每个模块与其他模块或外部样本特征进行关联,可用来筛选模块中的核心基因,而核心基因基本是位于调控网络中心的基因,是值得我们优先深入研究和挖掘的对象。本文把iPSCs向心肌细胞分化的4个阶段(D0、D2、D7和D14)的基因表达矩阵作为输入数据,在R中把基因表达量均值大于总体均值的基因作为高表达基因进行WGCNA分析,使用动态剪切树的方法将基因聚类并划分模块,设定最小模块大小的参数为30,将距离较近的模块合并成新的模块,选取剪切值(height)为0.25,将相似性大于0.75的模块合并。最后,使用Cytoscape 3.7.2里的cytoHubba插件把重要基因模块进行可视化,筛选出模块中的核心基因。

1.7 共同差异表达基因和重要模块的通路富集分析

在Cytoscape 3.7.2中安装上Reactome FIPlugIn插件,分析共同差异表达基因和重要模块中基因的相关通路,设定分析的标准为P≤0.01、FDR≤0.05。

2 结果

2.1 iPSCs向心肌细胞分化中差异表达的基因

使用FDR＜0.05和|log2(FC)|≥2作为差异基因的筛选标准,获取D0、D2、D7和D14样品之间的差异表达基因。结果显示:与D0样品相比,D2样品中鉴定出了2 723个差异表达基因(包括959个表达上调和1 764个表达下调的基因);D7样品中鉴定出了5 155个差异表达基因(包括2 707个表达上调和2 448个表达下调的基因);D14样品中鉴定出了6 777个差异表达基因(包括3 871个表达上调和2 906个表达下调的基因)。图1给出了D2、D7和D14样品中前20个表达上调或下调的差异基因的热图。

图1 前20个表达上调或下调的差异基因的热图(A)D2样品中的差异表达基因;(B)D7样品中的差异表达基因;(C)D14样品中的差异表达基因。横坐标为样品,D0样品为未分化的iPSCs(绿框),D2、D7和D14样品分别为中胚层细胞、早期心肌细胞和心肌细胞(红框)。纵坐标表示基因,蓝色表示下调基因,红色表示上调基因。Fig.1 Heatmap of top 20 differential genes with up-or down-regulated expression(A)Differentially expressed genes in D2 samples;(B)Differentially expressed genes in D7 samples;(C)Differentially expressed genes in D14 samples.The sample names are listed on the abscissa.D0 samples are undifferentiated iPSCs(green box),D2,D7 and D14 samples are mesoderm cells,early cardiomyocytes and cardiomyocytes(red box),respectively.The genes are shown on the ordinate.Blue represents down-regulated genes,and red represents up-regulated genes.

2.2 差异表达基因的GO功能分析

GO分析结果显示,D2样品中的差异表达基因有 999个 GO_BP项、58个 GO_CC项、94个GO_MF项(P＜0.05),其中前10项主要包括胚胎器官发育(GO_BP)、含胶原的细胞外基质(GO_CC)和G蛋白偶联受体结合(GO_MF)等(图2A);D7样品中的差异表达基因有1 677个GO_BP项、106个GO_CC项、144个GO_MF项(P＜0.05),其中前10项包括心脏形态发生(GO_BP)、跨膜转运复合体(GO_CC)和受体配体活性(GO_MF)等(图2B);D14样品中的差异表达基因有1 886个GO_BP项、123个GO_CC项、173个GO_MF项(P＜0.05),其中前10项包括心肌组织发育(GO_BP)、离子通道复合物(GO_CC)和信号受体激活剂的活性(GO_MF)等(图 2C)。

图2 差异表达基因的前10个GO功能分析(A)D2样品中差异表达基因的GO功能分析;(B)D7样品中差异表达基因的GO功能分析;(C)D14样品中差异表达基因的GO功能分析。圆的大小表示GO功能富集的基因数量,与圆的大小正相关。同样,GO功能中基因富集的显著性与圆的颜色深度正相关。Qvalue表示用FDR方法校正后的P值。Fig.2 Analysis of top 10 GO functions of differentially expressed genes(A)GO function analysis of differentially expressed genes in D2 samples;(B)GO function analysis of differentially expressed genes in D7 samples;(C)GO function analysis of differentially expressed genes in D14 samples.The circle size represents the number of genes enriched in GO function,which is positively related to the circle size.Similarly,the significance of gene enrichment in GO function is positively related to the color depth of the circle.Qvalue represents the P-value corrected by the FDR method.

2.3 差异表达基因的KEGG信号通路分析

基于KEGG信号通路的分析,D2样品中的差异表达基因在20条信号通路富集(P＜0.05),包括神经活性配体-受体相互作用、细胞因子与细胞因子受体的相互作用、调节干细胞多能性的信号通路、趋化因子信号通路和TGF-β信号通路等。在图3A显示的前20条信号通路中,神经活性配体-受体相互作用的信号通路富集最为显著,有68个差异表达基因富集到此信号通路。D7样品中的差异表达基因在57条信号通路富集(P＜0.05),包括PI3K-Akt信号通路、cGMP-PKG信号通路、心肌细胞中的肾上腺素信号转导、心肌收缩和cAMP信号通路等。在图3B显示的前20条信号通路中,神经活性配体-受体相互作用的信号通路同样富集得最为显著,有113个差异表达基因富集到此信号通路。在D14样品中,差异表达基因在77条信号通路富集(P＜0.05),包括Ras信号通路、肌动蛋白细胞骨架的调节、肥厚型心肌病、MAPK信号通路和扩张型心肌病等。图3C显示了其中的前20条信号通路,还是以神经活性配体-受体相互作用的信号通路富集最为显著,有138个差异表达基因富集到此信号通路。

图3 差异表达基因的前20条KEGG信号通路分析(A)D2样品中差异表达基因的KEGG信号通路分析;(B)D7样品中差异表达基因的KEGG信号通路分析;(C)D14样品中差异表达基因的KEGG信号通路分析。纵轴表示信号通路,横轴表示信号通路中富集的基因数量,颜色表示富集的显著性,Qvalue表示用FDR方法校正后的P值。Fig.3 Analysis of top 20 KEGG signaling pathways of differentially expressed genes(A)KEGG signaling pathway analysis of differentially expressed genes in D2 samples;(B)KEGG signaling pathway analysis of differentially expressed genes in D7 samples;(C)KEGG signaling pathway analysis of differentially expressed genes in D14 samples.The signaling pathways are listed on the vertical axis,and the numbers of genes enriched in the signaling pathway are on the horizontal axis.The color represents the significance of enrichment.Qvalue represents the P-value corrected by the FDR method.

2.4 共同差异表达基因的分析

通过Venn图分析,从3个样品中总共得到917个共同的差异表达基因,包括331个上调差异表达基因(图4A)和586个下调差异表达基因(图4B)。3个样品中的差异表达基因都是与D0样品比较得到的。

图4 Venn图分析得到共同差异表达的基因(A)3个样品中上调差异表达基因的Venn图分析;(B)3个样品中下调差异表达基因的Venn图分析。Fig.4 Common differentially expressed genes obtained by Venn diagram analysis(A)Venn diagram analysis of up-regulated differentially expressed genes in three samples;(B)Venn diagram analysis of downregulated differentially expressed genes in three samples.

2.5 PPI网络的构建和分析

利用STRING数据库对3个样品中共同差异表达的基因构建PPI网络(图5A),然后通过Cytoscape 3.7.2中的cytoHubba插件筛选出PPI网络中的前10个基因(图5B),分别为FPR2、ADCY2、CXCR2、CXCR4、PF4、GALR1、GAL、CXCL5、CXCL6和HCAR3。

图5 共同差异表达基因的PPI分析(A)共同差异表达基因的PPI网络,红色代表上调,绿色代表下调;(B)PPI网络筛选出的前10个基因,其中圆形代表下调基因,菱形代表上调基因,红色、橙色、橘黄色和黄色代表基因连接程度依次递减。Fig.5 PPI analysis of common differentially expressed genes(A)PPI network of common differentially expressed genes.Red represents up-regulation and green represents down-regulation;(B)The top 10 genes screened out in the PPI network.Circles represent down-regulated genes,and diamonds represent up-regulated genes.The colors red,orange,jacinth and yellow represent the decreasing connectivity of genes.

2.6 加权基因共表达网络分析

把所有样品中表达量较低的基因过滤后,选取9 006个基因进行WGCNA分析,将相似性大于0.75的模块合并后最终得到了9个基因模块(图6)。不同的颜色代表不同的模块,树图中每个叶节点代表一个基因,其中密集连接的分支代表接近的基因。表1列出了不同颜色模块所包含的基因数量。其中,绿色基因模块(Green)聚类最显著,包含的基因数目为2 175个;暗绿色基因模块(Darkturquoise)包含的基因数目最少,为133个。

图6 高表达基因的聚类树和模块划分Height为剪切值;Cluster dendrogram为聚类树;Dynamic tree cut为最初得到的模块;Merged dynamic为合并后最终得到的模块。Fig.6 Cluster dendrogram and module division of high expression genesHeight is the cut value;Cluster dendrogram is the cluster tree;Dynamic tree cut is the initial module;Merged dynamic is the final module after the merger.

表1 不同模块中包含的基因数量Table 1 The number of genes contained in different modules

通过WGCNA分析,我们可得到每个模块中直接相互作用的基因。文中选取聚类最显著的绿色基因模块(Green)进行Cytoscape的可视化分析,筛选模块中相互作用的基因(weight＞0.695),结果显示:绿色基因模块(Green)中筛选到399组相互作用的基因。进一步利用Cytoscape软件中的cytoHubba插件对模块中相互作用最多的前50个基因进行可视化,结果显示:在绿色基因模块(Green)中,与其他基因相互作用最多的基因是BMP5(图7)。与其他基因相互作用越多,表明该基因在模块中越接近核心地位。

图7 绿色基因模块中基因的相互作用关系红色表示与其他基因相互作用最多的基因。Fig.7 Interaction between genes in green gene moduleRed indicates the gene that interacts most with other genes.

2.7 共同差异表达基因和重要模块的通路富集

应用Cytoscape软件中的Reactome FIPlugIn插件对共同差异表达基因和绿色模块(Green)基因进行通路分析,结果显示:绿色模块(Green)的基因富集的通路总共有48条,其中5条为翻译后蛋白质修饰、膜转运、脂质代谢、天冬酰胺N-连接糖基化和细胞对压力的反应(表2);共同差异表达基因的富集通路总共有10条,其中5条为GPCR配体结合、肽配体结合受体、RUNX2的转录调控、多能干细胞的转录调控以及趋化因子受体结合趋化因子(表 2)。

表2 绿色基因模块和共同差异表达基因的通路富集Table 2 Pathway enrichment of green gene module and common differentially expressed genes

3 讨论

目前,为探索心血管疾病的细胞治疗途径,获得一种与人体生理状态相接近的心肌细胞模型极其重要。与已建立的人源心肌细胞模型相比,从iPSCs中诱导分化的心肌细胞具有很大优势和前景[11～12]。因此,探究人iPSCs向心肌细胞分化的作用机制以及筛选关键调控基因非常有必要。

本研究在D2、D7和D14三个样品中筛选出的共同差异表达基因有917个,其中包括331个上调基因和586个下调基因(图4)。根据PPI网络中共同差异表达基因的连接程度,我们筛选得到10个关键基因:FPR2、ADCY2、CXCR2、CXCR4、PF4、GALR1、GAL、CXCL5、CXCL6 和 HCAR3(图 5)。其中,FPR2、CXCR2、CXCR4 和 PF4 可能在 iPSCs向心肌细胞分化的过程中发挥重要作用[13～20]。

甲酰基肽受体(formyl peptide receptor,FPR)在人体中包括3个成员,分别为FPR1、FPR2和FPR3,它们是G蛋白偶联趋化因子受体家族的成员[13]。研究证实,FPR2和FPR1通过F-肌动蛋白聚合促进神经干细胞迁移并且向神经元分化[13]。另有研究发现,FPR2和FPR1在活性氧和PI3K-AKT信号通路的介导下,诱导神经干细胞向神经元分化[14]。FPR在骨髓间充质干细胞中表达,而FPR1与N-甲酰-甲硫氨酰-亮氨酰-苯丙氨酸(N-formylmethionyl-leucyl-phenylalanine,fMLP)相互作用,会诱导人的骨髓间充质干细胞向成骨细胞分化[15]。在本研究中,GO和KEGG富集分析显示,FPR2主要富集在第二信使介导的信号、调节细胞因子对生长因子刺激的反应、调控细胞趋化作用以及调控ERK1和ERK2级联反应等,体现了FPR2可能在iPSCs向心肌细胞分化过程中发挥作用。CXC趋化因子受体2(CXC chemokine receptor 2,CXCR2)是一种G蛋白偶联受体。研究发现,在hPSCs中抑制CXCR2表达,会抑制细胞向外胚层分化,导致细胞向中胚层和内胚层分化,随后hPSCs的特性逐渐消失,这表明CXCR2在维持hPSCs的多能性和增殖分化中扮演重要的角色[7,16]。此外,CXCR2及其下游信号在间充质干细胞诱导内皮祖细胞迁移、血管形成等运动中也发挥了关键作用[17]。本研究的GO和KEGG富集分析显示,CXCR2主要富集在心肌细胞凋亡、细胞因子与细胞因子受体的相互作用、血管生成调节、磷脂酶C激活G蛋白偶联受体信号通路等,表明CXCR2可能在iPSCs向心肌细胞分化中发挥作用。CXCR4也是一种G蛋白偶联受体。研究发现,间充质干细胞利用SDF-1/CXCR4信号转导可以诱导c-kit+心脏干细胞的分化以及增殖,该信号通路也可以控制神经干细胞以维持它的干细胞特性,而且被激活的CXCR4能促进人ESCs向神经干细胞分化[18]。在小鼠的iPSCs中增强CXCR4的表达可以提高细胞迁移,并且不改变细胞的干细胞特性[19]。本研究的GO和KEGG富集分析显示,CXCR4主要富集在心脏过程、心脏收缩和调节生长发育等,说明CXCR4可能促进iPSCs向心肌细胞分化。血小板第4因子(platelet factor 4,PF4)是CXC趋化因子家族的成员,是iPSCs分化为心肌细胞的一种生物标志物。研究报道,PF4可以促进iPSCs向心脏分化,这可能与其调节FGF信号有关[20]。在本研究中,GO和KEGG富集分析显示,PF4主要在血管发育调节、cAMP介导信号、造血正向调节和趋化因子的信号通路等富集,提示PF4可能参与iPSCs心肌细胞分化过程。综合以上分析可知,FPR2、CXCR2和CXCR4都是趋化因子的受体,PF4是趋化因子的成员,它们都与细胞的趋化作用有着密切的关系,而且趋化因子的受体几乎都是G蛋白偶联受体,它们都在G蛋白偶联受体配体结合的通路中富集。已有研究证明,趋化因子及其受体可以在多能干细胞的增殖、分化以及多能性维持中发挥重要作用[5～7]。例如:CXCL8(趋化因子)-CXCR2信号通路可以使hPSCs的分化能力提高;CXCL12(趋化因子)-CXCR4信号通路不仅参与hPSCs的迁移活动,还能够增强细胞多能性[21]。因此,我们推测趋化因子及其受体在iPSCs向心肌细胞分化的过程中具有重要作用。

通过对WGCNA的结果进行分析我们发现,最显著的基因模块为绿色基因模块(Green)(图6),在该基因模块中,BMP5与其他基因相互作用的程度最高(图7)。BMP5是BMPs家族的成员,而BMPs是转化生长因子-β超家族中的成员,它可以从细胞内分泌到细胞间隙或者血浆中,从而和靶细胞相应的抗体结合,发挥生物学作用[22]。BMPs及其受体被认为是胚胎发育和器官形成中的重要调节因子,对心血管结构和功能的调节也具有重要作用[23]。其中,BMP4与FGF2协同作用诱导人ESCs向中胚层细胞分化[24],BMP2与FGF8相互作用可能会诱导中胚层向心肌细胞分化[4]。现有研究发现,BMP5参与许多生命活动,如细胞的增殖和凋亡、软骨发育、人类干细胞的分化诱导、人体胚胎发育以及肿瘤的发生和预后[25]。BMP5/7不但可以在体外明显地促进人iPSCs的分化能力,而且可以使神经干细胞分化为中脑多巴胺神经元的能力提高3倍[26]。另外,BMP5还可以调控心肌细胞分化相关基因的表达[27]。本研究发现,BMP5在iPSCs分化为中胚层细胞、早期心肌细胞和心肌细胞的3个阶段都差异表达,而且BMP5主要富集于心肌组织发育、胚胎器官发育、上皮向间充质转换过程以及TGF-β信号通路。因此,我们推测BMP5具有促进iPSCs向心肌细胞分化的功能。

文献研究显示,在本文预测的关键基因中,BMP5和CXCR4已经被证明在干细胞向心肌细胞分化中起重要作用。谭金童[28]成功地利用BMP5重组腺病毒诱导间充质干细胞向心肌样细胞分化。相关研究发现,SDF-1(又称CXCL12)/CXCR4信号通路调节骨髓间充质干细胞分泌血管内皮生长因子(vascular endothelial growth factor,VEGF)和碱性成纤维细胞生长因子(basic fibroblast growth factor,bFGF)[29]。VEGF可以促进iPSCs向心肌细胞分化[30],bFGF能促进ESCs来源心肌细胞的早期成熟[31]。另有研究证实,CXCL12-CXCR4信号通路还参与多种心血管的发育过程[32]。此外,有研究报道,SDF-1/CXCR4信号通路在BMP(BMP2或BMP9)介导的成骨分化中发挥重要的作用[33～35],说明趋化因子及其受体介导的信号通路与BMP可以相互作用,但它们是否在心肌细胞分化过程中发挥作用目前尚不清楚。

本文通过生物信息学分析,虽然筛选得到了参与人iPSCs向心肌细胞分化的几个关键基因,但在该过程中起作用的基因不止于此。在被筛选出来的基因中,FPR2、CXCR2和PF4在iPSCs向心肌细胞分化的过程中可能会发挥重要作用,其中BMP5和CXCR4介导iPSCs向心肌细胞分化的可能性最大,虽然它们已经被证明参与干细胞向心肌细胞分化,但仍需要用实验来进一步证明它们是否参与iPSCs向心肌细胞分化的过程。