应用转录组学筛选全氟辛烷磺酸毒性枢纽基因和关键信号通路的方法

技术领域

本发明涉及生物信息学领域，具体涉及应用转录组学筛选全氟辛烷磺酸毒性枢纽基因和关键信号通路的方法。

背景技术

全氟辛烷磺酸(Perfluorooctane Sulfonate，PFOS)因其极低的表面张力以及疏水、疏油的独特物理性质，曾被广泛应用于民用和工业产品，如表面活性剂、防火涂层和光阻剂等。然而鉴于此类化合物具有毒性、环境持久性、生物蓄积性和长距离迁移等特性，联合国环境规划署于2009年将PFOS及其盐类和前体化合物列为“斯德哥尔摩公约”中的新型持久性有机污染物。已有研究表明PFOS对动物表现出肝毒性、睾丸毒性、胰腺毒性、神经毒性、免疫毒性、生殖和发育毒性等多种毒性作用，并且会诱导相关肿瘤的发生。

转录组学是在整体水平研究生物体中的基因转录情况和转录调控规律的技术。它从RNA水平研究基因表达的情况，广泛测定出基因的调控情况，并运用生物信息学方法，揭示基因调控的关键信号通路。当生物体受到外界干扰刺激或处于不同的生理状态时，基因表达受到调控，转录组能够直观地反应生物体的应答状况。

发明内容

本发明所要解决的技术问题是如何筛选PFOS毒性枢纽基因和关键信号通路和/或如何鉴定PFOS对动物或人类造成的影响。

为了解决上述技术问题，本发明提供了一种筛选或辅助筛选PFOS毒性枢纽基因和关键信号通路的方法。

所述方法包括如下步骤：

A1)将实验动物分为对照组和PFOS暴露组，对所述对照组进行空白溶剂给药处理，对所述PFOS暴露组进行PFOS溶液给药处理，分别得到给药处理后的对照组和PFOS暴露组实验动物；所述PFOS溶液的溶质为PFOS，溶剂为所述空白溶剂；

A2)对A1)中得到的所述给药处理的对照组和PFOS暴露组实验动物取血，提取全血RNA后进行测序，得到对照组和PFOS暴露组的转录组原始数据；

A3)对所述转录组原始数据进行处理，获得对照组和PFOS暴露组(所有样本)的基因表达量矩阵；对所述对照组和PFOS暴露组(所有样本)的基因表达矩阵进行处理和筛选，得到所述对照组和PFOS暴露组实验动物的差异表达基因；

A4)对A3)中所述差异表达基因进行基因本体论(GO)富集分析和KEGG通路富集分析，得到GO富集分析结果和/或KEGG通路富集分析结果；

A5)对A3)中所述差异表达基因构建蛋白相互作用网络，获得PFOS毒性相关枢纽基因；

A6)对A3)中所述基因表达量矩阵进行基因集富集分析(GSEA)分析，得到GSEA分析结果；

A7)选取A4)中得到的所述KEGG通路富集分析结果和A6)中得到的所述GSEA分析结果中的共有通路为PFOS毒性相关的关键信号通路。

上述测序可为高通量测序。

上述方法A1)步骤中所述空白溶剂可为体积比为2～5％的Tween 20或Tween 80的水溶液。

上述方法A1)步骤中所述给药处理条件可为连续灌胃给药7～28天。PFOS给药量可为2.5～10mg/kg/天。

上述方法A1)步骤中所述给药处理条件可为连续灌胃给药28天。PFOS给药量可为2.5mg/kg/天。

上述方法A3)步骤中所述筛选的筛选条件可为：Fold Change>2且Q value<0.001。所述Fold Change可为所述差异表达基因在PFOS暴露组中的表达量与其在对照组中的表达量的比值。所述Q value可为衡量错误发现率的指标。

上述方法A3)步骤中所述处理的方法可为使用FPKM作为单位，计算所述对照组和所述PFOS暴露组中基因的表达水平。所述筛选的方法可为使用DEGseq分析得出所述对照组和所述PFOS暴露组的差异表达基因。

