生信小白的福音，仅仅几分钟完全掌握DEseq2多组差异分析/SCI论文/科研/研究生/生信分
在进行GEO和TCGA数据库转录组数据挖掘时，差异分析是不可或缺的一部分，一般进行差异分析的主流软件有三款DEseq2,limma,edgeR,今天小果为大家带来的分享是通过DEseq2进行多组差异分析， 通过该推文将完全掌握利用DEseq2包进行差异分析，值得小伙伴阅读学习奥！话不多说，和小果一起开启今天的学习之旅吧！
1. 如何实现DEseq2多组差异分析？
该如何利用DEseq2实现多组差异分析，其实没那么难，小伙伴只需要准备好基因count矩阵文件和样本分组信息文件，可以基于分组信息文件进行多组的差异分析，小伙伴们只需要掌握DEseq2 R包参数使用方法，就可以顺利快速的进行分析，小云为大家介绍了是通过批量操作的方式进行多组差异分析，只需要掌握基础的R语言知识就可以进行自己数据的处理，很适合小白奥，那就和小云一起开启今天的实操吧！
2. 准备需要的R包
DEseq2包直接可以通过Bioconductor 安装就可以的，非常简单，小云为小伙伴们附上网址：
#安装需要的R包
BiocManager::install("DESeq2")
#载入需要的R包
library(DESeq2)
3. 数据准备
input_counts.txt
#基因count矩阵文件,行名为Gene,列为样本名。
input_subtype.txt
#样本分组信息文件，第一列为样本名，第二列为分组信息。
3. DEseq2进行多组差异分析
#读取基因count矩阵文件
expr <- read.table("input_counts.txt",sep = "\t",header = T,check.names = F,stringsAsFactors = F,row.names = 1)
# 读取样本分组信息文件
subt <- read.table("input_subtype.txt", sep = "\t", check.names = F, stringsAsFactors = F, header = T, row.names = 1)
# 亚型名称
n.sub.label <- unique(subt$TCGA_Subtype)
# 亚型个数
n.sub <- length(table(subt$TCGA_Subtype))
#创建配对比较的列表信息
group <- subt$TCGA_Subtype
names(group) <- rownames(subt)
# 创建需要配对比较的列表函数,创建了三个分组。
createList <- function(group=NULL) {
tumorsam <- names(group)
sampleList = list()
treatsamList =list()
treatnameList <- c()
ctrlnameList <- c()
#A-1: 类1 vs 其他
sampleList[[1]] = tumorsam
treatsamList[[1]] = intersect(tumorsam, names(group[group=="immune"])) # 亚型名称需要根据情况修改
treatnameList[1] <- "immune" # 该亚型的命名
ctrlnameList[1] <- "Others" # 其他亚型的命名
#A-2: 类2 vs 其他
sampleList[[2]] = tumorsam
treatsamList[[2]] = intersect(tumorsam, names(group[group=="keratin"]))
treatnameList[2] <- "keratin"
ctrlnameList[2] <- "Others"
#A-3: 类3 vs 其他
sampleList[[3]] = tumorsam
treatsamList[[3]] = intersect(tumorsam, names(group[group=="MITF-low"]))
treatnameList[3] <- "MITF-low"
ctrlnameList[3] <- "Others"
#如果有更多类，按以上规律继续写
return(list(sampleList, treatsamList, treatnameList, ctrlnameList))
}
complist <- createList(group=group)
# 配对DESeq2函数
twoclassDESeq2 <- function(res.path=NULL, countsTable=NULL, prefix=NULL, complist=NULL, overwt=FALSE) {
sampleList <- complist[[1]]
treatsamList <- complist[[2]]
treatnameList <- complist[[3]]
ctrlnameList <- complist[[4]]
allsamples <- colnames(countsTable)
options(warn=1)
for (k in 1:length(sampleList)) { # 循环读取每一次比较的内容
samples <- sampleList[[k]]
treatsam <- treatsamList[[k]]
treatname <- treatnameList[k]
ctrlname <- ctrlnameList[k]
compname <- paste(treatname, "_vs_", ctrlname, sep="") # 生成最终文件名
tmp = rep("others", times=length(allsamples))
names(tmp) <- allsamples
tmp[samples]="control"
tmp[treatsam]="treatment"
outfile <- file.path( res.path, paste(prefix, "_deseq2_test_result.", compname, ".txt", sep="") )
if (file.exists(outfile) & (overwt==FALSE)) { # 因为差异表达分析较慢，因此如果文件存在，在不覆盖的情况下（overwt=F）不再次计算差异表达
cat(k, ":", compname, "exists and skipped;\n")
next
}
saminfo <- data.frame("Type"=tmp[samples],"SampleID"=samples,stringsAsFactors = F)
cts <- countsTable[,samples]
coldata <- saminfo[samples,]
# 差异表达过程，具体参数细节及输出结果解释，请参阅相关document
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = as.formula("~ Type")) # 设计矩阵仅包含亚型信息
dds$Type <- relevel(dds$Type,ref = "control")
dds <- DESeq(dds)
res <- results(dds, contrast=c("Type","treatment","control"))
#将分析结果转化为数据框
resData <- as.data.frame(res[order(res$padj),])
#将行名作为id列
resData$id <- rownames(resData)
#提取想要的列数据
resData <- resData[,c("id","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj")]
#修改列名
colnames(resData) <- c("id","baseMean","log2FC","lfcSE","stat","PValue","FDR")
#输出到文件
write.table(resData, file=outfile, row.names=F, col.names=T, sep="\t", quote=F)
cat(k, ",")
}
options(warn=0)
}
# 差异表达分析过程比较慢请耐心等待
twoclassDESeq2(res.path = ".", #所有配对差异表达结果都会输出在res.path路径下
countsTable = expr[,intersect(colnames(expr),rownames(subt))],
prefix = "SKCM", #文件名以SKCM开头
complist = complist,
overwt = F)
3. 结果文件
1. SKCM_deseq2_test_result.immune_vs_Others.txt
该结果文件为计算的immune与其他分组的差异分析结果文件,第一列为基因名，第三列为log2FC,第六列为pvalue值，第七列为FDR值(矫正后的pvalue).
1. SKCM_deseq2_test_result.keratin_vs_Others.txt
该结果文件为该结果文件为计算的keratin与其他分组的差异分析结果文件,第一列为基因名，第三列为log2FC,第六列为pvalue值，第七列为FDR值(矫正后的pvalue).
1. SKCM_deseq2_test_result.MITF-low_vs_Others.txt
该结果文件为该结果文件为计算的MITF-low与其他分组的差异分析结果文件,第一列为基因名，第三列为log2FC,第六列为pvalue值，第七列为FDR值(矫正后的pvalue).
今天小云的分享就到这里啦！如果小伙伴有其他数据分析需求，可以尝试使用本公司新开发的生信分析小工具云平台，零代码完成分析，非常方便奥，云平台网址为：http://www.biocloudservice.com/home.html,欢迎大家和小云一起讨论学习哈！！！！
