基于组织特异性基因集权重的分析方法研究

田 甜, 王化琨

(黑龙江大学 数学科学学院, 哈尔滨 150080)

0 引 言

基因集检验或通路分析是一种功能强大且被广泛采用的基因分析方法，可用于分析和解释高维基因组数据[1-2]。基因集检验可使研究人员从研究单个基因变量的水平拓展到基因集合的多变量水平，进而探索具有生物学意义的基因组关联情况，如：参与特定代谢通路的关联基因。 基于特定功能相关的基因组变量进行的统计学检验比基于单个基因变量的检验有更多优势，包括改善的统计功效、更直观的生物学解释等[3-5]。鉴于这些优势，研究人员在过去的15年间投入了大量的精力，开发出很多有效的基因集检验方法[6-9]，如ROAST方法，它是自检验方法，运用了一种基于蒙特卡洛(Monte Carlo)多变量技术原理，对数据进行旋转来替代传统的交换排列阵的方法。再如CAMERA方法，它是竞争性检验方法，利用线性模型估算出基因集中平均的基因关联。在改善基因集检验方法的同时，研究人员也在建立大型基因集公共存储库方面取得了不错的进展, 如最常用的GO数据库[10-12]。但是基因集检验在实际应用中仍然受到一些限制，包括基因注释的质量、统计能力和组织特异性等方面。

组织特异性的存在是取决于构成该组织主要成分的细胞性质。随着微阵列技术的成熟运用，通过基因表达的全基因组分析，人们发现基因表达的功能与组织密切相关，并且一些普遍存在的生物过程也发生在特定的组织当中。随着人类蛋白图谱(HPA)[13]和基因组织表达工程(GTEx)[14]等项目的开展，人们对特定组织内的基因活动的认识也逐渐加深。现在，通过质谱分析和免疫组织化学等技术，可以对完整的人类基因组序列分析中鉴定出的约20 687个蛋白质编码基因的组织特异性活动进行研究。如在HPA中所述，在所有蛋白质编码基因中，有34%的蛋白质在至少一种组织中表达升高，其中有17%在特定组织中的mRNA水平至少是所有组织中平均水平的5倍，还有约44%的蛋白编码基因在所有组织中表达，并且在任何组织中都没有升高。

尽管对人类基因的组织特异性活动的描述花费了大量精力，但是在基因集检验中却很少利用组织特异性信息，如分子签名数据库[12](MSigDB)完全缺乏有关基因集的组织特异性和基因集的注释信息。由于组织特定版本的基因集不易创建，通常使用的作法为：在执行基因集检验时使用通用的基因集进行，而不去考虑实验组织的类型。如果基因在组织中都是无所不在的表达，忽视组织特异性进行基因集检验对结果影响不大。但除了持家基因是在所有组织中具有类似的表达水平之外，大多数基因在不同组织中表达水平是有差异性的，而在此条件下进行基因集检验，会在一定程度上提高Ⅰ型和Ⅱ型错误率。为了克服这种基因集检验的缺点，本研究使用来自人类蛋白图谱的组织特异性基因活性信息和分子签名数据库中的所有过滤后的基因集合生成了组织特异性基因集权重的集合[15]，并且利用这些组织特性的基因集权重对p值进行加权，以这种方式进行组织特异性的基因集检验。

1 材料与方法

1.1 数据来源

(1)人类蛋白质图谱HPA[13]

人类蛋白质图谱是一项基于瑞典的科研计划，它始于2003年，旨在利用多种组学技术的整合来绘制细胞、组织和器官中所有人类蛋白质的图谱，包括基于抗体的成像，基于质谱的蛋白质组学、转录组学和系统生物学。本文主要是使用组织图谱，该图谱是基于对37种主要的人体正常组织类型的RNA(RNA-seq) 进行深度测序，并在包含44种不同组织类型的组织微阵列上进行免疫组织化学分析，它包含人类基因在mRNA和蛋白质水平上的表达谱信息。本文关于人类蛋白质编码基因的组织特异性活性信息是从v18.1版本的人类蛋白质图谱中下载的。

① IHC蛋白丰度数据：是基于免疫组化学和组织微阵列的蛋白在正常人体组织的表达谱，该数据包括基因标识符、组织名称、注释细胞类型、表达值和表达值的基因可靠性，根据表中的蛋白表达值计算数值得分记为ProteinScore(其中未检测到=0，低=0.5，中等=1.0,高=2.0)。

② RNA-seq数据：是基于RNA-seq的37个组织的RNA水平，该数据包括基因标识符、分析样本(组织)和每百万转录本。根据表中的TPM值，计算折叠后的TPM值，记为FoldAboveMean, 该组织中的平均TPM值表示为mean(TPM),则FoldAboveMean=TPM/mean(TPM)。然后根据基因名和组织名将两个列表结合成一个表。

(2) 分子签名数据库MSigDB[12]

