表达差异分析

1、载入需要的包

library(limma)

library(edgeR)

2、加载数据

gse119382.express<- read.table("GSE119382expr.txt",header = T)

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)

按p值排序输出

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

过滤NA值

theDEG <- na.omit(topTableOutput)

8、显示表达差异显著的基因数（adj.p-Value < 0.05）

> results <- decideTests(fit2)

> vennDiagram(results)