上述处理的方法也可以使用RPKM等作为单位计算所述对照组和所述PFOS暴露组中基因的表达水平。

上文所述FPKM代表每千个碱基的转录每百万映射读取的碎片(fragments)。

上述方法A4)步骤中所述进行GO富集分析和KEGG通路富集分析的方法可为使用R语言clusterProfiler软件包中的enrichGO函数和enrichKEGG函数进行分析。

所述GO富集分析还包括将所述enrichGO函数的分析结果使用dotplot函数进行可视化。所述KEGG通路富集分析还包括将所述enrichKEGG函数分析结果使用emapplot函数进行可视化。

上述方法A4)步骤中所述进行GO富集分析和KEGG通路富集分析的方法在差异基因数量小于2000个时也可使用DAVID数据库进行分析。所述使用DAVID数据库进行分析具体可为将A3)中所述差异表达基因导入DAVID数据库，选择导入基因的类型和物种信息后，进行富集分析。所述富集分析结果可导入至R语言中进行可视化。所述可视化可为使用dotplot函数对所述GO富集结果进行可视化。所述可视化也可为使用emapplot函数对所述KEGG富集结果进行可视化。

上述方法中，A5)步骤中所述构建蛋白相互作用网络的方法可为使用R语言STRINGdb软件包进行构建；然后使用CytoHubba插件选择所述枢纽基因。

所述枢纽基因可为MNC值前10个基因。所述MNC值为Maximum NeighborhoodComponent(最大邻域分量)。所述使用STRINGdb包进行构建的过程可为：设置参数为：version＝11.0，species＝10116，score_threshold＝900，通过get_interactions函数获得蛋白相互作用关系。所述构建过程还包括使用Cytoscape软件将所述蛋白相互作用关系可视化。

上述方法中，当A3)步骤中差异基因数量在2000个以下时，A5)步骤中所述构建蛋白相互作用网络的方法也可为使用STRING数据库构建。所述使用STRING数据库进行分析具体可为将A3)中所述差异表达基因导入STRING数据库得到所述蛋白相互作用网络。

上述方法A6)步骤中所述进行GSEA分析的方法可为使用GSEA软件进行分析。在使用所述GSEA软件进行分析前需要先确定基因集。所述基因集可以直接利用在GSEA软件中给出的预定基因集。也可以通过R语言中的biomatR包建立基因集。

所述通过biomaRt包建立基因集的步骤包括抓取实验动物的KEGG数据库数据，建立各个KEGG通路的基因集本地数据库，将数据库导出并转为GSEA软件所学的gmt格式文件，同时将基因表达矩阵转为gct格式文件，并建立cls分组信息文件。

所述GSEA软件分析时参数设置具体可为Number of permutations＝1000，Permutation type＝gene set，Metric for ranking genes＝Signal2Noise，进行分析，筛选FDR<0.25，pvalue<0.1的结果得到显著富集通路。

上文所述GSEA软件可在网站http://www.gsea-msigdb.org/gsea/index.jsp进行下载。

上文所述实验动物可为大鼠或小鼠。

为了解决上述技术问题，本发明还提供了一种存储有计算机程序的计算机可读存储介质。所述计算机程序可使计算机执行如上述方法中A3)-A7)所述的步骤。

上述方法和计算机可存储介质的目的是非疾病诊断目的、非疾病预后目的和非疾病治疗目的。

本发明应用转录组学筛选PFOS毒性枢纽基因和关键信号通路。设定空白对照组和PFOS暴露组，提取生物样品的转录组，通过高通量测序仪进行转录组测序，并将获得的测序数据进行生物信息学分析，得到PFOS毒性相关的枢纽基因和关键信号通路。

附图说明

图1为基因本体论富集分析图。p.adjust代表调整后的p值；Count代表基因数量。

图2为KEGG富集分析图。p.adjust代表调整后的p值；size代表基因数量。

具体实施方式

下面结合具体实施方式对本发明进行进一步的详细描述，给出的实施例仅为了阐明本发明，而不是为了限制本发明的范围。以下提供的实施例可作为本技术领域普通技术人员进行进一步改进的指南，并不以任何方式构成对本发明的限制。