注释齐全的基因集代表了生物过程的整体，对于解释大规模基因组数据至关重要，分子签名数据库是此类集合中使用最广泛的存储库之一。分子签名数据库是将人类基因从位置、功能、代谢途径和靶标结合等多种角度出发，构建了许多基因集合，其中的一个基因集合中包含了具有相近位置或类似功能的许多基因。分子签名数据库涵盖了很多种类的基因集合和更广泛的基因集来源和类型，包括从原始研究出版物中提取的签名以及从GO[10]和KEGG[11]等专业资源中提取的整套集合。该数据库中的基因集是通过人工筛选和自动计算两种方法获得的。最初的分子签名数据库是在2005年发布在GSEA软件上，共有1 325套，并且该数据库是不断在更新的。

本文使用的基因集是从分子签名数据库(MSigDB)的v7.1版本下载的，包含了13个不同类别的集合，共有25 824个基因集。

1.2 分析方法

本文使用了来自人类蛋白质图谱的组织特异性基因活性信息来计算分子签名数据库中集合的组织特异性基因集权重，该方法如图1所示。由图可知，来自分子签名数据库的集合被表示为一个矩阵，行代表基因集，列代表基因，若基因和基因集之间存在注释，则元素记为1，然后使用人类蛋白质图谱中在不同组织的基因活性信息，根据公式(1)计算出组织特异性基因权重，最后将分子签名数据库中的基因集合作为输入，利用下述步骤，得到组织特异性基因集权重。

图1 组织特异性基因集权重计算的方法示意图

1.2.1 组织特异性基因权重的计算

(1)

式中：ei,t表示基因i在组织t中的RNA-seq数据折叠后的TPM值，单位为每千碱基转录本片段数/百万片段映射，若基因i在组织t中的RNA-sep数据缺失，则ei,t记为1；ai,t表示基因i在组织t中的IHC蛋白丰度数据，若基因i在组织t中的蛋白质丰度数据缺失，则ai,t记为1。式(1)生成了组织特异性基因权重，需要在蛋白质和RNA水平的证据才能产生非零值。

1.2.2 组织特异性基因集权重的计算

(2)

式中m表示注释到基因集j的基因，m={i=1,2,…,p}且Gj,i=1，其中Gj,i=1表示基因i注释到基因集j中，|m|为此集合的尺寸。mc为补集，mc={1,2,…,p}且Gj,i=0，其中Gj,i=0表示基因i未注释到基因集j中，|mc|为补集的尺寸。

采用单边双样本的t检验，得到了p值，则有：

(3)

式中pvalj,t是t检验中得到的p值，是以多大的误差拒绝原假设，原假设为组织t中注释到基因集j的基因的平均权重等于不在基因集j中的平均权重。取p值的负对数就得到了组织特异性的基因集权重,并且本文还在假设检验前对MSigDB数据库中的13个不同类别的基因集都进行了过滤，保留了基因集中基因个数在10～200之间的基因集合。

2 结果与分析

2.1 组织特异性基因集权重

在假设检验之前对分子签名数据库中的13个不同类别基因集都进行了过滤，如表1所示。这是为了解决某些基因集的组织特异性权重相对于其他基因集的组织特异性权重较大的问题，在这种情况下，一些不重要的p值会在基因集检验下生成显著的FDR值，从而影响结果的准确性，所以将基因集权重离散化，即对基因集进行过滤。

表1 MSigDB v7.1版本的所有13个不同集合

使用上述分析方法，为分子签名数据库中的13个不同类别基因集合和人类蛋白质图谱所支持的37种人类组织类型，生成了组织特异性基因集权重。人类蛋白质图谱支持的组织类型有脂肪组织、大脑皮层、骨髓、肝脏、肾脏等37种组织类型。这些组织特异性的基因集权重，在下文中还会使用到。

2.2 组织特异性基因集检验

当仅检验一个基因集时，p值是统计显著性的适当度量，但是当检验包含数千个基因集的大型基因集时，在基因集中可能会出现看似很高的假阳性的p值，这就被称为多重假设检验的问题。

2.2.1 白血病数据

GSE131184数据集是来源美国国家生物技术信息中心(NCBI)网站中的基因表达数据库(GEO),平台号是GPL570，该数据集是从急性髓系白血病(AML)或T急性淋巴细胞白血病(T-all)患者的骨髓样本中获取的总RNA，样品分为两个批次运行(标记为L或M)，共有125个样本。

2.2.2 重度抑郁症数据

GSE54563数据集是来源美国国家生物技术信息中心(NCBI)网站中的基因表达数据库(GEO),平台号是GPL6947，该数据集是关于重度抑郁症的，表达数据是从人类大脑前扣带皮层组织中获取的，包含了25对共50个样本的对照样本和重度抑郁症样本。

2.2.3 二型糖尿病数据

GSE73034数据集是来源美国国家生物技术信息中心(NCBI)网站中的基因表达数据库(GEO),平台号是GPL6480，该数据集是关于二型糖尿病的，该数据比较了瘦、肥胖胰岛素敏感(OIS)、肥胖胰岛素抵抗(OIR)和肥胖T2D患者肌肉活检的基因表达差异。表达数据从人类骨骼肌中获取，每组有7个样本，共28个样本。

2.2.4 结果

