首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >使用maSigPro进行时间序列数据的差异分析

使用maSigPro进行时间序列数据的差异分析

作者头像
生信修炼手册
发布2020-05-08 16:58:18
2.9K0
发布2020-05-08 16:58:18
举报
文章被收录于专栏:生信修炼手册生信修炼手册

欢迎关注”生信修炼手册”!

对于转录组的差异分析而言,case/control的实验设计是最为常见,也最为基础的一种,有很多的R包可以处理这种类型的数据分析。在很多时候,还会有非常复杂的实验设计,比如时间序列, 时间序列与不同实验条件同时存在等情况,对于这种类型的差异分析而言,最常见的分析策略就是回归分析,将基因的表达量看做因变量,将时间和实验条件等因素看自变量,通过回归分析来构建一个合适的模型。

maSigPro是一个用于分析时间序列数据的R包,不仅支持只有时间序列的实验设计,也支持时间序列和分组同时存在的复杂设计,网址如下

https://www.bioconductor.org/packages/release/bioc/html/maSigPro.html

这个R包首先基于多元线性回归模型来拟合时间,实验条件等因素和基因表达量之间的关系,然后运用逐步回归法寻找最佳的自变量组合,具体步骤示意如下

通过5个函数即可实现整个分析流程。

1. makeDesignMatrix

在分析之前,我们需要提供基因的表达量和样本对应的时间序列,实验分组这两种信息。对于表达量而言,需要提供归一化之后的表达量,每一行是一个基因,每一列代表一个样本,这种格式在很多软件中都有介绍,这里就不展开了,对于样本的分组信息,格式如下

sample

Time

Replicate

Control

Case

sample1

1

1

1

0

sample2

1

1

1

0

sample3

1

1

1

0

sample4

2

2

0

1

sample5

2

2

0

1

sample6

2

2

0

1

Time代表样本对应的时间点,Replicate代表生物学重复,属于生物学重复的样本其编号相同,后面的列代表样本对应的不同实验条件,采用01这种类似二进制的表示法,1表示样本属于这一组,0代表样本不属于这一样,每个实验条件对应一列。

对于只有时间因素的实验,除了TimeReplicate外,只需要再添加一列就可以了,取值全部为1,意味着所有实验条件相同。

上述的样本信息用数据框来保存,示例如下

> ss.edesign
       Time Replicates Group
Array1     1          1     1
Array2     1          1     1
Array3     1          1     1

行名称为样本名,列对应上述的分组信息。准备好样本信息后,第一步就是构建回归模型的表达式,代码如下

design <- make.design.matrix(
sample.group ,
degree = length(unique(sample.group$Time)) - 1)

degree代表自由度,取值为时间点减去1。

2. p.vector

回归模型的表达式有了,第一步就是先运用最小二乘法求得每个自变量的系数值,同时,还会用F检验来评估回归方程的显著性,代码如下

fit <- p.vector(
count,
design,
Q = 0.05,
MT.adjust = "BH",
min.obs = 20)

p.vector函数中,包括以下几个操作步骤

第一个参数count代表基因的表达量矩阵,在运行分析前,默认对基因有一个过滤机制,包括以下两点

  1. 去除NA值太多的基因,代码如下 # Removing rows with many missings: count.na <- function(x) (length(x) - length(x[is.na(x)])) dat <- dat[apply(dat, 1, count.na) >= min.obs, ]
  2. 去除在所有样本中表达量都为0的基因,代码如下 # Removing rows with all ceros: sumatot <- apply(dat, 1, sum) counts0 <- which(sumatot == 0) if (length(counts0) > 0) { dat <- dat[-counts0,] }

基因过滤之后,就可以求解回归方程,对于每个基因都进行回归分析,代码如下

dis <- as.data.frame(design$dis)
model.glm<- glm(y~.,data=dis , family=family, epsilon=epsilon)

确定了回归表达式之后,用F检验确定回归方程的显著性,代码如下

model.glm<- glm(y~.,data=dis , family=family, epsilon=epsilon)
model.glm.0<-glm(y~1, family=family, epsilon=epsilon)
test<-anova(model.glm.0,model.glm,test="F")

F检验之后,采用BH的校正方式,校正多重假设检验的p值,代码如下

p.adjusted <- p.adjust(p.vector, method = MT.adjust, n = length(p.vector))

校正后的p值小于p.vector的参数Q的基因就作为候选基因,进行下一步的分析。通过fit$SELEC可以得到候选基因的表达量信息。

3. T.fit

上述的回归方程是基于所有的自变量的组合构建的,接下来就是通过逐步回归法确定最佳的自变量组合,代码如下

tstep <- T.fit(
fit,
step.method = "backward",
alfa = 0.05)

逐步回归就是通过在先前建立好的回归方程的基础上,去除其中的某些自变量之后,再次进行回归分析。对于每种自变量的组合,都会有一个对应的回归方程,而且也都会用F检验来评估其显著性。

在挑选最佳的自变量组合时,通过每种自变量组合对应的回归模型的拟合优度值R2来进行判断,R2取值范围为0到1,数值越大,越接近1,回归模型的效果越好。

4. get.siggenes

对于每个基因,根据其自变量的组合,是有对应的多个回归模型的。通过get.siggenes可以查看其中显著性的基因,这个函数有两个关键参数

  1. rsq rsq指定拟合优度的阈值,如果一个基因的回归模型的拟合优度值小于该阈值,会被过滤掉
  2. vars vars的取值有3种,取值为all时每个基因直接给出一个最佳的回归模型,取值为groups时,只给出不同实验条件下相比control组中的差异基因,取值为each时,会给出时间点和实验条件的所有组合对应差异基因列表。

示例用法如下

sigs <- get.siggenes(
tstep,
rsq = 0.6,
vars = "groups")

sigs对象是一个list, 可以通过names(sigs$sig.genes)查看其中包含的组别,然后通过名称来提取其中的差异基因。

对于多个集合的差异基因列表,还可以方便的绘制venn图,代码如下

suma2Venn(sigs$summary[, c(2:4)])
5. see.genes

对于上述提取到的基因列表,通过see.genes函数可以对其在各个样本中的表达模式进行聚类,代码如下

see.genes(
 sigs$sig.genes$ColdvsControl,
 show.fit = T,
 dis =design$dis,
 cluster.method="hclust" ,
 cluster.data = 1,
 k = 8,
 newX11 = F)

首先是在所有样本中的表达模式的聚类,示意如下

其次是在不同时间点的表达模式,示意如下

maSigPro同时支持芯片和NGS数据的分析,注意表达量必须是归一化之后的表达量。

·end·

—如果喜欢,快分享给你的朋友们吧—

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2018-11-09,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信修炼手册 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1. makeDesignMatrix
  • 2. p.vector
  • 3. T.fit
  • 4. get.siggenes
  • 5. see.genes
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档