单细胞-富集分析-GO

开启掘金成长之旅！这是我参与「掘金日新计划 · 12 月更文挑战」的第2天，点击查看活动详情

library(Seurat)
library(ggplot2)
library(Dimplot)
library(tidyverse)
macrophage <- subset(x = scRNA2, subset = labels == "Macrophages")#提出来巨噬细胞亚群,也可以提取不同的多个亚群
#不用跑流程，直接算差异基因，如果是亚群分析，需要重新跑流程

macrophage$celltype.group <- paste(macrophage$labels, macrophage$orig.ident, sep = "_") #添加细胞_组标签
Idents(macrophage) <- "celltype.group"#将新建的标签作为设置标准

mydeg <- FindMarkers(macrophage,ident.1 = 'Macrophages_ill',ident.2 = 'Macrophages_Control', verbose = TRUE, test.use = 'wilcox',min.pct = 0.1)
head(mydeg)#这是疾病与Control之间的差异基因 ident.1是疾病，ident.2是对照，表示疾病与对照相比变化的基因

mydeg_1 <- FindMarkers(macrophage,ident.1 = 'Macrophages_treat',ident.2 = 'Macrophages_ill', verbose = TRUE, test.use = 'wilcox',min.pct = 0.1)
head(mydeg_1)#这是治疗与疾病之间的差异基因

##Control和ill####

#上调的基因
sig_dge.all_ill <- subset(mydeg, p_val_adj<0.05&abs(avg_log2FC)>0.15)#按p校正值<0.05且|abs(log2FC)|>0.15取差异基因，由于10xGenomics单细胞转录组的基因表达较smart-seq2的数据较低，因此一般进行差异分析的时候，其阈值不能直接按照smart-seq2的阈值设置，**seurat一般设为差异倍数的log值为0.25**。
sig_dge.up_ill <- subset(sig_dge.all_ill, p_val_adj<0.05&avg_log2FC>0.15)#疾病状态下相对于对照组上调的基因，默认按照显著性排序
sig_dge.up_ill <- sig_dge.up_ill[order(sig_dge.up_ill$avg_log2FC,decreasing = T),]#按照avg_log2FC降序排列
sig_dge.up_ill_TOP30 <- rownames(sig_dge.up_ill[1:30,])#TOP30

#下调的基因
sig_dge.down_ill <- subset(sig_dge.all_ill, p_val_adj<0.05&avg_log2FC< -0.15)
sig_dge.down_ill <- sig_dge.down_ill[order(sig_dge.down_ill$avg_log2FC,decreasing = T),]
sig_dge.down_ill_TOP30 <- rownames(sig_dge.down_ill[1:30,])


##ill和treat#####

#上调的基因
sig_dge.all_treat <- subset(mydeg_1, p_val_adj<0.05&abs(avg_log2FC)>0.15)
sig_dge.up_treat <- subset(sig_dge.all_treat, p_val_adj<0.05&avg_log2FC>0.15)
sig_dge.up_treat <- sig_dge.up_treat[order(sig_dge.up_treat$avg_log2FC,decreasing = T),]
sig_dge.up_treat_TOP30 <- rownames(sig_dge.up_treat[1:30,])

#下调的基因
sig_dge.down_treat <- subset(sig_dge.all_treat, p_val_adj<0.05&avg_log2FC< -0.15)
sig_dge.down_treat <- sig_dge.down_treat[order(sig_dge.down_treat$avg_log2FC,decreasing = T),]
sig_dge.down_TOP30_treat <- rownames(sig_dge.down_treat[1:30,])

write.table(sig_dge.up_ill,"sig_dge.up_ill.csv",row.names=TRUE,col.names=TRUE,sep=",")
write.table(sig_dge.down_ill,"sig_dge.down_ill.csv",row.names=TRUE,col.names=TRUE,sep=",")
write.table(sig_dge.down_treat,"sig_dge.down_treat.csv",row.names=TRUE,col.names=TRUE,sep=",")
write.table(sig_dge.up_treat,"sig_dge.up_treat.csv",row.names=TRUE,col.names=TRUE,sep=",")


#在自己本地文件中找基因的交集，交集的公式是"=IF (COUNTIF (A:A,G2)>0,G2,0)”来求交集，提取基因名，添加列名为 SYMBOL 
macrophage_genes <- read.csv("macrophage_genes.csv",header=T,col.names = 1)
colnames(macrophage_genes) <- "SYMBOL"
#如果是改多列的名字就用c(),这种命名是按照顺序给各列命名的
colnames(macrophage_genes) <- c(""SYMBOL"","ENTREZID")


#GO分析
BiocManager::install("topGO",force = TRUE)#这个是读取GO数据的R包
BiocManager::install("clusterProfiler",force = TRUE) #clusterProfiler是业界大神Y叔写的R包，专门用于通路富集分析
BiocManager::install("GO.db",force = TRUE)#这个是GO数据集
library(topGO)
library(clusterProfiler)
library(org.Mm.eg.db)
keytypes(org.Mm.eg.db)#查看支持哪些数据类型

#添加一列，内容和列名相同
macrophage_genes<- as.data.frame(macrophage_genes)#把导入的基因文件转化为数据框

##添加ID栏并进行ID转化
ids <- bitr(macrophage_genes$SYMBOL,  #我的的数据框中基因的名字
            fromType = "SYMBOL", #我的ID的数据类型
            toType = "ENTREZID", #转化的数据类型
            OrgDb = org.Mm.eg.db)  
            
##ids	内容像下面写的这个样子
#    SYMBOL   ENTREZID
# 1  Fcna     14133
# 2  Cd5l     11801
# 3  Sdc3     20970
#GO通路富集所有的BP,CC,MF，如果是富集别的内容直接改ont这项。
ego_ALL <- enrichGO(gene          = row.names(sig_dge.all),
                    #universe     = row.names(dge.celltype),
                    OrgDb         = 'org.Mm.eg.db',
                    keyType       = 'SYMBOL',
                    ont           = "ALL",  #设置为ALL时BP, CC, MF都计算
                    pAdjustMethod = "BH",
                    pvalueCutoff  = 0.01,
                    qvalueCutoff  = 0.05)
ego_all <- data.frame(ego_ALL) #后面便于用ggplot进行修改，因为ggplot只接受数据框
write.csv(ego_ALL,'enrichGO_all.csv')#保存GO分析数据
#将柱状图按照P值进行排序，否则长长短短的很丑
ego_bp <- ego_bp[order(ego_bp$p.adjust),] #按照P值排序
ego_bp_top30 <- ego_bp[1 : 30,]#取1-30行
#画图方式1
barplot(ego_bp_top30, showCategory=20,title="EnrichmentGO_BP")#默认的图形颜色，不推荐
#画图方式2
GO_plot<-ggplot(data=ego_bp_top30, aes(x=Description,y=Count)) + 
  geom_bar(stat="identity", width=0.8,fill='salmon1') + 
  coord_flip() +  xlab("Num of Genes") + ylab("GO term") + 
  theme_bw() #推荐
ggsave(file="文件名.png",plot=GO_plot,width=5,height=4)#保存图片
#下面我主要是富集的巨噬细胞BP
macrophage_genes_go_BP <- enrichGO(gene = ids$ENTREZID,
                   OrgDb= org.Mm.eg.db,
                   ont = "BP", ##这里可以更改
                   pAdjustMethod = "BH",
                   minGSSize = 1,
                   pvalueCutoff = 0.01,
                   qvalueCutoff = 0.01,
                   readable = TRUE)