墨滴
J

Jay2Echo

V1

2023/02/15阅读：137主题：默认主题

[单细胞|R]单细胞之富集分析

背景

前面介绍了单细胞常规分析中的降维聚类与Cluster注释，现在介绍常见的富集分析。前期内容见下述链接。

[单细胞|R]单细胞降维聚类

[单细胞|R]单细胞自动注释与手动注释

Just Do It！

软件包安装

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("pathview")
install.packages("msigdbr")

library(Seurat)
library(tidyverse)
library(clusterProfiler)
library(pathview)
library(enrichplot)
library(msigdbr)

寻找差异表达基因

path_to_pbmc_rds = "hg19/pbmc3k_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值阈值即可得到。
log2对数函数图像
log2对数函数图像
  • 差异表达分析结果解读：
    • 对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差异表达分析结果
Cluster1和2差异表达分析结果
差异表达基因阈值过滤
差异表达基因阈值过滤

GO富集分析

# 加载物种包
# 根据物种来指定
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"))
# 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))
  • “gene_convert <- bitr(markers$gene, fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = OrgDb)”作用：
    • 基因Symbol转化为基因ENTREZID
  • “markers = markers%>%left_join(gene_convert,by=c("gene"="SYMBOL"))”作用：
    • “gene_convert”数据框中包含基因Symbol和基因ENTREZID信息，markers数据框中包含差异表达基因信息，两者均含有基因但是列名不同。通过管道符及左联函数将两个数据框信息合并，见下图。
基因名转换及合并
基因名转换及合并

  • GO富集分析可分为三个层次BP, CC, MF，实际使用enrichGO（）函数进行富集分析时可通过指定BP, CC, MF,all四个参数。此处仅使用BP举例说明。

    • 气泡图：“dotplot(go.results)”
气泡图
气泡图
  • 柱形图：“barplot(go.results, showCategory = 10)”
柱形图
柱形图

KEGG富集分析

organism = "hsa" # 对应KEGG数据库里的物种缩写
kegg.results <- enrichKEGG(markers$ENTREZID, organism = organism)
kegg.results
kegg.results <- setReadable(kegg.results, OrgDb = OrgDb, keyType='ENTREZID')

dotplot(kegg.results)
barplot(kegg.results, showCategory = 10)
  • KEGG富集分析过程与GO类似。
气泡图
气泡图
柱形图
柱形图

GSEA富集分析

  • GSEA富集分析与GO/KEGG通路富集分析的区别：
    • GSEA使用的背景基因称为Gene set，GO/KEGG使用的背景基因称为Pathway。
      • 两者的区别在于Pathway包含基因间互作信息，Gene set只是一个单纯的彼此独立的集合。
      • 就像新生第一次组成一个班，彼此不认识，相互没联系，就是Gene set，经过一段时间相处彼此熟络起来，就是Pathway。 Pathway与Gene Set的区别
    • 另外，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富集分析

GSEA原理
GSEA原理
MSigDB数据库可选基因集
MSigDB数据库可选基因集 
m_df = msigdbr(species = 'Homo sapiens' , category = "C2")
mf_df= m_df %>% dplyr::select(gs_name,gene_symbol)

gsea.results <- GSEA(Markers_genelist, TERM2GENE = mf_df)
gsea.results

gseaplot(gsea.results, gsea.results@result$ID[3])
  • msigdbr()作用：
    • 获取接下来用于注释的MSigDB背景数据集，可以通过“category”参数指定数据集。
  • GSEA（）函数作用：
    • 进行GSEA富集分析，参数“Markers_genelist”中包含按log2FC排序好的基因，参数“TERM2GENE”中包含基因集与基因集中所有的基因信息。
  • gseaplot（）函数作用：
    • 可视化“gsea.results@result$ID[3]”该条基因集的富集结果。
    • 结果图解读：
      • 结果图有两个部分，上半部分横坐标为待富集的基因，纵坐标是log2FC，可以看到基因是按log2FC降序排序好的。
      • 下半部分，横坐标反应基因在基因集中命中的情况，如果命中坐标轴上就出现一天竖线，纵坐标表示富集得分（ES,Enrich Score）,命中就+某值，没命中就-某值。
      • 结合上下两部分来看，我们知道此处观察Cluster1与2之间的差异表达情况，当log2FC大于0时，差异基因列表中的基因更多出现在C2参考基因集中，此时FC比值大于1，即基因在Cluster1中更显著表达。 GSEA结果图

数据及All_code

数据

数据继承上一条推文，详见本文顶端链接。

All_code

############# 软件安装 ########################

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("pathview")
install.packages("msigdbr")

library(Seurat)
library(tidyverse)
library(clusterProfiler)
library(pathview)
library(enrichplot)
library(msigdbr)

############# 寻找差异表达基因 ##################
## 读取上次课程中经过seurat分析的rds对象，例如：path_to_pbmc_rds="pbmc3k.rds"

path_to_pbmc_rds = "hg19/pbmc3k_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)

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)

############# GO & KEGG 富集分析 ##############
## ID 转换http://127.0.0.1:11035/graphics/plot_zoom_png?width=1920&height=1017
OrgDb = "org.Hs.eg.db" # 根据物种来指定
gene_convert <- bitr(markers$gene, fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = OrgDb)
markers = markers%>%left_join(gene_convert,by=c("gene"="SYMBOL"))

# GO
go.results <- enrichGO(markers$ENTREZID, keyType="ENTREZID",ont="BP",OrgD = OrgDb, readable = TRUE)
go.results
dotplot(go.results)
barplot(go.results, showCategory = 10)


# KEGG
organism = "hsa" # 对应KEGG数据库里的物种缩写
kegg.results <- enrichKEGG(markers$ENTREZID, organism = organism)
kegg.results
kegg.results <- setReadable(kegg.results, OrgDb = OrgDb, keyType='ENTREZID')

dotplot(kegg.results)
barplot(kegg.results, showCategory = 10)
emapplot(pairwise_termsim(kegg.results))

# 可视化其中一个通路
diff_genes_avg_logFC = markers$avg_log2FC
names(diff_genes_avg_logFC) = markers$ENTREZID
pathview(gene.data = diff_genes_avg_logFC, species = organism, pathway.id = kegg.results@result$ID[1])

# 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)

# 导入MSigDB


m_df = msigdbr(species = 'Homo sapiens' , category = "C2")
mf_df= m_df %>% dplyr::select(gs_name,gene_symbol)

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")
  • 推文多平台同步发布，公众号内容食用更佳
  • 更多内容，请关注微信公众号“生信矿工”
  • 如有意见或建议可以在评论区讨论
公众号二维码
公众号二维码

分类：

后端

标签：

后端

作者介绍

J
Jay2Echo
V1