在進行GEO和TCGA數(shù)據(jù)庫轉(zhuǎn)錄組數(shù)據(jù)挖掘時，差異分析是不可或缺的一部分，一般進行差異分析的主流軟件有三款DEseq2,limma,edgeR,今天小果為大家?guī)淼姆窒硎峭ㄟ^DEseq2進行多組差異分析， 通過該推文將完全掌握利用DEseq2包進行差異分析，值得小伙伴閱讀學習奧！話不多說，和小果一起開啟今天的學習之旅吧！

1.?如何實現(xiàn)DEseq2多組差異分析？

該如何利用DEseq2實現(xiàn)多組差異分析，其實沒那么難，小伙伴只需要準備好基因count矩陣文件和樣本分組信息文件，可以基于分組信息文件進行多組的差異分析，小伙伴們只需要掌握DEseq2 R包參數(shù)使用方法，就可以順利快速的進行分析，小云為大家介紹了是通過批量操作的方式進行多組差異分析，只需要掌握基礎的R語言知識就可以進行自己數(shù)據(jù)的處理，很適合小白奧，那就和小云一起開啟今天的實操吧！

2.?準備需要的R包

DEseq2包直接可以通過Bioconductor 安裝就可以的，非常簡單，小云為小伙伴們附上網(wǎng)址：

#安裝需要的R包

BiocManager::install("DESeq2")

#載入需要的R包

library(DESeq2)

3.?數(shù)據(jù)準備

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)

# 亞型個數(shù)

n.sub <- length(table(subt$TCGA_Subtype))

#創(chuàng)建配對比較的列表信息

group <- subt$TCGA_Subtype

names(group) <- rownames(subt)

# 創(chuàng)建需要配對比較的列表函數(shù),創(chuàng)建了三個分組。

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"])) # 亞型名稱需要根據(jù)情況修改

??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"

??#如果有更多類，按以上規(guī)律繼續(xù)寫

??return(list(sampleList, treatsamList, treatnameList, ctrlnameList))

}

complist <- createList(group=group)

# 配對DESeq2函數(shù)

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)) { # 循環(huán)讀取每一次比較的內(nèi)容

????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,]

????# 差異表達過程，具體參數(shù)細節(jié)及輸出結(jié)果解釋，請參閱相關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"))

????#將分析結(jié)果轉(zhuǎn)化為數(shù)據(jù)框

resData <- as.data.frame(res[order(res$padj),])

????#將行名作為id列

resData$id <- rownames(resData)

????#提取想要的列數(shù)據(jù)

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 = ".", #所有配對差異表達結(jié)果都會輸出在res.path路徑下

???????????????countsTable = expr[,intersect(colnames(expr),rownames(subt))],

???????????????prefix = "SKCM", #文件名以SKCM開頭

???????????????complist = complist,

???????????????overwt = F)

3.?結(jié)果文件

1.?SKCM_deseq2_test_result.immune_vs_Others.txt

該結(jié)果文件為計算的immune與其他分組的差異分析結(jié)果文件,第一列為基因名，第三列為log2FC,第六列為pvalue值，第七列為FDR值(矯正后的pvalue).

?

1.?SKCM_deseq2_test_result.keratin_vs_Others.txt

該結(jié)果文件為該結(jié)果文件為計算的keratin與其他分組的差異分析結(jié)果文件,第一列為基因名，第三列為log2FC,第六列為pvalue值，第七列為FDR值(矯正后的pvalue).

1.?SKCM_deseq2_test_result.MITF-low_vs_Others.txt

該結(jié)果文件為該結(jié)果文件為計算的MITF-low與其他分組的差異分析結(jié)果文件,第一列為基因名，第三列為log2FC,第六列為pvalue值，第七列為FDR值(矯正后的pvalue).

今天小云的分享就到這里啦！如果小伙伴有其他數(shù)據(jù)分析需求，可以嘗試使用本公司新開發(fā)的生信分析小工具云平臺，零代碼完成分析，非常方便奧，云平臺網(wǎng)址為：http://www.biocloudservice.com/home.html,歡迎大家和小云一起討論學習哈?。。?！