IMF是评价牛肉品质的指标，可以决定牛肉的风味、多汁性和营养价值等特性[1]，富含IMF的牛肉具有更好的经济效益[2]。日本和牛在24月龄时肌内脂肪含量达24.4%[3]。安格斯牛是欧美地区最具雪花牛肉生产潜力的牛种之一。与和牛和安格斯牛相比，荷斯坦牛肌内脂肪含量匮乏。3种肉牛IMF含量差异较大，但IMF的沉积机制相似。IMF是一种高度可遗传的性状，育种者可以通过遗传改良增加IMF含量[4]。随着基因芯片技术的快速发展，高通量研究特定组织间差异表达基因成为热点，为GEO数据库积累大量有价值的公共测序数据[5]。在医学领域，挖掘GEO数据发现很多与人类疾病相关的基因[6-7]。但是，在畜牧领域只有很少的测序数据被二次挖掘。姚潮刚[8]从GEO数据库筛选3套与猪相关的基因芯片数据，通过整合分析出AMPK、PPAR和脂肪酸代谢3个信号通路与猪IMF沉积相关。何侃等[9]从GEO数据库中选取牛早期胚胎发育过程的基因表达谱数据，显著性及聚类分析显示8-细胞期和16-细胞期在牛的整个早期胚胎发育过程中至关重要。稳健排序整合算法（robust rank aggregation，RRA）是一种功能强大的排序算法，对从每个测序数据集中筛选出已经排序的差异表达基因均分配一个P值来估计显著性概率，P值的大小对应于最终排名列表中相应基因的重要性，P值小、排名位置高、重要性大[10]。王俊[11]应用这种方法从8套胃癌基因芯片数据中挖掘出346个稳健的差异表达基因。Zhang等[12]通过这种方法从5个小鼠转录组数据中，筛选出成脂相关的潜在关键基因。此外，RRA在致病关键基因筛选的研究中应用更加广泛。本研究从GEO数据库中下载两套基因芯片数据集，从包含的3组肌内脂肪和皮下脂肪测序数据中筛选出差异表达基因。研究通过RRA算法获得稳健的差异表达基因，进行GO和KEGG富集分析，构建蛋白质互作网络，筛选关键基因，为进一步研究肌内脂肪沉积提供参考。1资料与方法1.1基因芯片数据集于美国国立生物技术信息中心GEO数据库（https://www.ncbi.nlm.nih.gov/geo/）下载两套基因芯片数据集（GSE 136981和GSE 39006）。GSE 136981数据集包括4头安格斯牛的肌内和皮下脂肪测序数据，平台为GPL 11649；GSE 39006数据集包括3头荷斯坦牛和3头和牛的肌内和皮下脂肪测序数据，平台为GPL 2112。1.2方法1.2.1DEGs的筛选及可视化利用GEO数据库自带分析平台GEO 2R对3组数据测序。以肌内脂肪为试验组、皮下脂肪为对照组进行初步分析。筛选出校正后P0.05且logFC2的基因作为显著上调基因，校正后P0.05且logFC-2的基因作为显著下调基因。利用GraphPad Prism 8.0.2绘制DEGs火山图。1.2.2稳健排序整合算法分析对初步分析获得的3组DEGs按照|logFC|降序排序，应用R语言中的Robust Rank Aggregation包对3组DEGs整合排序，获得稳健DEGs。DEGs定义标准为RRA包默认值。利用GraphPad Prism 8.0.2绘制稳健DEGs热图。1.2.3GO和KEGG分析应用在线分析工具KOBAS 3.0（http://kobas.cbi.pku.edu.cn/kobas3/）和David（https://david.ncifcrf.gov/）对RRA排序的DEGs进行GO和KEGG分析。P0.05表示显著富集。1.2.4蛋白互作网络构建和关键基因筛选将稳健DEGs导入String数据库（https://string-db.org/），构建蛋白互作网络。设定互作评分0.4为筛选阈值。将构建的蛋白互作网络数据导入Cytoscape v3.8.0，利用Cytohubba插件查找蛋白互作网络上的前十关键基因。2结果与分析2.1筛选DEGs（见图1）由图1可知，对GSE 136981数据集处理筛选后，从安格斯牛肌内和皮下脂肪中共筛选出1 086个DEGs，其中，上调基因985个，下调基因101个；对GSE 39006数据集处理筛选后，从和牛中筛选出136个DEGs，其中，上调基因123个，下调基因13个；从荷斯坦牛中筛选出26个DEGs，其中，上调基因23个，下调基因3个。图1差异表达基因火山图10.13557/j.cnki.issn1002-2813.2021.12.016.F1a1（a）安格斯牛10.13557/j.cnki.issn1002-2813.2021.12.016.F1a2（b）和牛10.13557/j.cnki.issn1002-2813.2021.12.016.F1a3（c）荷斯坦牛2.2RRA算法筛选稳健DEGs（见图2）由图2可知，通过RRA算法整合上述3组差异表达基因数据，最终获得18个稳健DEGs，其中显著上调基因17个，分别为骨骼肌肌动蛋白α1（ACTA1）、碳酸酐酶3（CA3）、细胞肌钙蛋白T3（TNNT3）、肌红蛋白（MB）、肌原调节蛋白1（MYOZ1）、肌球蛋白调节轻链（MYPLF）、连接蛋白1（JPH1）、原肌球调节蛋白4（TMOD4）、肌球蛋白结合蛋白C（MYBPC1）、辅肌动蛋白3（ACTN3）、烯醇化酶3（ENO3）、肌球蛋白轻链2（MYL2）、解耦联蛋白3（UCP3）、类平滑肌样蛋白2（SMTNL2）、肌型磷酸果糖激酶（PFKM）、肌节蛋白（MYOT）和伴肌动蛋白（NEB），显著下调基因1个为CD163蛋白样1（CD163L1）。10.13557/j.cnki.issn1002-2813.2021.12.016.F002图2差异表达基因热图注：红色表示上调，蓝色表示下调，空白处表示在相应组别中差异不显著。2.3稳健DEGs的GO和KEGG分析（见图3、图4）由图3可知，上调稳健差异表达基因GO富集分析显示，稳健DEGs主要参与肌肉收缩、肌原纤维组装生物过程，分布在Z盘、横纹肌细丝、肌原纤维和肌球蛋白复合体组分，行使肌动蛋白结合、原肌球蛋白结合、肌肉结构组成的分子功能。10.13557/j.cnki.issn1002-2813.2021.12.016.F003图3DEGs的GO富集结果由图4可知，KEGG通路分析显示，稳健DEGs显著富集在糖酵解/糖异生、氨基酸的生物合成、RNA降解、HIF-1信号通路、碳代谢和黏着斑等14个通路。10.13557/j.cnki.issn1002-2813.2021.12.016.F004图4DEGs的KEGG通路分析结果2.4PPI网络构建和关键基因鉴定（见图5）由图5（a）可知，将RRA运算得到的18个基因导入String数据库构建蛋白互作网络，该网络共有18个节点，47条边。由图5（b）可知，将String数据库下载构建的结果导入Cytoscape软件，利用Cytohubba插件寻找关键基因，选择具有高连通性的10个DEGs作为关键基因。图5DEGs的PPI分析10.13557/j.cnki.issn1002-2813.2021.12.016.F5a1（a）10.13557/j.cnki.issn1002-2813.2021.12.016.F5a2（b）3讨论肌内脂肪含量是评定牛肉肉质优劣的标准，受遗传、营养和管理等多种因素影响[13]。肌内脂肪含量的遗传力相对较高[14]。因此，有必要讨论遗传因素在提高肌内脂肪沉积中的作用。近年来，畜牧工作者通过基因芯片和RNA-seq技术对肌内脂肪沉积的分子机制进行一系列研究，取得了可观成果。但测序平台和样本容量不足使得测序数据不能得到充分的利用，仅使用简单的方法筛选差异基因难免产生偏倚结果。GEO数据库的扩充为不同数据集之间的整合分析提供参考依据[15]。对单个或多个GEO数据集重新分析和集成是对研究成果的进一步扩大。对多个GEO数据集整合分析是将几个数据集所得差异表达基因取其交集，这种模式会忽略在某一个数据集中差异不显著，但是在其他数据集中差异都显著的重要基因，RRA算法可以解决这个问题。RRA算法是一种对多个差异表达基因原有排名再次进行整合排名获得综合性排名的算法[16]。本研究中，对3组DEGs运行RRA算法获得了18个稳健的DEGs，如果运用简单的取交集的方法，MYL2、MYOT、NEB和CD163L1这4个基因将会被遗漏。本研究从GSE 136981和GSE 39006两个数据集3组样本之间筛选出各自分组中的DEGs后，利用RRA算法整合排序获得18个稳健的DEGs，上调基因17个，下调基因1个。GO富集分析显示，上调的稳健DEGs注释到7个条目中，有5个基因分配至Z盘。脂质是脊椎动物Z盘的主要成分[17]。KEGG通路分析显示，上调的稳健DEGs显著富集在糖酵解/糖异生、黏着斑和氨基酸的生物合成等14个信号通路。Huang等[18]比较两个高低肌内脂肪含量品种牛转录组，发现DEGs显著富集在糖酵解/糖异生信号通路。肌内脂肪生成和黏着斑等信号通路有关[19-20]。氨基酸的生物合成通路和鸡胸肉肌内脂肪沉积有关[21]。DEGs与多个信号通路有关，表明这些通路之间存在某种联系。以上信号通路和肌内脂肪生成有关，与本研究结果一致。将稳健DEGs导入String数据库构建PPI互作网络，利用Cytoscape软件筛选到10个高连通性的关键基因。Shin等[22]研究发现，ACTA1基因的表达和韩牛大理石花纹显著相关。田万年等[23]研究发现，ACTA1基因显著影响延边黄牛肌内脂肪含量。Myoz1基因[24]和PFKM基因[25]与猪肌内脂肪含量有关。Dunner等[26]从15个欧洲牛种中筛选16个对肌内脂肪有显著作用的基因，其中包括Myoz1基因，Reverter等[27]也有相似的结论。TMOD4基因在脂肪生成中上调，可能通过调节脂肪生成因子来促进脂肪的生成，并且可能在肌生成和脂肪生成之间起转换作用[28]；Tong等[29]研究发现，MYBPC1高表达可能与和牛大理石花纹有关。UCP3基因和牛[30]、羊[31]和猪[32]肌内脂肪含量有关；MYOT基因的两种突变与和牛肌内脂肪含量显著相关[33]。以上结果均印证了本研究思路和结果的可靠性。4结论本研究基于GEO数据库，通过生物信息学方法筛选出可能参与调控牛肌内脂肪形成的10个稳健DEGs，为牛肌内脂肪形成机制的进一步研究奠定基础。