# 转录组差异表达分析小实战（二）

#### 差异基因表达分析

1. 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数)
diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene <- 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）

2. 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）

3. 取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"``````
4. 其他两个R包（DESeq和limma）就不在这尝试了，我之前做过对于这4个R包的简单使用笔记，可以参考下：
简单使用DESeq做差异分析
简单使用DESeq2/EdgeR做差异分析
简单使用limma做差异分析

#### GO&&KEGG富集分析

1. 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)``````
2. 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)``````

RNA-seq analysis，有时间可以都尝试下。