示例数据
#BiocManager::install('maSigPro') library(maSigPro) # 载入示例数据 data(data.abiotic) data.abiotic[1:5,1:5] data(edesign.abiotic) head(edesign.abiotic) > data.abiotic[1:5,1:5] Control_3H_1 Control_3H_2 Control_3H_3 Control_9H_1 Control_9H_2 STMDF90 0.13735714 -0.3653065 -0.15329448 0.44754535 0.287476796 STMCJ24 NA NA NA NA NA STMJH42 0.07864449 0.1002328 -0.17365488 -0.25279484 0.184855409 STMDE66 0.22976991 0.4740975 0.46930716 0.37101059 -0.004992029 STMIX74 0.14407618 -0.4801864 -0.07847999 0.05692331 0.013045420 > head(edesign.abiotic) Time Replicate Control Cold Heat Salt Control_3H_1 3 1 1 0 0 0 Control_3H_2 3 1 1 0 0 0 Control_3H_3 3 1 1 0 0 0 Control_9H_1 9 2 1 0 0 0 Control_9H_2 9 2 1 0 0 0 Control_9H_3 9 2 1 0 0 0
注意：data.abiotic是已经标准化过的基因表达矩阵。
建立回归模型
生成回归矩阵(makeDesignMatrix)
design <- make.design.matrix(edesign.abiotic, degree = 2) design$groups.vector
示例数据有三个时间的，故考虑二次回归模型(degree = 2)。
> design$groups.vector [1] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control" [5] "ColdvsControl" "HeatvsControl" "SaltvsControl" "Control" [9] "ColdvsControl" "HeatvsControl" "SaltvsControl"
寻找重要基因(p.vector)
F检验确定回归方程的显著性，采用BH的校正方式，校正多重假设检验的p值。
校正后的p值小于p.vector的参数Q的基因就作为候选基因，进行下一步的分析。通过fit$SELEC可以得到候选基因的表达量信息。
fit <- p.vector(data.abiotic, # 标准化的表达矩阵 design, # 实验设计的矩阵 make.design.matrix 生成 Q = 0.05, # 显著性水平 MT.adjust = "BH", min.obs = 20 # 最低表达样本数 不应小于(degree+1)xGroups+1 )
fit$i # 显著性基因的数量 fit$SELEC # 显著性基因表达矩阵
寻找显著性差异(T.fit)
上述的回归方程是基于所有的自变量的组合构建的，接下来就是通过逐步回归法确定最佳的自变量组合。
tstep <- T.fit(fit, # p.vector结果 step.method = "backward", alfa = 0.05) # 在逐步回归中用于变量选择的显著性水平
在挑选最佳的自变量组合时，通过每种自变量组合对应的回归模型的拟合优度值R-squared来进行判断，R-squared取值范围为0到1，数值越大，越接近1，回归模型的效果越好。
获取显著性基因列表(get.siggenes)
sigs <- get.siggenes(tstep, # T.fit结果 rsq = 0.6, # 逐步回归中的R-squared截至值 vars = "groups") # vars参数有3种 # all 每个基因直接给出一个最佳的回归模型 # groups 只给出不同实验条件下相比control组中的差异基因 # each 会给出时间点和实验条件的所有组合对应差异基因列表
names(sigs) names(sigs$sig.genes$ColdvsControl) sigs$sig.genes$ColdvsControl$sig.profiles # 查看cold vs control的结果
结果可视化
韦恩图(suma2Venn)
suma2Venn(sigs$summary[, c(2:4)]) # 左图 suma2Venn(sigs$summary[, c(1:4)]) # 右图 # 这个韦恩图面积大小与数量不成比例 较普通
see.genes()
see.genes(sigs$sig.genes$ColdvsControl, # 差异基因表达矩阵 show.fit = T, # 是否显示回归拟合线(虚线) dis =design$dis, # 回归设计矩阵 cluster.method="hclust" , # 聚类方法 cluster.data = 1, k = 9) # 聚类数目
这一步生成两个图，如图可分别查看。注意调整图片显示区域大小，以免报错。
Cluster Analysis ColdvsControl significant genes
Expression Profiles ColdvsControl significant genes
PlotGroups()
选择某一特定genes的表达进行可视化。
# 选取STMDE66基因 STMDE66 <- data.abiotic[rownames(data.abiotic)=="STMDE66", ] PlotGroups (STMDE66, edesign = edesign.abiotic)
# 添加回归拟合线 PlotGroups (STMDE66, edesign = edesign.abiotic, show.fit = T, dis = design$dis, groups.vector = design$groups.vector)