下述实施例中的实验方法，如无特殊说明，均为常规方法，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。

下述实施例中所用试剂的配制方法如下：

含2％Tween 20的水溶液：将Tween 20(购自Sigma-Aldrich)溶于无菌水中，得到Tween 20体积比为2％的水溶液。

质量浓度为1mg/mL的PFOS溶液的配制：将PFOS(购自Sigma-Aldrich)溶于含2％Tween 20的水溶液中，得到质量浓度为1mg/mL的PFOS溶液。

实施例一、应用转录组学评鉴PFOS毒性机制方法的建立和结果分析

1.转录组学数据的获得

1.1 PFOS处理

下述实验包含三个重复。

(1)选择6-8周龄的Sprague-Dawley大鼠(购自北京维通利华实验动物技术有限公司)，共8只，雌雄各4只。8只大鼠出生日期相同或者相近，同性间体重相近，记录年龄(雄鼠6周，雌鼠6周)和体重(雄鼠体重281±13g，雌鼠体重190±6g)。将8只大鼠随机分为2组，依次为对照组和PFOS暴露组，每组雌雄各2只。对照组和PFOS暴露组大鼠不可养在一个笼子中。

(2)大鼠饲养1周后进行实验。对照组每只大鼠每天给药一次含2％Tween 20的水溶液，给药体积为当天暴露组大鼠给药体积的平均值；PFOS暴露组每只大鼠每天给药一次质量浓度为1mg/mL的PFOS溶液(每天新鲜配置)，给药剂量为2.5mg PFOS/kg大鼠体重。两个处理组每天上午同一时间灌胃给药，连续给药28天。

1.2转录组学数据获得

在第29天，于同一时间将对照组和PFOS暴露组的8只大鼠使用二氧化碳处死，然后分别对每只大鼠进行心脏采血，得到对照组和PFOS暴露组各个样本的全血。对对照组和PFOS暴露组各个样本的全血进行RNA提取后，使用BGISEQ-500进行高通量测序，分别得到对照组和PFOS暴露组各个样本的转录组数据。

2.筛选差异基因

