一、准备工作
1. 软件包安装
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("pathview")
install.packages("msigdbr")
install.packages("ggnewscale")
library(Seurat)
library(tidyverse)
library(clusterProfiler)
library(pathview)
library(enrichplot)
library(msigdbr)
2. 寻找差异表达基因
读取数据
path_to_pbmc_rds = "sce_final.rds"
pbmc <- readRDS(path_to_pbmc_rds)
筛选差异基因
## 在寻找差异基因之前，把默认的assay切换为RNA。
DefaultAssay(pbmc) <- 'RNA'
## 定义好你想要在哪一个分群基础上找差异表达基因
Idents(pbmc) <- 'seurat_clusters'
## 在不同cluster/或者celltype中找差异表达基因
only.pos = TRUE # 只找上调的差异表达基因
logfc.threshold = 0.25 # 差异基因的avg_log2FC必须要大于0.25
# 对所有Cluster依次进行差异表达分析，将所有Cluster分为两类，我与非我
# all.clusters.markers <- FindAllMarkers(pbmc, only.pos = T, logfc.threshold = logfc.threshold)
# head(all.clusters.markers)
# 只对指定的Cluster进行差异表达分析
markers <- FindMarkers(pbmc, ident.1 = 1, ident.2 = 2, only.pos = T)
head(markers)
markers = markers %>% rownames_to_column('gene') %>% filter(p_val_adj < 0.05)
head(markers)
注：
readRDS()函数读取*.rds对象。.rds与.rdata都是R语言特有的数据保存格式，前者只能保存单一对象，此处用于获取上次教程保存的Seurat对象；后者可以保存某次R运行过程中所有的环境变量，可以下次重复使用，不用重复从头执行代码。DefaultAssay()函数常见的参数有"RNA"和"integrated"。在寻找差异基因之前，把默认的assay切换为RNA，意味着后续的处理将使用原始Count值。设置为integrated意味着使用整合后的数据进行后续操作，此类数据可以认为是经过批处理过的。-
FindAllMarkers()与FindMarkers()函数：-
FindAllMarkers()对所有Cluster依次进行差异表达分析，将所有Cluster分为两类，我与非我。比如对Cluster0进行分析时，其余Cluster1~7被认为是另一个整体。 -
FindMarkers()仅对指定的Cluster进行分析。
-
-
差异表达分析理论：
- FoldChange(FC，姑且翻译为差异倍数)，相同基因在不同样本组（可以是疾病组与正常组，也可以是这里的两个Clusters）中平均表达量差异倍数，是A/B的倍数关系。
- 常见的差异倍数是将上述平均表达量的倍数关系进行取log2（FC）。
- 根据对数图像可以发现，当x趋近于0时，y趋近于-∞，因此看到的最终差异倍数往往是log2[(A+1)/(B+1)]。
- 最终差异表达基因根据提前设定好的logFC和p值阈值即可得到。
差异表达分析结果解读：
- 对Cluster1和2进行差异表达分析，行名是差异表达基因，列名分别为p值，log2FC平均值，pct.1与pct.2表示相应基因在两个Cluster中的表达比例（即多少细胞表达了该基因），矫正p值（是一种更严格的p值，可替代p值作为筛选条件）。
- “markers = markers %>% rownames_to_column('gene') %>% filter(p_val_adj < 0.05)”作用依次为将行名作为列内容添加到dataframe中，相应的数据维度增加一列，该新增列列名为'gene'，最后根据矫正p值小于0.05筛选差异基因。
- “%>%”符号类似于Linux中管道符的作用，将数据从上游传到下游。
Cluster1和2差异表达分析结果
差异表达基因阈值过滤
ID 转换
# 加载物种包
# 根据物种来指定
OrgDb = "org.Hs.eg.db"
# 基因Symbol转化为基因ENTREZID
gene_convert <- bitr(markers$gene, fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = OrgDb)
# 合并两个表格
markers = markers%>%left_join(gene_convert,by=c("gene"="SYMBOL"))
head(markers )
二、富集分析
1. GO富集分析
注意：进行GO分析时需要使用转换后的基因 ENTREZID
# GO
# posiible value: BP, CC, MF, all
ont = "BP"
go.results <- enrichGO(markers$ENTREZID, keyType="ENTREZID",
ont="BP",OrgD = OrgDb, readable = TRUE)
go.results
dotplot(go.results)
barplot(go.results, showCategory = 10)
emapplot(pairwise_termsim(go.results))
GO富集分析可分为三个层次BP, CC, MF，实际使用enrichGO() 函数进行富集分析时可通过指定BP, CC, MF,all四个参数。此处仅使用BP举例说明。
markers = markers%>%left_join(gene_convert,by=c("gene"="SYMBOL")) 作用：
- “gene_convert”数据框中包含基因Symbol和基因ENTREZID信息，markers数据框中包含差异表达基因信息，两者均含有基因但是列名不同。通过管道符及左联函数将两个数据框信息合并，见下图。
三个层次结果图片一起画
参考
org.Hs.eg.db
#GO富集分析
go.results <- enrichGO(gene = markers$ENTREZID,
OrgDb = OrgDb ,
pvalueCutoff =0.05,
qvalueCutoff = 0.05,
ont="all",
readable =T)
write.table(go.results,file="GO.txt",sep="\t",quote=F,row.names = F) #保存富集结果
#柱状图
# tiff(file="barplot.tiff",width = 26,height = 20,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="barplot.pdf",width = 10,height = 8)
barplot(go.results, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
#气泡图
# tiff(file="dotplot.tiff",width = 26,height = 20,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="dotplot.pdf",width = 10,height = 8)
dotplot(go.results,showCategory = 10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
#热图
tiff(file="heatplot.tiff",width = 40,height = 20,units ="cm",compression="lzw",bg="white",res=600)
heatplot(go.results,showCategory =20, foldChange=2,color = "pvalue")
dev.off()
柱状图：这张图可以看出在CC,MF,BP中各个功能的terms富集的基因程度。横坐标是富集在GO term中的基因数左边的是GO的功能，右边是GO属于什么数据库以及可以看出颜色所代表的含义，越红代表越显著.
气泡图：可以到CC,BP,MF下的各个terms富集的基因数量和表达情况。横坐标代表基因所占的比例，右边可以看出点的大小所代表的含义，点越大，富集的基因越多，颜色越红代表富集越显著。
热图：可以看到每个基因在不同的GO的terms的表达程度，这可以选择表达程度最好的几个基因进行分析。
2. KEGG富集分析
KEGG富集分析过程与GO类似
organism = "hsa" # 对应KEGG数据库里的物种缩写
kegg.results <- enrichKEGG(markers$ENTREZID, organism = organism)
write.table(kegg.results,file="KEGG.txt",sep="\t",quote=F,row.names = F) # 保存富集结果
kegg.results <- setReadable(kegg.results, OrgDb = OrgDb, keyType='ENTREZID')
#柱状图
#tiff(file="barplot.tiff",width = 20,height = 12,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="barplot.pdf",width = 10,height = 7)
barplot(kegg.results, drop = TRUE, showCategory = 30)
dev.off()
#气泡图
#tiff(file="dotplot.tiff",width = 20,height = 12,units ="cm",compression="lzw",bg="white",res=600)
pdf(file="bubble.pdf",width = 10,height = 7)
dotplot(kegg.results, showCategory = 30)
dev.off()
#热图
tiff(file="heatplot.tiff",width = 25,height = 15,units ="cm",compression="lzw",bg="white",res=600)
heatplot(kegg.results, showCategory =20, foldChange=cor)
dev.off()
#通路图
library("pathview")
keggxls=read.table("KEGG.txt",sep="\t",header=T)
for(i in keggxls$ID){
pv.out <- pathview(gene.data = cor, pathway.id = i, species = "hsa", out.suffix = "pathview")
}
kegg的通路里有的只是基因的id却不是基因的名字，所以需要转化基因的id，参考。
可视化其中一个通路：
diff_genes_avg_logFC = markers$avg_log2FC
head(diff_genes_avg_logFC)
names(diff_genes_avg_logFC) = markers$ENTREZID
head(diff_genes_avg_logFC)
pathview(gene.data = diff_genes_avg_logFC, species = organism, pathway.id = kegg.results@result$ID[3])
3. GSEA富集分析
GSEA富集分析与GO/KEGG通路富集分析的区别：
- 两者的区别在于Pathway包含基因间互作信息，Gene set只是一个单纯的彼此独立的集合。就像新生第一次组成一个班，彼此不认识，相互没联系，就是Gene set，经过一段时间相处彼此熟络起来，就是Pathway。
- GSEA使用的背景基因称为Gene set，GO/KEGG使用的背景基因称为Pathway。
- 另外，GSEA富集分析不用进行差异表达基因的寻找，但是需要将基因按照log2FC进行排序。
将基因按照log2FC进行排序
# markers for GSEA
## 至少多少比例的细胞表达这个基因，过滤一些只在极少数细胞中有表达的基因
min.pct = 0.01
## 过滤掉在两组中几乎没有差异的基因
logfc.threshold = 0.01
markers.for.gsea <- FindMarkers(pbmc, ident.1 = 1, ident.2 = 2, min.pct = min.pct, logfc.threshold=logfc.threshold)
# GSEA 要求输入的是一个排好序的列表
Markers_genelist <- markers.for.gsea$avg_log2FC
names(Markers_genelist)= rownames(markers.for.gsea)
Markers_genelist <- sort(Markers_genelist, decreasing = T)
虽然说GSEA分析不用进行差异表达基因的寻找，但是为了减少计算量还是提前过滤只在极少数细胞中表达的基因。
- 通过“min.pct = 0.01”和“logfc.threshold = 0.01”参数指定。
- 如果担心过滤掉有意义的基因此处可以将两个参数设为0，此处进行FindMarkers（）操作主要是想得到基因的log2FC用于排序。
进行GSEA富集分析
导入MSigDB
m_df = msigdbr(species = 'Homo sapiens' , category = "C2")
mf_df= m_df %>% dplyr::select(gs_name,gene_symbol)
colnames(mf_df)<-c("term","gene")
gsea.results <- GSEA(Markers_genelist, TERM2GENE = mf_df)
gsea.results
gseaplot(gsea.results, gsea.results@result$ID[3])
# 保存富集结果
write_csv(gsea.results %>% data.frame, "gsea_results.csv")
msigdbr()作用：获取接下来用于注释的MSigDB背景数据集，可以通过“category”参数指定数据集。
GSEA（）函数作用：进行GSEA富集分析，参数“Markers_genelist”中包含按log2FC排序好的基因，参数“TERM2GENE”中包含基因集与基因集中所有的基因信息。
gseaplot（）函数作用：结果图有两个部分，上半部分横坐标为待富集的基因，纵坐标是log2FC，可以看到基因是按log2FC降序排序好的。
下半部分，横坐标反应基因在基因集中命中的情况，如果命中坐标轴上就出现一天竖线，纵坐标表示富集得分（ES,Enrich Score）,命中就+某值，没命中就-某值。
结合上下两部分来看，我们知道此处观察Cluster1与2之间的差异表达情况，当log2FC大于0时，差异基因列表中的基因更多出现在C2参考基因集中，此时FC比值大于1，即基因在Cluster1中更显著表达。[图片上传中...(image-e53601-1679814952079-0)]
可视化“gsea.results@result$ID[3]”该条基因集的富集结果。
结果图解读：
这个图分为三个部分
第一部分：为基因 Enrichment Score 的折线图，横轴为该基因下的每个基因，纵轴为对应的 Running ES，在折线图中有个峰值，该峰值就是这个基因集的 Enrichemnt score，峰值之前的基因就是该基因集下的核心基因，核心基因就是在单细胞很重要很关键的基因，如果基因在顶部富集，我们可以说，从总体上看，该基因集是上调趋势，反之，如果在底部富集，则是下调趋势；
第二部分：也就是hit值，也就是图中的黑线出现的位置，这个黑线可以看出第一部分每次折现点的基因的位置
第三部分：这是rank值图，采用了 Signal2Noise 算法，我对他的理解是，纵轴它反映了基因集中的每个基因的相对表达量，横轴可以看出这是这个基因在基因数据集中出现的位置。
三、GO/KEGG富集分析与GSEA区别
1. 概念
GO：是基因本体联合会所建立的数据库，旨在建立一个适用于各种物种的，对基因和蛋白质功能进行限定和描述的，并能随着研究不断深入而更新的语义词汇标准。GO 提供了一系列的语义用于描述基因功能的概念/类，以及这些概念之间的关系。
KEGG：（京都基因与基因组百科全书）是了解高级功能和生物系统，从分子水平信息，尤其是大型分子数据集生成的基因组测序和其他高通量实验技术的实用程序数据库资源，是国际最常用的生物信息数据库之一，以“理解生物系统的高级功能和实用程序资源库”著称。
GSEA：基因集富集分析，用于确定先验基因集是否在两种生物状态（例如表型）之间差异显著。
2. 区别
GO/KEGG差异基因的一刀切法——仅关注少数几个显著上调或下调的基因，容易遗漏部分差异表达不显著却有重要生物学意义的基因，忽略一些基因的生物特性、基因调控网络之间的关系及基因功能和意义等有价值的信息
GSEA不需要指定明确的差异基因阈值，算法根据实际整体趋势分析。
参考：
https://mp.weixin.qq.com/s/ov-tLoxxK2laBm4NHjUS0Q
https://mp.weixin.qq.com/s/kP7qORDcLNFg4IDoy7m0PQ
https://mp.weixin.qq.com/s/oyj45DLIb-MNAw0LV1-VIA
https://mp.weixin.qq.com/s/pBvhxcTv0crH2VgwJRUN-Q
https://mp.weixin.qq.com/s/RvltqYd4PeECgW6yLoB2Lw
https://mp.weixin.qq.com/s/IowGCKgy-re2tzWDk35lJw
https://mp.weixin.qq.com/s/9jLQFuVATuXcvH-bVvwsRQ
https://mp.weixin.qq.com/s/JhFJJiQL-9Z6uDq4DjsgVA
https://mp.weixin.qq.com/s/QbSgYG_y1wsrMqdzGBRrQg
