表达差异分析
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)
转载本文请联系原作者获取授权，同时请注明本文来自陈新科学网博客。
链接地址：https://wap.sciencenet.cn/blog-62701-1150107.html?mobile=1
收藏
