本文是參考學(xué)習(xí)單細(xì)胞轉(zhuǎn)錄組基礎(chǔ)分析七：差異基因富集分析
的學(xué)習(xí)筆記?？赡芨鶕?jù)學(xué)習(xí)情況有所改動。

此前的分析我們按轉(zhuǎn)錄特征把細(xì)胞分成了很多類別，例如seurat聚類分析得到的按cluster分類，singleR分析得到的按細(xì)胞類型分類，monocle分析得到的按擬時狀態(tài)（state）分類。不同的細(xì)胞類型之間，有哪些表達(dá)差異基因呢，這些差異基因有特別的意義嗎？

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

library(Seurat)

差異基因GO富集分析

#差異基因GO富集分析

圖片

差異基因kegg富集分析

genelist <- bitr(row.names(sig_dge.celltype), fromType="SYMBOL",

圖片

> library(Seurat)
> library(tidyverse)
> library(patchwork)
> library(monocle)
> library(clusterProfiler)
> library(org.Hs.eg.db)
> rm(list=ls())
> dir.create("enrich")
Warning message:
In dir.create("enrich") : 'enrich' already exists
> scRNA <- readRDS("scRNA4.rds")
> mycds <- readRDS("mycds.rds")
> #比較cluster0和cluster1的差異表達(dá)基因
> dge.cluster <- FindMarkers(scRNA,ident.1 = 0,ident.2 = 1)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s  
> sig_dge.cluster <- subset(dge.cluster, p_val_adj<0.01&abs(avg_logFC)>1)
> #比較B_cell和T_cells的差異表達(dá)基因
> #B T 未定義
> dge.celltype <- FindMarkers(scRNA, ident.1 = 'B_cell', ident.2 = 'T_cells', group.by = 'celltype')
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
> sig_dge.celltype <- subset(dge.celltype, p_val_adj<0.01&abs(avg_logFC)>1)
> #比較擬時State1和State3的差異表達(dá)基因
> p_data <- subset(pData(mycds),select='State')
> scRNAsub <- subset(scRNA, cells=row.names(p_data))
> scRNAsub <- AddMetaData(scRNAsub,p_data,col.name = 'State')
> dge.State <- FindMarkers(scRNAsub, ident.1 = 1, ident.2 = 3, group.by = 'State')
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05s  
> sig_dge.State <- subset(dge.State, p_val_adj<0.01&abs(avg_logFC)>1)
> #差異基因GO富集分析
> ego_ALL <- enrichGO(gene          = row.names(sig_dge.celltype),
+                     #universe     = row.names(dge.celltype),
+                     OrgDb         = 'org.Hs.eg.db',
+                     keyType       = 'SYMBOL',
+                     ont           = "ALL",
+                     pAdjustMethod = "BH",
+                     pvalueCutoff  = 0.01,
+                     qvalueCutoff  = 0.05)
> write.csv(ego_all,'enrich/enrichGO.csv')
Error in is.data.frame(x) : object 'ego_all' not found
> ego_all <- data.frame(ego_ALL)
> write.csv(ego_all,'enrich/enrichGO.csv')
> ego_CC <- enrichGO(gene          = row.names(sig_dge.celltype),
+                    #universe     = row.names(dge.celltype),
+                    OrgDb         = 'org.Hs.eg.db',
+                    keyType       = 'SYMBOL',
+                    ont           = "CC",
+                    pAdjustMethod = "BH",
+                    pvalueCutoff  = 0.01,
+                    qvalueCutoff  = 0.05)
> ego_MF <- enrichGO(gene          = row.names(sig_dge.celltype),
+                    #universe     = row.names(dge.celltype),
+                    OrgDb         = 'org.Hs.eg.db',
+                    keyType       = 'SYMBOL',
+                    ont           = "MF",
+                    pAdjustMethod = "BH",
+                    pvalueCutoff  = 0.01,
+                    qvalueCutoff  = 0.05)
> ego_BP <- enrichGO(gene          = row.names(sig_dge.celltype),
+                    #universe     = row.names(dge.celltype),
+                    OrgDb         = 'org.Hs.eg.db',
+                    keyType       = 'SYMBOL',
+                    ont           = "BP",
+                    pAdjustMethod = "BH",
+                    pvalueCutoff  = 0.01,
+                    qvalueCutoff  = 0.05)
> ego_CC@result$Description <- substring(ego_CC@result$Description,1,70)
> ego_MF@result$Description <- substring(ego_MF@result$Description,1,70)
> ego_BP@result$Description <- substring(ego_BP@result$Description,1,70)
> p_BP <- barplot(ego_BP,showCategory = 10) + ggtitle("barplot for Biological process")
> p_CC <- barplot(ego_CC,showCategory = 10) + ggtitle("barplot for Cellular component")
> p_MF <- barplot(ego_MF,showCategory = 10) + ggtitle("barplot for Molecular function")
> plotc <- p_BP/p_CC/p_MF
> p_BP
> p_CC
> p_MF
> plotc
> ggsave('enrich/enrichGO.png', plotc, width = 12,height = 10)
> genelist <- bitr(row.names(sig_dge.celltype), fromType="SYMBOL",
+                  toType="ENTREZID", OrgDb='org.Hs.eg.db')
'select()' returned 1:1 mapping between keys and columns
Warning message:
In bitr(row.names(sig_dge.celltype), fromType = "SYMBOL", toType = "ENTREZID",  :
  2.33% of input gene IDs are fail to map...
> genelist <- pull(genelist,ENTREZID)
> ekegg <- enrichKEGG(gene = genelist, organism = 'hsa')
Reading KEGG annotation online:

Reading KEGG annotation online:

> p1 <- barplot(ekegg, showCategory=20)
> p2 <- dotplot(ekegg, showCategory=20)
wrong orderBy parameter; set to default `orderBy = "x"`
> plotc = p1/p2
> p1
> p2
> plotc
> ggsave("enrich/enrichKEGG.png", plot = plotc, width = 12, height = 10)
> 
>