差異基因表達(dá)分析

我按照前面的流程轉(zhuǎn)錄組差異表達(dá)分析小實戰(zhàn)（一），將小鼠的4個樣本又重新跑了一遍，從而獲得一個新的count文件：mouse_all_count.txt，有需要的話，可以下載下來進(jìn)行后續(xù)的差異分析。

一般來說，由于普遍認(rèn)為高通量的read count符合泊松分布，所以一些差異分析的R包都是基于負(fù)二項式分布模型的，比如DESeq、DESeq2和edgeR等，所以這3個R包從整體上來說是類似的（但各自標(biāo)準(zhǔn)化算法是不一樣的）。

當(dāng)然還有一個常用的R包則是Limma包，其中的limma-trend和limma-voom都能用來處理RNA-Seq數(shù)據(jù)（但對應(yīng)適用的情況不一樣）

下面準(zhǔn)備適用DESeq2和edgeR兩個R包分別對小鼠的count數(shù)據(jù)進(jìn)行差異表達(dá)分析，然后取兩者的結(jié)果的交集的基因作為差異表達(dá)基因。

1. DEseq2

   library(DESeq2)
   ##數(shù)據(jù)預(yù)處理
   database <- read.table(file = "mouse_all_count.txt", sep = "\t", header = T, row.names = 1)
   database <- round(as.matrix(database))
   
   ##設(shè)置分組信息并構(gòu)建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函數(shù)估計離散度，然后差異分析獲得res對象
   dds <- DESeq(dds)
   res <- results(dds)
   
   #最后設(shè)定閾值，篩選差異基因，導(dǎo)出數(shù)據(jù)(全部數(shù)據(jù)。包括標(biāo)準(zhǔn)化后的count數(shù))
   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個差異基因（篩選標(biāo)準(zhǔn)為padj < 0.05, |log2FoldChange| > 1）

2. edgeR

   library(edgeR)
   ##跟DESeq2一樣，導(dǎo)入數(shù)據(jù)，預(yù)處理（用了cpm函數(shù)）
   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,]
   
   ##設(shè)置分組信息，并做TMM標(biāo)準(zhǔn)化
   exprSet <- DGEList(counts = exprSet, group = group_list)
   exprSet <- calcNormFactors(exprSet)
   
   ##使用qCML（quantile-adjusted conditional maximum likelihood）估計離散度（只針對單因素實驗設(shè)計）
   exprSet <- estimateCommonDisp(exprSet)
   exprSet <- estimateTagwiseDisp(exprSet)
   
   ##尋找差異gene(這里的exactTest函數(shù)還是基于qCML并且只針對單因素實驗設(shè)計)，然后按照閾值進(jìn)行篩選即可
   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個差異基因（篩選標(biāo)準(zhǔn)為FDR < 0.05, |log2FC| > 1）

3. 取DESeq2和edgeR兩者結(jié)果的交集

   diff_gene <- diff_gene_DESeq2[diff_gene_DESeq2 %in% diff_gene_edgeR]

   最終的差異基因數(shù)目為545個

   head(diff_gene)
   [1] "ENSMUSG00000003309.14" "ENSMUSG00000046323.8"  "ENSMUSG00000001123.15"
   [4] "ENSMUSG00000023906.2"  "ENSMUSG00000044279.15" "ENSMUSG00000018569.12"

4. 其他兩個R包（DESeq和limma）就不在這嘗試了，我之前做過對于這4個R包的簡單使用筆記，可以參考下：
   簡單使用DESeq做差異分析
   簡單使用DESeq2/EdgeR做差異分析
   簡單使用limma做差異分析

GO&&KEGG富集分析

以前一直沒有機會用Y叔寫的clusterProfiler包，這次正好看說明用一下。

1. GO富集，加載clusterProfiler包和org.Mm.eg.db包（小鼠嘛），然后將ENSEMBL ID后面的版本號去掉，不然后面不識別這個ID，然后按照clusterProfiler包的教程說明使用函數(shù)即可。

   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)
   ##查看富集結(jié)果數(shù)據(jù)
   enrich_go <- as.data.frame(ego)
   ##作圖
   barplot(ego, showCategory=10)
   dotplot(ego)
   enrichMap(ego)
   plotGOgraph(ego)

2. KEGG富集，首先需要將ENSEMBL ID轉(zhuǎn)化為ENTREZ ID，然后使用ENTREZ ID作為kegg id，從而通過enrichKEGG函數(shù)從online KEGG上抓取信息，并做富集

   library(clusterProfiler)
   library(org.Mm.eg.db)
   
   ##ID轉(zhuǎn)化
   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)
   ##查看富集結(jié)果數(shù)據(jù)
   enrich_kegg <- as.data.frame(kk)
   ##作圖
   dotplot(kk)

到這里為止，轉(zhuǎn)錄組的差異表達(dá)分析算是做完了，簡單的來說，這個過程就是將reads mapping 到reference上，然后計數(shù)獲得count數(shù)，然后做差異分析，最后來個GO KEGG，over了。。。

對于mapping和計數(shù)這兩部還有其實還有好多軟件，具體可見文獻(xiàn)：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis，有時間可以都嘗試下。

至于GO && KEGG這兩步，對于人、小鼠等模式物種來說，不考慮方便因素來說，完全可以自己寫腳本來完成，數(shù)據(jù)可以從gene ontology官網(wǎng)下載，然后就是GO id與gene id相互轉(zhuǎn)化了。KEGG 也是一樣，也可以用腳本去KEGG網(wǎng)站上自行抓取，先知道URL規(guī)則，然后爬數(shù)據(jù)即可。