簡單使用DESeq2/EdgeR做差異分析
DESeq2和EdgeR都可用於做基因差異表達分析，主要也是用於RNA-Seq數據，同樣也可以處理類似的ChIP-Seq,shRNA以及質譜數據。
這兩個都屬於R包，其相同點在於都是對count data數據進行處理，都是基於負二項分布模型。因此會發現，用兩者處理同一組數據，最后在相同閾值下篩選出的大部分基因都是一樣的，但是有一部分不同應該是由於其估計離散度的不同方法所導致的。
DESeq2的使用方法：
-
輸入矩陣數據，行名為sample，列名為gene；DESeq2不支持無生物學重復的數據，因此我選擇了2個樣本，3個生物學重復的數據；並對count data取整（經大神指點，這里需要說明下，我的測試數據readcount是RSEM定量的結果，並不是常見的htseq-count的結果，所以count值會有小數點，而DESeq2包不支持count數有小數點，所以這里需要round取整）。
database_all <- read.table(file = "readcount", sep = "\t", header = T, row.names = 1) database <- database_all[,1:6] #type <- factor(c(rep("LC_1",3), rep("LC_2",3))) database <- round(as.matrix(database)) -
設置分組信息以及構建dds對象
condition <- factor(c(rep("LC_1",3), rep("LC_2",3))) coldata <- data.frame(row.names = colnames(database), condition) dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition) -
使用DESeq函數進行估計離散度，然后進行標准的差異表達分析，得到res對象結果
dds <- DESeq(dds) res <- results(dds) -
最后設定閾值，篩選差異基因，導出數據
table(res$padj <0.05) res <- res[order(res$padj),] resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE) write.csv(resdata,file = "LC_1_vs_LC_2.csv")
EdgeR的使用方法：
-
跟DESeq2一樣，EdgeR輸入矩陣數據，行名為sample，列名為gene；DESeq2不支持無生物學重復的數據，因此我選擇了2個樣本，3個生物學重復的數據。
exprSet_all <- read.table(file = "readcount", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE) exprSet <- exprSet_all[,1:6] group_list <- factor(c(rep("LC_1",3), rep("LC_2",3))) -
設置分組信息，去除低表達量的gene以及做TMM標准化
exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,] exprSet <- DGEList(counts = exprSet, group = group_list) exprSet <- calcNormFactors(exprSet) -
使用qCML（quantile-adjusted conditional maximum likelihood）估計離散度（只針對單因素實驗設計）
exprSet <- estimateCommonDisp(exprSet) exprSet <- estimateTagwiseDisp(exprSet) -
尋找差異gene(這里的exactTest函數還是基於qCML並且只針對單因素實驗設計)，然后按照閾值進行篩選即可
et <- exactTest(exprSet) tTag <- topTags(et, n=nrow(exprSet)) tTag <- as.data.frame(tTag) write.csv(tTag,file = "LC_1_vs_LC_2_edgeR.csv")
Summary
以上我主要針對單因素兩兩比較組進行差異分析，其實DESeq2和EdgeR兩個R包都可以對多因素進行差異分析。
-
DESeq2修改以上代碼的分組信息design參數以及在差異分析results函數中添加所選定的分組因素，其他代碼基本一樣，具體參照DESeq2手冊
-
EdgeR則需要用Cox-Reid profile-adjusted likelihood (CR)方法來估算離散度，y <- estimateDisp(y, design)或者分別使用三個函數（y <- estimateGLMCommonDisp(y, design)，y <- estimateGLMTrendedDisp(y, design)， ）y <- estimateGLMTagwiseDisp(y, design)；然后差異表達分析也跟單因素分析不同，主要使用generalized linear model (GLM) likelihood ratio test 或者 quasi-likelihood(QL) F-test，具體代碼可以參照EdgeR手冊。