一个跨物种研究关联基因表达模式的好方法

Original 生信阿拉丁生信阿拉丁 2022-05-16

收录于合集

基本概念

WGCNA（ Weighted correlation network analysis）译为加权基因共表达网络分析，是一种经常用于基因组和系统生物学研究中的网络分析方法。在对一般的数据挖掘分享中，用来发现大量变量的成对关系，例如利用基因表达数据，来构建基因共表达网络，得到基因互作的关系。

在许多WGCNA中，需要对网络中各个模块的特性，以及随着不同情况下的变化，模块的变化情况进行研究。比如，可以研究不同组织或不同物种之间的模块保存程度（module preservation）。模块保存程度统计信息（Preservation Statistic）旨在量化参考数据集和测试数据集之间模块保存情况。WGCNA R包中的函数modulePreservation用于计算不同数据集之间的模块保存统计信息。modulePreservation函数用法可以参考网页：https://rdrr.io/cran/WGCNA/man/modulePreservation.html。

以往的模块保存程度统计信息可以分为四大类：

交叉列表统计信息（Cross-tabulation）：会比较参考数据集和测试数据集中的聚类情况（即模块类别分配），因此它们需要将聚类过程也应用于测试数据。
密度统计信息（Density）：不需要在测试数据集中分配模块标签。
可分离性统计信息（separability）：不需要在测试数据集中分配模块标签。
稳定性统计信息（ stability）：当将一定数量的人工噪声添加到数据中时，稳定性统计通常研究数据集的稳定性，很少用在量化模块保存程度。

在以往分类的基础上，基于量化模块保存度，目前又细分为以下五类：

交叉列表统计信息（Cross-tabulation）；
密度统计信息（Density）；
可分离性统计信息（separability）；
连通性统计信息（Connectivity）：用于确定参考网络基因之间的连接模式是否类似于测试网络中的连接模式；
综合统计信息（Composite）：整合密度统计信息与连通性统计信息，综合评估模块保存度。

这些统计数据通常通过每个数据集中相似对象（通常是基因）的数量或每个数据集之间不相似对象的数量来评估数据集。

modulePreservation函数，一般会统计以下列表的统计信息，其中Type列为每类统计量的类别，Network列为每类统计量使用的网络类型，*input列表示需要输入的数据类型，输入数据需要以下类别：

Lbl：module label；
Adj：general network adjacency；
datX：numeric data from which a correlation network is constructed。

适用的场景

转录组或芯片数据在数据库中存在大量类似的研究结果，那么这些结果之间的相容性如何？WGCNA中modulePreservation函数用于阐述不同实验结果的基因表达模式的关联；
不同物种之间基因表达模式的分析也是多数文章的研究方向，WGCNA中modulePreservation函数同样也适用于阐述不同物种的基因表达模式的关联，比如人与黑猩猩基因表达模块之间的基因保存性，类似的分析有很多，比如人与牛，人与猪等。

modulePreservation 函数实践

以下命令用于演示在黑猩猩数据集中人脑共表达模块的保存情况。

配置R运行环境

启动R之后，设置工作目录并加载所需的包。

 1# Display the current working directory
 2getwd()
 3# If necessary, change the path below to the directory where the data files are stored.
 4# "." means current directory. On Windows use a forward slash / instead of the usual 
 5workingDir = "."
 6setwd(workingDir)
 7# Load the package
 8library(WGCNA);
 9# The following setting is important, do not omit.
10options(stringsAsFactors = FALSE)

导入输入数据

输入表达矩阵

首先下载数据，存放在本地目录中，表达矩阵数据存放在在Dataset 1 (network construction).csv.bz2中。

1# data input
2file = bzfile("Dataset 1 (network construction).csv.bz2")
3dat1 = read.csv(file, header=T)

Dataset 1 (network construction).csv数据展示：

注：为展示方便，该数据共39列，中间省略32列。

表达矩阵包含36个样本，每一行对应一个基因/探针，每一列对应一个样本或辅助信息，我们需要删除辅助信息，并转置表达矩阵用于后续分析。我们只保留BrainvariantH列>0的基因/探针。

1dim(dat1)
2names(dat1)
3datExpr=data.frame(t(dat1[dat1$Brain_variant_H>0,2:39]))
4indexHuman=c(19:36)
5indexChimp=c(1:18)

现在我们设置多组表达数据和相应的模块颜色。

 1# Number of data sets that we work with
 2nSets = 2
 3# Object that will contain the expression data
 4multiExpr = list()
 5multiExpr[[1]] = list(data = datExpr[indexHuman, ])
 6multiExpr[[2]] = list(data = datExpr[indexChimp, ])
 7# Names for the two sets
 8setLabels = c("Human", "Chimp")
 9# Important: components of multiExpr must carry identificating names
10names(multiExpr) = setLabels
11# Display the dimensions of the expression data (if you are confused by this construct, ignore it):
12lapply(multiExpr, lapply, dim)

最后一条命令的输出确认有两个表达数据集，每个数据集包含18个样本中4000个基因的表达测量值。

注意multiExpr的结构：它是一个列表，每个输入数据集必须包含一个称为data的组件并且包含表达数据。外部列表必须具有适当的名称，以便表达式数据可以与我们接下来创建的模块标签匹配。

加载模块标签

