R/Bioconductor与转录组数据分析（6）
2018-12-6 09:58

1、载入需要的包

library(limma)

library(edgeR)

2、加载数据

exprSet <- gse119382.express

rownames(exprSet) <- exprSet[,1]

exprSet <- exprSet[,-1]

3、描述分组

group_list <- c("Drought","Drought","Drought","Water","Water","Water")

4、生成实验设计矩阵：

design <- model.matrix(~0+factor(group_list))

colnames(design)=levels(factor(group_list))

rownames(design)=colnames(exprSet)

design

Drought Water

GSM3373992_D_1       1     0

GSM3373993_D_2       1     0

GSM3373994_D_3       1     0

GSM3373998_W_1       0     1

GSM3373999_W_2       0     1

GSM3374000_W_3       0     1

attr(,"assign")

[1] 1 1

attr(,"contrasts")

attr(,"contrasts")\$factor(group_list)

[1] "contr.treatment"

5、制作差异比较矩阵

contrast.matrix <- makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)

6、根据分组对表达矩阵进行标准化（RNA-seq数据）

exprSet.norm <- voom(exprSet,design,normalize = "quantile",plot = T)

7、用limma进行表达差异分析

fit <- lmFit(exprSet.norm,design)

fit2 <- contrasts.fit(fit,contrast.matrix)

fit2 <- eBayes(fit2)

topTableOutput <- topTable(fit2,coef = 1, n = Inf)

theDEG <- na.omit(topTableOutput)

> results <- decideTests(fit2)

> vennDiagram(results)