2021-04-16 使用GEOquery和limma完成芯片数

今天复现了一下 https://www.jianshu.com/p/1537efae5be9如何使用GEOquery和limma完成芯片数据的差异表达分析。做一下笔记，虽然现在还是不太懂其中的原理，先记着吧哈哈。。。

rm(list=ls())
options(stringsAsFactors = F)
Sys.setenv(LANGUAGE='en')
library(ggplot2)
library(tidyverse)
library(dplyr)
library(GEOmirror)
gset <- geoChina('GSE13535')
save(gset,file='GSE13535.Rdata')
gset <- gset[[1]]
exp <- exprs(gset)
exp1 <- as.data.frame(exp)
pdata=pData(gset)

看一下样本分组信息 sample.png

其中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)
class(treat_time) #character
treat_type <- rep(rep(c('vehicle_control','PE1.3_embolized','PE2.0_embolized'),c(3,4,4)),times=2)
rm(treat_type)
design_df <- data.frame(sample,treat_time,treat_type)

构建单因素试验设计矩阵，进行线性拟合

TS <- paste(design_df$treat_time, design_df$treat_type, sep=".")
TS <- factor(TS,levels = unique(TS))
design <- model.matrix(~0+TS)
fit <- lmFit(exp,design) ##似乎表达矩阵是matrix或者矩阵都行
fit1 <- lmFit(exp1,design)  
class(fit)
# "MArrayLM"
#attr(,"package")
#[1] "limma"

然后根据我们要回答的问题，来设置比较对象。比如不同时间段下控制组哪些基因发生了差异表达，处理18小时后，处理组相对于对照组有哪些基因发生差异表达，也就是做多次t检验。

cont.matrix <- makeContrasts(
  vs1 = TS18h.vehicle_control-TS2h.vehicle_control,#对照组在前后的差异表达基因
  vs2 = TS18h.PE2.0_embolized-TS2h.PE2.0_embolized,#PE2.0处理前后的差异基因
  vs3 = TS18h.PE1.3_embolized-TS2h.PE1.3_embolized,#PE1.3在处理前后差异基因
  # 处理18小时后，PE2.0相对于对照变化的基因再与PE1.3与对照的差异比较
  diff = (TS18h.PE2.0_embolized-TS18h.vehicle_control)-(TS18h.PE1.3_embolized-TS18h.vehicle_control),
  levels = design
)
fit2 <- contrasts.fit(fit,cont.matrix)
results <- decideTests(fit2)
vennDiagram(results)