WGCNA分析中，通过构建共表达矩阵，计算基因间的相似性，确定基因模块，通常用不同的颜色来代表模块，即每个模块中的基因对应同一种color，方便后续的模块可视化。后续会分享该部分的相关内容。

HumanChimp-OldhamAnalysis-colorHuman-colorChimp-inNetwork.RData中包含黑猩猩表达矩阵中每个基因的模块颜色信息（colorChimp）以及人表达矩阵中每个基因的模块颜色信息（colorHuman）。

1# 加载从网络分析获得的模块标签，并创建一个列表。
2x = load("HumanChimp-OldhamAnalysis-colorHuman-colorChimp-inNetwork.RData") # Create an object (list) holding the module labels for each set:
3colorList = list(colorHuman, colorChimp)
4# Components of the list must be named so that the names can be matched to the names of multiExpr
5names(colorList) = setLabel

注意，multiExpr的名称和multiColor（即这里的colorList）的名称必须对应。

对模块保存程度进行统计分析

这里使用modulePreservation进行模块保存性计算，运行后将结果保存为Rdata，下次可直接加载，节省时间。

1system.time( {
2 mp = modulePreservation(multiExpr, colorList,
3             referenceNetworks = c(1:2),
4             loadPermutedStatistics = FALSE,
5             nPermutations = 200,
6             verbose = 3)
7} )
8
9save(mp, file = "HumanChimp-HumanSpecific-modulePreservation.RData")

该函数输出的是一个嵌套列表，包括quality，preservation，accracy，permutationDetails，referenceSeparability和testSeparability，每个列表分别包含4或5个组成部分，包含：

observed（观察值）；
Z（Z得分）；
log.p（以10为底的p值的对数）；
log.pBonf（Bonferoni校正后的p值的以10为底的对数）；
q（可选，q值）。

每个列表包含Z，log.p，log.pBonf（可选的q值），accuracy中包含的observedOverlapCounts和observedFisherPvalues被构造为2级列表，其中外部组件对应于参考集，内部组件对应于测试集。例如，preservation$observed[[1]] [[2]]包含用于保存测试集中参考集的密度和保存性等保存统计信息，1是参考集，2是测试集。

对模块保存程度进行结果展示

1ref = 1 # Select the human data as reference
2test = 2 # Select the chimp data as test
3statsObs = cbind(mp$quality$observed[[ref]][[test]][, -1], mp$preservation$observed[[ref]][[test]][, -1])
4statsZ = cbind(mp$quality$Z[[ref]][[test]][, -1], mp$preservation$Z[[ref]][[test]][, -1])

我们看主要输出Zsummary评分。

1print(signif(statsZ[, "Zsummary.pres", drop = FALSE],2))

一般文献中认为：

Zsummary>10：高度保存；
2<Zsummary<10：弱保存；
Zsummary<2：无保存证据。

查看不同物种的不同模块之间的overlap基因数与overlap pvalue。

1overlap = cbind(mp$accuracy$observedCounts[[1]][[2]], mp$accuracy$observedFisherPvalues[[1]][[2]])
2print(overlap)

输出结果的前8列为overlap基因数，后8列为 overlap pvalue。

参考文献

Jiang Z, Sun J, Dong H, Luo O, Zheng X, Obergfell C, Tang Y, Bi J, O'Neill R, Ruan Y, Chen J, Tian XC. Transcriptional profiles of bovine in vivo pre-implantation development. BMC Genomics. 2014 Sep 4;15(1):756. doi: 10.1186/1471-2164-15-756. PMID: 25185836; PMCID: PMC4162962.
Langfelder P, Luo R, Oldham MC, Horvath S. Is my network module preserved and reproducible? PLoS Comput Biol. 2011 Jan 20;7(1):e1001057. doi: 10.1371/journal.pcbi.1001057. PMID: 21283776; PMCID: PMC3024255.
Khan FA, Liu H, Zhou H, Wang K, Qamar MTU, Pandupuspitasari NS, Shujun Z. Analysis of Bos taurus and Sus scrofa X and Y chromosome transcriptome highlights reproductive driver genes. Oncotarget. 2017 Apr 13;8(33):54416-54433. doi: 10.18632/oncotarget.17081. PMID: 28903352; PMCID: PMC5589591.
Li Y, Sun J, Ling Y, Ming H, Chen Z, Fang F, Liu Y, Cao H, Ding J, Cao Z, Zhang X, Bondioli K, Jiang Z, Zhang Y. Transcription profiles of oocytes during maturation and embryos during preimplantation development in vivo in the goat. Reprod Fertil Dev. 2020 Apr;32(7):714-725. doi: 10.1071/RD19391. PMID: 32317096.
Xue Z, Huang K, Cai C, Cai L, Jiang CY, Feng Y, Liu Z, Zeng Q, Cheng L, Sun YE, Liu JY, Horvath S, Fan G. Genetic programs in human and mouse early embryos revealed by single-cell RNA sequencing. Nature. 2013 Aug 29;500(7464):593-7. doi: 10.1038/nature12364. Epub 2013 Jul 28. PMID: 23892778; PMCID: PMC4950944.
Miller JA, Horvath S, Geschwind DH. Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways. Proc Natl Acad Sci U S A. 2010 Jul 13;107(28):12698-703. doi: 10.1073/pnas.0914257107. Epub 2010 Jun 25. PMID: 20616000; PMCID: PMC2906579.