对步骤1产生的原始转录组数据文件中低质量的reads进行过滤，得到cleandata。使用HISAT 2软件(http://daehwankimlab.github.io/hisat2)对clean data进行参考基因组比对，使用Bowtie2(http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)将clean data比对到大鼠参考基因序列(版本：rn6)，使用RESM(http://deweylab.biostat.wisc.edu/rsem/rsem-calculate-expression.html)计算基因表达量，基因表达量单位用FPKM表示，然后通过计算对照组和PFOS暴露组各个样本的基因表达水平得到对照组和PFOS暴露组的基因表达量矩阵。

通过DEGseq(http://bioconductor.org/packages/release/bioc/html/DEGseq.html)进行差异基因(PFOS暴露组相对于对照组的转录组数据中存在差异表达的基因)分析，差异基因筛选设置条件为Fold Change(差异表达倍数)>2且Q value<0.001。通过筛选总共得到6838个PFOS暴露组相对于对照组的转录组数据中存在差异表达的基因，部分差异基因如表1所示。

表1部分差异表达基因信息

注：“Entrz ID”为NCBI Entrez数据库中的基因标识符；“Symbol ID”为基因名称；“Control Group FPKM”为对照组基因表达水平；“Experimental Group FPKM”为PFOS暴露组基因表达水平；“Fold Change”为差异表达倍数；“Q value”为衡量错误发现率的指标。

3.差异基因分析

3.1差异基因基因本体论(GO)富集分析和KEGG通路富集分析。

将上述步骤2中得到的6838个差异表达基因，导入至R语言中，使用clusterProfiler包中的enrichGO函数对差异基因进行GO富集分析，选择物种信息为大鼠，获得GO的富集结果(包括生物过程、细胞组件和分子功能三部分富集结果)，设置参数pvalueCutoff＝0.05，ont＝ALL，pAdjustMethod＝fdr；随后使用dotplot函数对GO富集结果进行可视化，GO富集可视化结果如图1所示。

使用enrichKEGG函数对差异基因进行KEGG通路富集，选择物种信息为大鼠，设置参数pvalueCutoff＝0.05，organism＝rno，pAdjustMethod＝fdr；使用emapplot函数对通路富集结果进行可视化，KEGG通路富集可视化结果如图2所示，显著富集的通路有为B细胞受体信号通路、Fcγ受体诱导的吞噬、剪接体、NF-κB信号通路、急性髓系白血病、内质网中的蛋白质加工、乙型肝炎、RNA转运、DNA修复、核苷酸切除修复、内吞和结肠癌。

3.2.差异基因蛋白相互作用网络构建分析

将上述步骤2中共得到6838个差异表达基因，使用R语言STRINGdb包，通过get_interactions函数获得蛋白相互作用关系，设置参数为version＝11.0，species＝10116，score_threshold＝900，将结果导入至Cytoscape软件中进行可视化，得到可视化蛋白相互作用网络。利用CytoHubba插件筛选蛋白相互作用网络中的枢纽基因，具体为利用CytoHubba提供的拓扑分析方法，对筛选出的基因进行从大到小排序，选取排名靠前的多个基因作为枢纽基因。筛选MNC值(Maximum Neighborhood Component)排在前10个的基因作为枢纽基因。枢纽基因筛选结果如表2所示。

表2 MNC值排在前10的枢纽基因信息

3.3.差异基因基因集富集分析(Gene Set Enrichment Analysis,GSEA)

为了更准确地筛选出关键通路，利用GSEA软件(http://www.gsea-msigdb.org/gsea/index.jsp)进行基因集富集分析，通过R语言biomaRt包抓取大鼠的KEGG数据库数据，建立各个KEGG通路的基因集本地数据库，将数据库导出并转为GSEA软件所需要的gmt格式文件，同时将步骤2中得到的6838个差异基因表达矩阵转为gct格式文件，并建立cls分组信息文件，将三个文件(gmt格式文件、gct格式文件和cls分组信息文件)导入GSEA软件后进行分析得到GSEA富集分析结果，设置参数Number of permutations＝1000，Permutationtype＝gene set，Metric for ranking genes＝Signal2Noise。获得的GSEA富集分析结果中，认为FDR<0.25，pvalue<0.1的结果认为显著富集。

3.4差异基因关键信号通路分析

将3.3中得到的GSEA富集分析结果与3.1中得到的KEGG通路富集分析结果取交集，得到6条关键通路，分别为B细胞受体信号通路、Fcγ受体诱导的吞噬、剪接体、NF-κB信号通路、急性髓系白血病、内质网中的蛋白质加工和乙型肝炎。6条关键通路的信息如表3所示。

表3差异基因分析得到的关键信号通路

注：FDR代表错误发现率；ES代表富集得分。

由上述分析结果可知，PFOS毒性对大鼠影响较为重要的枢纽基因为Ncbp1、Ncbp2、Cdc5l、Polr2b、Polr2j、Polr2g、Polr2d、Gtf2f1、Gtf2f2和Dync1h1；PFOS毒性对大鼠影响的关键信号通路为B细胞受体信号通路、Fcγ受体诱导的吞噬、剪接体、NF-κB信号通路、急性髓系白血病、内质网中的蛋白质加工和乙型肝炎。因此应用本发明中所提供的应用转录组学筛选PFOS毒性枢纽基因和关键信号通路的方法，可筛选得到PFOS对于动物模型毒性关键的枢纽基因和信号通路，为鉴定PFOS对动物或人类造成的严重影响和/或PFOS对动物或人类的毒性机制以及开发针对PFOS毒性的相关药物提供帮助。

以上对本发明进行了详述。对于本领域技术人员来说，在不脱离本发明的宗旨和范围，以及无需进行不必要的实验情况下，可在等同参数、浓度和条件下，在较宽范围内实施本发明。虽然本发明给出了特殊的实施例，应该理解为，可以对本发明作进一步的改进。总之，按本发明的原理，本申请欲包括任何变更、用途或对本发明的改进，包括脱离了本申请中已公开范围，而用本领域已知的常规技术进行的改变。按以下附带的权利要求的范围，可以进行一些基本特征的应用。