时间序列研究的是基因表达的动态行为,测量的是一系列和时间点之间有强烈相关性的过程。和针对某一时间点的基因表达进行差异分析不同,时间序列更加关注是发现基因表达的趋势,以有助于理解生物学动态变化过程(比如对刺激的反应、发育过程、周期行为等)。也就是说,时间序列关注的是整体变化趋势而不是某特异表达。这样就有几个挑战,一是要分析的数据量会很大,二是实验条件变多,三是要发掘实验动态变化本质,传统的统计学方法比如t-tests就无能为力了,需要运用新的统计学方法,四是样本间的时间间隔并不总是相等。
maSigPro的全称是Microarray Significant Profiles,采用2步回归策略。第一步选择基因,首先,可以使用统计学程序来鉴定显著表达变化的基因;第二步选择变量,把随时间变化发生显著表达变化的基因进行聚类并且可视化。这个可以通过回归来解决,其中,时间被视为数量变量,实验条件视为分类变量,这样就可以分析趋势变化。
下面简要概述maSigPro分析的步骤、参数及原理。
1. 读入表达数据以及实验设计
data("data.abiotic")
head(data.abiotic)
data("edesign.abiotic")
head(data.abiotic)
data.abiotic为一个芯片表达数据集,每一行为一个基因,每一列为一个样品/条件/重复。maSigPro同样可以处理RNA-seq数据,但是maSigPro没有集合数据标准化的算法,数据读入之前需要首先经过标准化的处理,可以将edge R得到的cpm作为maSigPro的读入数据集。
edesign.abiotic为实验设计表,行名为data.abiotic(读入数据)的列名,Time为时间信息,Replicate为生物学重复信息,生物学重复之间使用相同的数字表示。后面的列比较灵活,用于表示实验设计,有3种常见设计:
1. 考察单一时间变量下基因表达变化(无处理)
所在列均为1即可,如
Time = rep(c(1,5,10,24),each=3)
Replicates = rep(c(1:4),each=3)
Group = rep(1,12)
edesign = cbind(Time,Replicates,Group)
rownames(edesign) = paste("Array",c(1:12), sep= "")
2. 不同处理具有相同的时间起点
data("edesignCT")
head(edesignCT)
不同样品在相同时间点均为1即可。
3. 多种处理组合
在用maSigPro进行分析时,一般情况都是两个或两个以上的感兴趣的变量,其中一个典型的就是时间变量,另外一个通常都是分类变量,代表实验组别(比如不同的处理,细胞株,组织等)。
edesign.abiotic
对应处理列为1,非对应为0。
2. 确定回归模型
d = make.design.matrix(edesign.abiotic, degree = 2)
maSigPro采用多项式回归,参数degree设置多项式使用的次数,在本示例中有3个时间点,degree设置为2,多项式次数设置过高会导致过拟合,一般在能够解释自变量和因变量关系的前提下,次数应该越低越好。
3. 寻找显著基因
这时maSigPro两步法的第一步,计算得到每个基因的回归拟合,并挑选显著基因用于后续分析。
fit = p.vector(data.abiotic, d, Q=0.05, MT.adjust = "BH",counts = FALSE)
Q为FDR阈值,MT.adjust为多重检验方法。counts为布尔值,当为False时表示采用正太分布对芯片数据进行回归拟合,当为True时表示采用负二项分布对Counts数据进行回归拟合。
这一步将筛选出至少在一个group下表现出差异的基因。
4. 寻找显著差异
tstep = T.fit(fit,step.method="backward")
T.fit将采用逐步回归法,查找基因在哪些group下具有显著差异。
5. 获得显著基因集
sigs = get.siggenes(tstep, rsq = 0.7, vars = "groups")
rsq是R-squared的阈值,默认是0.7,作者文中认为当有三个生物学重复时,0.5已经可以得到非常收敛的结果。
vars可以接受3个参数:"groups","all"和"each"。
groups:每一个实验设计groups产生一个列表,包括reference groups以及处理和control之间的差异。
all: 只有一个列表,包含所有的模型变异。
each: 每一种可能的变异为一个列表,包括control和处理之间,以及处理或control和时间尺度之间的差异。
运行masigpro 主要有四步:
5. 结果获取
1. 冷处理条件下,基因表达随时间的变化
data("data.abiotic")
head(data.abiotic)
data.control = data.abiotic[,c(10:18)] ## 提取冷处理数据
data("edesign.abiotic")
head(edesign.abiotic)
edesign.control = edesign.abiotic[c(10:18),c(1:2,4)] ##构建冷处理实验设计
edesign.control
d= make.design.matrix(edesign.control,degree = 2)
fit = p.vector(data.control,d,Q=0.05,counts=FALSE)
fit$i ## 差异基因数
tstep = T.fit(fit)
sigs = get.siggenes(tstep,rsq=0.5,vars="groups")
names(sigs$sig.genes)
pdf("Figure 1.pdf")
cl = see.genes(sigs$sig.genes$Cold,newX11 = F,alfa = 0.05, k = 5)## 对差异结果进行聚类分析
dev.off()
cluster = as.data.frame(cl$cut) ## 提取聚类信息
head(cluster)
2. 多组合处理
d = make.design.matrix(edesign.abiotic,degree = 2)
fit = p.vector(data.abiotic,d,Q=0.05,counts=FALSE)
fit$i
tstep = T.fit(fit)
sigs = get.siggenes(tstep,rsq=0.5,vars="groups")
names(sigs$sig.genes)
pdf("Figure 2.pdf")
cl = see.genes(sigs$sig.genes$ColdvsControl,newX11 = F,alfa = 0.05,k=5)
dev.off()