欢迎关注”生信修炼手册”!
对于转录组的差异分析而言，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代表生物学重复，属于生物学重复的样本其编号相同，后面的列代表样本对应的不同实验条件，采用0和1这种类似二进制的表示法，1表示样本属于这一组，0代表样本不属于这一样，每个实验条件对应一列。
对于只有时间因素的实验，除了Time和Replicate外，只需要再添加一列就可以了，取值全部为1，意味着所有实验条件相同。
上述的样本信息用数据框来保存，示例如下
代码语言：javascript复制> ss.edesign
Time Replicates Group
Array1 1 1 1
Array2 1 1 1
Array3 1 1 1行名称为样本名，列对应上述的分组信息。准备好样本信息后，第一步就是构建回归模型的表达式，代码如下
代码语言：javascript复制design <- make.design.matrix(
sample.group ,
degree = length(unique(sample.group$Time)) - 1)degree代表自由度，取值为时间点减去1。
2. p.vector
回归模型的表达式有了，第一步就是先运用最小二乘法求得每个自变量的系数值，同时，还会用F检验来评估回归方程的显著性，代码如下
代码语言：javascript复制fit <- p.vector(
count,
design,
Q = 0.05,
MT.adjust = "BH",
min.obs = 20)在p.vector函数中，包括以下几个操作步骤
第一个参数count代表基因的表达量矩阵，在运行分析前，默认对基因有一个过滤机制，包括以下两点
- 去除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, ] - 去除在所有样本中表达量都为0的基因，代码如下
# Removing rows with all ceros: sumatot <- apply(dat, 1, sum) counts0 <- which(sumatot == 0) if (length(counts0) > 0) { dat <- dat[-counts0,] }
基因过滤之后，就可以求解回归方程，对于每个基因都进行回归分析，代码如下
代码语言：javascript复制dis <- as.data.frame(design$dis)
model.glm<- glm(y~.,data=dis , family=family, epsilon=epsilon)确定了回归表达式之后，用F检验确定回归方程的显著性，代码如下
代码语言：javascript复制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
上述的回归方程是基于所有的自变量的组合构建的，接下来就是通过逐步回归法确定最佳的自变量组合，代码如下
代码语言：javascript复制tstep <- T.fit(
fit,
step.method = "backward",
alfa = 0.05)逐步回归就是通过在先前建立好的回归方程的基础上，去除其中的某些自变量之后，再次进行回归分析。对于每种自变量的组合，都会有一个对应的回归方程，而且也都会用F检验来评估其显著性。
在挑选最佳的自变量组合时，通过每种自变量组合对应的回归模型的拟合优度值R2来进行判断，R2取值范围为0到1，数值越大，越接近1，回归模型的效果越好。
4. get.siggenes
对于每个基因，根据其自变量的组合，是有对应的多个回归模型的。通过get.siggenes可以查看其中显著性的基因，这个函数有两个关键参数
- rsq rsq指定拟合优度的阈值，如果一个基因的回归模型的拟合优度值小于该阈值，会被过滤掉
- vars
vars的取值有3种，取值为all时每个基因直接给出一个最佳的回归模型，取值为groups时,只给出不同实验条件下相比control组中的差异基因，取值为each时，会给出时间点和实验条件的所有组合对应差异基因列表。
示例用法如下
代码语言：javascript复制sigs <- get.siggenes(
tstep,
rsq = 0.6,
vars = "groups")sigs对象是一个list, 可以通过names(sigs$sig.genes)查看其中包含的组别，然后通过名称来提取其中的差异基因。
对于多个集合的差异基因列表，还可以方便的绘制venn图，代码如下
代码语言：javascript复制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·
—如果喜欢，快分享给你的朋友们吧—
