DESeq2 multiple group comparison (biostars.org)
说来也神奇，中文是基本没有关于DESeq2的多组比较是怎么搞的，都说去看文档，搞得像是商业机密一样，但是搜了一下英文是有人做过介绍的,基本运行一下下面的代码就知道怎么搞了老实说不依靠这种，直接弄一个循环也是可以搞的，但是我就是嫌麻烦。。）
# create random data of 12 samples and 4800 genes
x <- round(matrix(rexp(480 * 10, rate=.1), ncol=12), 0)
rownames(x) <- paste("gene", 1:nrow(x))
colnames(x) <- paste("sample", 1:ncol(x))
# fake metadata
coldata <- data.frame(
condition = factor(c(
rep("ctl", 3),
rep("A", 3),
rep("B", 3),
rep("C", 3))))
x
coldata
# set `ctl` as reference level
coldata$condition <- relevel(coldata$condition, ref = "ctl")
str(coldata)
library(DESeq2)
dds <- DESeqDataSetFromMatrix(
countData = x,
colData = coldata,
design= ~ condition)
dds <- DESeq(dds, betaPrior = FALSE)
resultsNames(dds)
res <- results(dds, name="condition_A_vs_ctl")
# same as above but with lfc shrinkage
res <- lfcShrink(dds, coef = 'condition_A_vs_ctl', type = 'apeglm', res = res)
res
results(dds, c("condition", "A", "C"))
此外DESeq2文档里面的确是有案例告诉你怎么搞的，只需要help(results)就有案例了，同样是运行了才知道是怎么搞的,顺便也贴在这里了。
## Example 1: two-group comparison
dds <- makeExampleDESeqDataSet(m=4)
dds <- DESeq(dds)
res <- results(dds, contrast=c("condition","B","A"))
# with more than two groups, the call would look similar, e.g.:
# results(dds, contrast=c("condition","C","A"))
# etc.
## Example 2: two conditions, two genotypes, with an interaction term
dds <- makeExampleDESeqDataSet(n=100,m=12)
dds$genotype <- factor(rep(rep(c("I","II"),each=3),2))
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)
# the condition effect for genotype I (the main effect)
results(dds, contrast=c("condition","B","A"))
# the condition effect for genotype II
# this is, by definition, the main effect *plus* the interaction term
# (the extra condition effect in genotype II compared to genotype I).
results(dds, list( c("condition_B_vs_A","genotypeII.conditionB") ))
# the interaction term, answering: is the condition effect *different* across genotypes?
results(dds, name="genotypeII.conditionB")
## Example 3: two conditions, three genotypes
# ~~~ Using interaction terms ~~~
dds <- makeExampleDESeqDataSet(n=100,m=18)
dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)
# the condition effect for genotype I (the main effect)
results(dds, contrast=c("condition","B","A"))
# the condition effect for genotype III.
# this is the main effect *plus* the interaction term
# (the extra condition effect in genotype III compared to genotype I).
results(dds, contrast=list( c("condition_B_vs_A","genotypeIII.conditionB") ))
# the interaction term for condition effect in genotype III vs genotype I.
# this tests if the condition effect is different in III compared to I
results(dds, name="genotypeIII.conditionB")
# the interaction term for condition effect in genotype III vs genotype II.
# this tests if the condition effect is different in III compared to II
results(dds, contrast=list("genotypeIII.conditionB", "genotypeII.conditionB"))
# Note that a likelihood ratio could be used to test if there are any
# differences in the condition effect between the three genotypes.
# ~~~ Using a grouping variable ~~~
# This is a useful construction when users just want to compare
# specific groups which are combinations of variables.
dds$group <- factor(paste0(dds$genotype, dds$condition))
design(dds) <- ~ group
dds <- DESeq(dds)
resultsNames(dds)
# the condition effect for genotypeIII
results(dds, contrast=c("group", "IIIB", "IIIA"))
用上面的代码写自己的多重比较代码，把整个过程写成函数，只需要简单输入组别，每个组别有几个样本，表达矩阵,比较的结果就全部生成到这个文件夹里面
data<-read.csv("data2.csv")
rownames(data)<-data[,1]
data<-data[,-1]
data#就是普通的表达矩阵
group=c("A","B","C","D") #就是分成多少组
num=c(4,3,4,4) #这几组每组分别有多少个
num
group_list=c("A","B","C","D") #就是分成多少组
num=c(3,3,3,2) #这几组每组分别有多少个
num
group_list
num
Batch_Deseq_differnece<-function(exprSet,group,num,save_dir="Alldiffenece",save_dir2="NEW_MA"){
##create a folder
save_dir<-paste0(save_dir,"/")
dir.create(save_dir)
## creat a group
group_list= factor(rep(group,num))
group_list
colData=data.frame(row.names = colnames(exprSet),
group=group_list)
#dat<-data.frame()
## use the Deseq2 to have Diffence analyse
for (i in 1:length(group)){
name=unique(group)[i]
print(name)
colData$group<-relevel(colData$group,ref=name)
dds=DESeq2::DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~group)
dds <- dds[ rowSums(DESeq2::counts(dds)) > 10, ]
dds <- DESeq2::DESeq(dds)
for (j in 2:length(DESeq2::resultsNames(dds))){
resname=DESeq2::resultsNames(dds)[j]
res=DESeq2::results(dds, name=resname)
res_lfc <- lfcShrink(dds, coef=j, res=res, type="apeglm")
res_lfc
#res=res_lfc
summary(res_lfc)
summary(res)
dir.create(save_dir)
write.csv(res,paste0(save_dir,resname,".csv"))
save_dir2=paste0(save_dir2,"/")
dir.create(save_dir2)
save_dir_MA=paste0(save_dir2,"/",resname)
dir.create(save_dir_MA)
write.csv(res,paste0(save_dir_MA,"/",resname,"_res.csv"))
write.csv(res_lfc,paste0(save_dir_MA,"/",resname,"_reslfc.csv"))
png(paste0(save_dir_MA,"/",resname,"_MA.png"),width=600*3,height=3*600,res=72*3)
plotMA(res, ylim=c(-3,3),main=paste0(resname," MA"))
dev.off()
png(paste0(save_dir_MA,"/",resname,"_MAlfc.png"),width=600*3,height=3*600,res=72*3)
xlim <- c(1,1e5); ylim<-c(-3,3)
plotMA( res_lfc, xlim=xlim, ylim=ylim, main=paste0(resname," apeglm"))
dev.off()
}
}
}
Batch_Deseq_differnece(data2,group=group_list,num,save_dir = "New",save_dir2="NEW_MA")
一个简单的示例
# Load necessary package
if (!requireNamespace("DESeq2", quietly = TRUE)) {
install.packages("DESeq2")
}
library(DESeq2)
# Set seed for reproducibility
set.seed(123)
# Generate a matrix of gene expression data
# Let's assume we have 100 genes and 11 samples
genes <- paste0("gene", 1:100)
samples <- paste0("sample", 1:11)
exprSet <- matrix(sample(1:500, 1100, replace=TRUE), nrow=100, ncol=11,
dimnames=list(genes, samples))
# Define group list and num
group_list <- c("A", "B", "C", "D")
num <- c(3, 3, 3, 2)
# Run the function
Batch_Deseq_differnece(exprSet, group=group_list, num, save_dir = "New", save_dir2="NEW_MA")
网友评论