以上三个数据集的基因表达谱不能直接使用，所以利用R语言中的Bioconductor进行数据预处理，将表达谱中的探针集注释到对应的基因上，得到处理好后的归一化基因表达数据，然后使用两阶段的竞争性基因集检验方法CAMERA进行基因集检验[19]，该方法可以随基因集成员之间的相关性进行调整，这种竞争性的基因集检验形式假定了基因水平权重的独立性，并且适用于检验包含许多基因集的大型数据库中哪些基因集相对于其它基因集整体变化更为显著，也适用于找出具有意义的基因集，此方法由Limma包中的Camera功能实现[20]。本次基因集检验使用到的基因集合是分子签名数据库中过滤后的C2.CP集合，并使用了Camera中的默认设置，然后将使用Camera方法后得到的FDR值与使用BH方法控制的加权FDR值进行比较，该结果见表3。这三种类型的疾病在q值小于0.2的情况下的加权的FDR分析中产生了更多的发现，如表中所示白血病基因表达数据在未加权的情况下有32个发现，加权后有44个发现；重度抑郁症基因表达数据在未加权的情况下有1个发现，而加权后有18个发现；二型糖尿病基因表达数据集在未加权的情况下产生了1个发现，而加权后有2个发现。在加权后wFDR分析产生了额外的基因集发现，一般来说，这个结果在生物学上是合理的。

p值加权在以下两种情况下最有效，一是当基因集检验的目的是识别目标组织中起重要作用的基因集，二是在基因集检验中因变量与被分析组织的正常活动和组织特异性过程显著关联。由于权重较大的基因集更能反映正常活动和组织特异性，在这两种情况下，p值加权能更好地提高统计功效。

表3 组织特异性基因集检验结果

表4列出了使用上述方法进行wFDR分析骨骼肌相对于二型糖尿病产生的两个基因集发现。

表4 在骨骼肌中与二型糖尿病有关的通路

表5 肺(左)、肾脏(右)在C2.CP中权重前十的基因集

基因集权重REACTOME_ORGANIC_CATION_ANION_ZWITTERION_TRANSPORT128REACTOME_TRANSPORT_OF_BILE_SALTS_AND_ORGANIC_ACIDS_METAL_IONS_AND_AMINE_COMPOUNDS89KEGG_GLYCINE_SERINE_AND_THREONINE_METABOLISM77REACTOME_GLYOXYLATE_METABOLISM_AND_GLYCINE_DEGRADATION71REACTOME_TRANSPORT_OF_INORGANIC_CATIONS_ANIONS_AND_AMINO_ACIDS_OLIGOPEPTIDES58REACTOME_SLC_TRANSPORTER_DISORDERS49KEGG_ARGININE_AND_PROLINE_METABOLISM39REACTOME_NA_CL_DEPENDENT_NEUROTRANSMITTER_TRANSPORTERS38REACTOME_PYRIMIDINE_CATABOLISM33KEGG_PROXIMAL_TUBULE_BICARBONATE_RECLAMATION29

2.3 人类组织的功能表征

组织特异性基因集权重可以用来分析相关人类组织的功能表征。具体说，对分子签名数据库中的某类集合如C2.CP集合，按照1.2节所述步骤，对分配给给定组织的每个基因集合的权重进行排序，权重大的基因集更能反映组织中活跃的主要生物过程。如对分子签名数据库中的典型通路(C2.CP)中的肺组织和肾脏组织进行分析，表5显示了给定肺组织、肾脏组织在分子签名数据库中的C2.CP集合中排名前十的基因集合。对给定的分子签名数据库中的集合名可以通过GSEA软件查找，分子签名数据库中对该集合的描述，或者通过外部链接到其他数据库(如Reactome信号通路数据库)的描述信息。Reactome是一个免费开源的通路数据库，提供直观的生物信息学工具，用于可视化、解释和分析通路相关知识，以支持基础研究、基因组分析、建模和系统生物学研究等。根据这些数据库中的信息可知，权重高的基因集合捕捉到了与该组织有关的已知生物学特性和过程，识别出了肺组织中的代谢通路，如REACTOME_SURFACTANT_METABOLISM是与肺组织表面活性物质有关的新陈代谢通路，也识别出与肾脏有关的运输通路，如REACTOME_ORGANIC_CATION_ANION_ZWITTERION_TRANSPORT是有机阳离子/阴离子/两性离子转运通路，这些转运蛋白在肾脏组织中表达。

3 结 论

使用人类蛋白质图谱的组织特异性基因活性信息和过滤后的分子签名数据库中的基因集生成了组织特异性的基因集权重。为了避免权重较大的基因集影响结果的准确性，将分子签名数据中的基因集合进行过滤。利用组织特异性基因集权重对p值进行加权并用BH方法提供FDR控制，进行wFDR分析，在不同组织的三种疾病下加权的FDR一共有64个发现，而未加权的FDR一共有34个发现。在p值小于0.2的假设下，使用组织特异性基因集权重对p值进行加权能产生更多的发现。使用组织特异性基因集权重可以提高基因集检验的统计功效，并且也可以提高在高维基因组数据的实验中识别生物学上有关的基因集关联的信息。