例:查看pData
其中title部分告诉了我们分组信息,2小时和18小时,每个时间段又有vehicle control, PE1.3 embolized, PE2.0 embolized,也就是2x2的双因素试验设计, 我们可以现在R语言里构建实验设计的数据框。
sample <- pData$geo_accession
treat_time <- rep(c("2h","18h"),each=11)
treat_type <- rep(rep(c("vehicle_control","PE1.3_embolized","PE2.0_embolized"), c(3,4,4)), times=2)
design_df <- data.frame(sample, treat_time, treat_type)
根据Limma的使用手册的"9.5 Interaction Models: 2 X 2 Factorial Design"进行手续的分析。这里仅仅展示单个因素的分析过程,多个因素看文档依葫芦画瓢就行。
构建单因素试验设计矩阵,进行线性拟合:
TS <- paste(design_df$treat_time, design_df$treat_type, sep=".")
TS <- factor(TS, levels = unique(TS))
后续可有两种方法:
design <- model.matrix(~0+TS)
fit <- lmFit(exprSet, design)
然后设置比较对象。比如不同时间段控制组哪些基因发生了差异表达、处理18小时后处理组相对于对照组有哪些基因发生差异表达,也就是做多次t检验。
cont.matrix <- makeContrasts(
vs1 = TS18.vehicle_control-TS2.vehicle_control,
vs2 = TS18.PE2.0_embolized-TS2.PE2.0_embolized,
vs3 = TS18.PE1.3_embolized-TS2.PE1.3_embolized,
diff = (TS18.PE2.0_embolized-TS18.vehicle_control)-(TS18.PE1.3_embolized-TS18.vehicle_control), levels = design)
# {{vs1:对照组在前后的差异表达基因;
#vs2: PE2.0处理前后的差异基因;
# vs3: PE1.3在处理前后差异基因;
# 处理18小时后,PE2.0相对于对照变化的基因再与PE1.3与对照的差异比较}}
fit2 <- contrasts.fit(fit, cont.matrix)
results <- decideTests(fit2)
design=model.matrix(~factor(TS))
fit=lmFit(exprSet,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
没有实际操作之前,觉得简单的更清爽、适用,但使用后是第一种更可取。差异基因fold change可能搞反!!!在没有显示指定的情况下,我们难以真正确定我们比对的结果是High-Low还是Low-High。另外,前一种方法更利于将差异的比较过程程序化。