轉錄組差異表達分析小實戰（二）
轉錄組差異表達分析小實戰（二）
差異基因表達分析
我按照前面的流程轉錄組差異表達分析小實戰（一），將小鼠的4個樣本又重新跑了一遍，從而獲得一個新的count檔案：mouse_all_count.txt，有需要的話，可以下載下來進行後續的差異分析。
一般來說，由於普遍認為高通量的read count符合泊松分佈，所以一些差異分析的R包都是基於負二項式分佈模型的，比如DESeq、DESeq2和edgeR等，所以這3個R包從整體上來說是類似的（但各自標準化演算法是不一樣的）。
當然還有一個常用的R包則是Limma包，其中的limma-trend和limma-voom都能用來處理RNA-Seq資料（但對應適用的情況不一樣）
下面準備適用DESeq2和edgeR兩個R包分別對小鼠的count資料進行差異表達分析，然後取兩者的結果的交集的基因作為差異表達基因。
-
DEseq2
library(DESeq2) ##資料預處理 database <- read.table(file = "mouse_all_count.txt", sep = "\t", header = T, row.names = 1) database <- round(as.matrix(database)) ##設定分組資訊並構建dds物件 condition <- factor(c(rep("control",2),rep("Akap95",2)), levels = c("control", "Akap95")) coldata <- data.frame(row.names = colnames(database), condition) dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition) dds <- dds[ rowSums(counts(dds)) > 1, ] ##使用DESeq函式估計離散度，然後差異分析獲得res物件 dds <- DESeq(dds) res <- results(dds) #最後設定閾值，篩選差異基因，匯出資料(全部資料。包括標準化後的count數) res <- res[order(res$padj),] diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) diff_gene_DESeq2 <- row.names(diff_gene) resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE) write.csv(resdata,file = "control_vs_Akap95.csv",row.names = F)最終獲得572個差異基因（篩選標準為padj < 0.05, |log2FoldChange| > 1）
-
edgeR
library(edgeR) ##跟DESeq2一樣，匯入資料，預處理（用了cpm函式） exprSet<- read.table(file = "mouse_all_count.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE) group_list <- factor(c(rep("control",2),rep("Akap95",2))) exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,] ##設定分組資訊，並做TMM標準化 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)) diff_gene_edgeR <- subset(tTag$table, FDR < 0.05 & (logFC > 1 | logFC < -1)) diff_gene_edgeR <- row.names(diff_gene_edgeR) write.csv(tTag$table,file = "control_vs_Akap95_edgeR.csv")最終獲得688個差異基因（篩選標準為FDR < 0.05, |log2FC| > 1）
-
取DESeq2和edgeR兩者結果的交集
diff_gene <- diff_gene_DESeq2[diff_gene_DESeq2 %in% diff_gene_edgeR]最終的差異基因數目為545個
head(diff_gene) [1] "ENSMUSG00000003309.14" "ENSMUSG00000046323.8" "ENSMUSG00000001123.15" [4] "ENSMUSG00000023906.2" "ENSMUSG00000044279.15" "ENSMUSG00000018569.12" -
其他兩個R包（DESeq和limma）就不在這嘗試了，我之前做過對於這4個R包的簡單使用筆記，可以參考下：
簡單使用DESeq做差異分析
簡單使用DESeq2/EdgeR做差異分析
簡單使用limma做差異分析
GO&&KEGG富集分析
以前一直沒有機會用Y叔寫的clusterProfiler包，這次正好看說明用一下。
-
GO富集，載入clusterProfiler包和org.Mm.eg.db包（小鼠嘛），然後將ENSEMBL ID後面的版本號去掉，不然後面不識別這個ID，然後按照clusterProfiler包的教程說明使用函式即可。
library(clusterProfiler) library(org.Mm.eg.db) ##去除ID的版本號 diff_gene_ENSEMBL <- unlist(lapply(diff_gene, function(x){strsplit(x, "\\.")[[1]][1]})) ##GOid mapping + GO富集 ego <- enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db, keytype = "ENSEMBL", ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05) ##檢視富集結果資料 enrich_go <- as.data.frame(ego) ##作圖 barplot(ego, showCategory=10) dotplot(ego) enrichMap(ego) plotGOgraph(ego) -
KEGG富集，首先需要將ENSEMBL ID轉化為ENTREZ ID，然後使用ENTREZ ID作為kegg id，從而通過enrichKEGG函式從online KEGG上抓取資訊，並做富集
library(clusterProfiler) library(org.Mm.eg.db) ##ID轉化 ids <- bitr(diff_gene_ENSEMBL, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = "org.Mm.eg.db") kk <- enrichKEGG(gene = ids[,2], organism = "mmu", keyType = "kegg", pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.1) ##檢視富集結果資料 enrich_kegg <- as.data.frame(kk) ##作圖 dotplot(kk)
到這裡為止，轉錄組的差異表達分析算是做完了，簡單的來說，這個過程就是將reads mapping 到reference上，然後計數獲得count數，然後做差異分析，最後來個GO KEGG，over了。。。
對於mapping和計數這兩部還有其實還有好多軟體，具體可見文獻：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis，有時間可以都嘗試下。
至於GO && KEGG這兩步，對於人、小鼠等模式物種來說，不考慮方便因素來說，完全可以自己寫指令碼來完成，資料可以從gene ontology官網下載，然後就是GO id與gene id相互轉化了。KEGG 也是一樣，也可以用指令碼去KEGG網站上自行抓取，先知道URL規則，然後爬資料即可。