Seurat_StepV_DEG
此文章是通过学习瑞典国家生物信息学基础设施(NBIS) 所开放的单细胞分析教程加上网上所查找的资料，自身的理解所形成的，可能会有不足之处。该部分是对不同集群之间的差异基因进行分析。
参考来源：https://github.com/NBISweden/workshop-scRNAseq/blob/master/labs/seurat/seurat_05_dge.Rmd
感兴趣的话可以阅读原文。
差异基因
在该节教程中，使用不同的方法来对不同集群之间差异基因，在单细胞分析中，差异基因可以从不同角度进行查找，比比如在同一个样本中的不同细胞clustering之间，或者在不同实验条件样本的细胞clustering之间。除了计算差异基因以外，还将演示合并计算批次效应。
suppressPackageStartupMessages({
library(Seurat)
library(dplyr)
library(cowplot)
library(ggplot2)
library(pheatmap)
library(enrichR)
library(rafalib)
})
alldata <- readRDS("data/results/covid_qc_dr_int_cl.rds")
# 选择reslution为0.5的由louvain算法计算得到的集群信息
sel.clust = "CCA_snn_res.0.5"
alldata <- SetIdent(alldata, value = sel.clust)
table(alldata@active.ident)
#
# 0 1 2 3 4 5 6 7 8 9
# 1499 1021 576 460 422 406 405 254 252 237
# 绘制clustering散点图
plot_grid(ncol = 3, DimPlot(alldata, label = T) + NoAxes(), DimPlot(alldata, group.by = "orig.ident") +
NoAxes(), DimPlot(alldata, group.by = "type") + NoAxes())
细胞标记基因
首先，需要计算每个clustering中高度差异基因的排名，有许多不同参数的不同检验方法可以选择，可以恰当的选择合适自身数据的方法，从而改进结果。我们所寻找的标记基因，是在一种细胞类型中高度表达而在其他细胞类型中不表达的基因。
# 计算差异基因，使用wilcox检验，注意此处使用的还是未矫正批次效应的基因表达矩阵
markers_genes <- FindAllMarkers(alldata, log2FC.threshold = 0.2, test.use = "wilcox",
min.pct = 0.1, min.diff.pct = 0.2, only.pos = TRUE, max.cells.per.ident = 50,
assay = "RNA")
# 选择每个集群中25个最显著的基因进行展示
markers_genes %>%
group_by(cluster) %>%
top_n(-25, p_val_adj) -> top25
top25
展示这25个显著的基因在不同集群中的排名。
# 进行绘图
mypar(2, 5, mar = c(4, 6, 3, 1))
for (i in unique(top25$cluster)) {
barplot(sort(setNames(top25$avg_log2FC, top25$gene)[top25$cluster == i], F),
horiz = T, las = 1, main = paste0(i, " vs. rest"), border = "white", yaxs = "i")
abline(v = c(0, 0.25), lty = c(1, 2))
}
以下将针对每个clustering中排名靠前的差异基因使用不同的方式进行可视化展示。
热图
选择每个集群中最显著的前5个基因进行热图展示。
markers_genes %>%
group_by(cluster) %>%
top_n(-5, p_val_adj) -> top5
# 对所选择的基因的表达矩阵进行标准化
alldata <- ScaleData(alldata, features = as.character(unique(top5$gene)), assay = "RNA")
DoHeatmap(alldata, features = as.character(unique(top5$gene)), group.by = sel.clust,
assay = "RNA")
点图
另一种展示差异基因的方法是在点图中展现所有集群中的基因表达即表达率。
DotPlot(alldata, features = rev(as.character(unique(top5$gene))), group.by = sel.clust,
assay = "RNA") + coord_flip()
小提琴图
# 选择集群中排名中前3的基因
top5 %>%
group_by(cluster) %>%
top_n(-3, p_val) -> top3
# pt.size展示的是图中散点的大小
VlnPlot(alldata, features = as.character(unique(top3$gene)), ncol = 5, group.by = sel.clust, assay = "RNA", pt.size = 0)
上面的差异基因分析皆是建立在差异基因分析方法为wilcox秩和检验的基础上的，但差异基因的检验方法除了wilcox秩和检验外，还有“bimod”（Likelihood-ratio test）、“roc”（使用 ROC 分析识别标记基因））、“t”（ t 检验）、“negbinom”（负二项式广义线性模型）、“poisson”（泊松广义线性模型）、“LR”（逻辑回归）、“MAST”（hurdle model)） ），“DESeq2”（负二项分布）等等。
不同实验条件之间的差异分析
计算差异基因的第二种方法就是，在不同样本中的同一种细胞类型的基因，他们之间差异表达的基因都有什么。比如，在这个案例中，提供了一个来自患者和健康人群的文库，探索这不同条件下，在特定细胞类型中哪些基因受到了比较大的影响。
# 选择集群1中的所有细胞
cell_selection <- subset(alldata, cells = colnames(alldata)[alldata@meta.data[, sel.clust] ==
2])
cell_selection <- SetIdent(cell_selection, value = "type")
# Compute differentiall expression
DGE_cell_selection <- FindAllMarkers(cell_selection, log2FC.threshold = 0.2, test.use = "wilcox",
min.pct = 0.1, min.diff.pct = 0.2, only.pos = TRUE, max.cells.per.ident = 50,
assay = "RNA")
绘制不同类型样本之间的基因表达情况。
DGE_cell_selection %>%
group_by(cluster) %>%
top_n(-5, p_val) -> top5_cell_selection
VlnPlot(cell_selection, features = as.character(unique(top5_cell_selection$gene)),
ncol = 5, group.by = "type", assay = "RNA", pt.size = 0.1)
也可以绘制这些基因在其他clustering中的表达量，依旧是在不同样本之间的表达，用来检验这些基因在其他细胞类型是否同样也上/下调。
VlnPlot(alldata, features = as.character(unique(top5_cell_selection$gene)), ncol = 5,
split.by = "type", assay = "RNA", pt.size = 0)
基因富集分析
在得到了差异基因集之后，就可以使用超几何检验对这些差异基因进行基因富集分析，从而分析他们的功能差异。
一般用来分析的基因集有以下几种：
GO_Biological_Process_2017bKEGG_2019_HumanKEGG_2019_MouseWikiPathways_2019_HumanWikiPathways_2019_Mouse
# 加载基因富集分析R包
library(enrichR)
# 检查能够进行基因富集分析的基因集数量
enrichR::listEnrichrDbs()
# 进行基因富集分析
enrich_results <- enrichr(genes = DGE_cell_selection$gene[DGE_cell_selection$cluster ==
"Covid"], databases = "GO_Biological_Process_2017b")[[1]]
# 可视化富集结果
par(mfrow = c(1, 1), mar = c(3, 25, 2, 1))
barplot(height = -log10(enrich_results$P.value)[10:1], names.arg = enrich_results$Term[10:1],
horiz = TRUE, las = 1, border = FALSE, cex.names = 0.6)
abline(v = c(-log10(0.05)), lty = 2)
abline(v = 0, lty = 1)
基因集富集分析（GSEA）
除了超几何分布检验进行基因富集分析之外，还可以进行基因集富集分析Gene Set Enrichment Analysis，在差异基因集中对基因进行排序与评分（通常是基于变化倍数FC）并进行检验计算，查看是否特定的基因集是否较多的存在与上下调基因，或者无调节基因。
DGE_cell_selection <- FindMarkers(cell_selection, ident.1 = "Covid", log2FC.threshold = -Inf,
test.use = "wilcox", min.pct = 0.1, min.diff.pct = 0, only.pos = FALSE, max.cells.per.ident = 50,
assay = "RNA")
# 在基因集中基于变化倍数创建一个基因排序列表
gene_rank <- setNames(DGE_cell_selection$avg_log2FC, casefold(rownames(DGE_cell_selection),
upper = T))
决定完基因排序之后么，就可以使用从MSigDB中获取到的特征基因集，对获取到的基因集进行富集分析，这里选择KEGG富集通路展示。
# install.packages('msigdbr')
library(msigdbr)
# 下载基因集
msigdbgmt <- msigdbr::msigdbr("Homo sapiens")
msigdbgmt <- as.data.frame(msigdbgmt)
# 列出可用的基因集
unique(msigdbgmt$gs_subcat)
# [1] "MIR:MIR_Legacy" "TFT:TFT_Legacy" "CGP" "TFT:GTRD"
# [5] "" "VAX" "CP:BIOCARTA" "CGN"
# [9] "GO:BP" "GO:CC" "IMMUNESIGDB" "GO:MF"
# [13] "HPO" "CP:KEGG" "MIR:MIRDB" "CM"
# [17] "CP" "CP:PID" "CP:REACTOME" "CP:WIKIPATHWAYS"
# 选择想要用的基因子集
msigdbgmt_subset <- msigdbgmt[msigdbgmt$gs_subcat == "CP:WIKIPATHWAYS", ]
gmt <- lapply(unique(msigdbgmt_subset$gs_name), function(x) {
msigdbgmt_subset[msigdbgmt_subset$gs_name == x, "gene_symbol"]
})
names(gmt) <- unique(paste0(msigdbgmt_subset$gs_name, "_", msigdbgmt_subset$gs_exact_source))
接下来就可以使用GSEA，这将生成一个包含多个通路信息的表格，然后，我们可以对这些通路信息进行排序和过滤，仅保留最显著的通路。可以通过p值或标准化的富集分数normalized enrichemnet score (NES).进行过滤。
library(fgsea)
# 进行富集分析
fgseaRes <- fgsea(pathways = gmt, stats = gene_rank, minSize = 15, maxSize = 500)
fgseaRes <- fgseaRes[order(fgseaRes$NES, decreasing = T), ]
# 过滤基因集，过滤结果中只包含前10上调或下调的通路信息
top10_UP <- fgseaRes$pathway[1:10]
# 表格绘制
dev.off()
plotGseaTable(gmt[top10_UP], gene_rank, fgseaRes, gseaParam = 0.5)
保存数据
saveRDS(alldata, "data/3pbmc_qc_dr_int_cl_dge.rds")
write.csv(markers_genes)
扫一扫
1177
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
