3大差异分析r包:DESeq2、edgeR和limma
致敬说明书
2.多组比较
生信星球解释较清楚,多组搭配
library(limma)##加载包
load("GSE37761exp_groupfile-00.RData")#加载匹配过的基因表达矩阵
eset=exp ##将基因表达矩阵赋值给eset
targets<-read.csv("group_file.csv",row.names = 1)##读入样本数据,包括两列,
#第一列GSM号,第二列为样本分组#该数据集中实际样本分组中存在24h_1、24h_5、24h_21等,然而这不符合R的命名原则
targets$Target=gsub("_",".",targets$Target)##将"_"替换成“.”,也可以不替换
##该数据集中实际样本分组中存在24h_1、24h_5、24h_21等,然而这不符合R的命名原则,所以在没个分类前加一个“F”,具体自己定
targets$Target=c(paste0("F",c(targets$Target),collapse = NULL,sep=""))
colnames(targets)=c("FileName","Target")#更改列名,为了和limma包中的一致
lev<-unique(targets$Target)##使用unique()函数进行去重
f <- factor(targets$Target, levels=lev)
design <- model.matrix(~0+f) #样本矩阵
colnames(design) <- lev #更改列名为levels名
###两两之间的比较,求差异基因使用topTable函数
cont.wt <- makeContrasts("F12h.1-Fcontrol",#举例
levels=design)
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, cont.wt)
fit2 <- eBayes(fit2)
tT=topTable(fit2, adjust="BH",sort.by="logFC",n=Inf)
tT = subset(tT, select=c("adj.P.Val","P.Value","logFC"))
colnames(tT)=c("FDR","P.Value","logFC")
###
##后面可以设置fdr及logFC的阈值进行筛选,进而得到差异基因列表,此处略。
##两组以上的比较,寻找差异基因,使用topTableF函数
cont.wt <- makeContrasts("F12h.1-Fcontrol",##举例
"F24h.1-F12h.1",
"F5d.1-F24h.1",
"F10d.1-F5d.1",
levels=design)
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, cont.wt)
fit2 <- eBayes(fit2)
tT=topTableF(fit2, adjust="BH",sort.by="F",n=Inf)
##tT = subset(tT, select=c("adj.P.Val","P.Value","logFC")) 这一条和下一条代码需要根据自己需要进行更改,topTableF函数得到的结果中包括F值,而logFC则是没两组的比较,均得出了
##colnames(tT)=c("FDR","P.Value","logFC")
setwd("E://Rworkplace")
data<-read.table("GSE98793.txt")
group_list<-read.table("gene_list.txt")
new_group<- group_list[order(group_list[,1]),]
group<-new_group[,2]
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group))
colnames(design)=levels(factor(group))
rownames(design)=colnames(data)
contrast.matrix<-makeContrasts("HC-MDD",levels=design) # 这里好像顺序反了欸
##step1
fit <- lmFit(data,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
write.table(nrDEG,"limma_GSE98793.txt")