基因表达模式识别专题

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

2021-06-14  本文已影响0人  生信阿拉丁

作者：oct
审稿：童蒙
编辑：amethyst

一、基本概念

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

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

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

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

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

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

二、适用的场景

转录组或芯片数据在数据库中存在大量类似的研究结果，那么这些结果之间的相容性如何？WGCNA中modulePreservation函数用于阐述不同实验结果的基因表达模式的关联；

不同物种之间基因表达模式的分析也是多数文章的研究方向，WGCNA中modulePreservation函数同样也适用于阐述不同物种的基因表达模式的关联，比如人与黑猩猩基因表达模块之间的基因保存性，类似的分析有很多，比如人与牛，人与猪等。

三、modulePreservation 函数实践

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

1. 配置R运行环境

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

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

2. 导入输入数据

2.1 输入表达矩阵

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

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

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

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

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

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

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

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

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

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

2.2 加载模块标签

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

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

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

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

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

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

system.time( {
 mp = modulePreservation(multiExpr, colorList,
             referenceNetworks = c(1:2),
             loadPermutedStatistics = FALSE,
             nPermutations = 200,
             verbose = 3)
} )

save(mp, file = "HumanChimp-HumanSpecific-modulePreservation.RData")

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

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

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

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

我们看主要输出Zsummary评分。

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

一般文献中认为：

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

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

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

参考文献

上一篇下一篇

猜你喜欢

热点阅读