时间序列研究的是基因表达的动态行为，测量的是一系列和时间点之间有强烈相关性的过程。和针对某一时间点的基因表达进行差异分析不同，时间序列更加关注是发现基因表达的趋势，以有助于理解生物学动态变化过程（比如对刺激的反应、发育过程、周期行为等）。也就是说，时间序列关注的是整体变化趋势而不是某特异表达。这样就有几个挑战，一是要分析的数据量会很大，二是实验条件变多，三是要发掘实验动态变化本质，传统的统计学方法比如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)
Figure 1
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()
Figure 2
网友评论