人人都需掌握的单细胞分群聚类分析&Marker基因的可视化

2023-11-14 17:23:02 浏览数 (3)

一年内应该是大家不需要最新版的就是V5版本的SeuratObject 和Seurat，所以都回去吧，回到v4版本哈...

批量转录组测序可以为组织或细胞系的整体基因表达提供信息但无法解析不同细胞间的异质性问题；单细胞测序则可提供每个单个细胞的基因表达特征，允许研究者探索细胞类型、细胞状态、细胞亚群等更细致的分子特征。由于单细胞测序数据具有高度异质性，因此需要特定的数据分析方法来处理和解释数据，如聚类分析、降维算法等。这里意味着单细胞测序分析要比批量转录组的分析更为发杂。

在未成为生信技能树的学徒之前，我从老没有接触过单细胞测序分析，今天有幸和优秀的老师学习，也很庆幸能把我学习到的知识传播给更多的初学者，一起见证我们彼此的成长！好啦，言归正传，今天我们先跟着健明老师学习一下单细胞分群聚类分析&Marker基因的可视化。

以公共数据集GSE182434为例，代码和结果如下：

#安装软件

代码语言：javascript复制
rm(list=ls())
options(stringsAsFactors = F)
getwd()
library(Seurat)
library(ggplot2)
#install.packages(c('clustree'))
library(clustree)
library(cowplot)
library(dplyr)
library(patchwork)

#读入matix数据

代码语言：javascript复制
library(data.table)
dat=fread( "GSE182434_raw_count_matrix.txt.gz",data.table = F)
dim(dat)
dat[1:4,1:4]
rownames(dat)=dat[,1]
dat=dat[,-1]
dat[1:4,1:4]
annotation = fread( "GSE182434_cell_annotation.txt.gz",data.table = F)
annotation[1:4,1:4]
table(annotation$Patient,annotation$CellType)

#筛选相关案例病人样本（DLBCL002 、DLBCL007、 DLBCL008、DLBCL111）

代码语言：javascript复制
#筛选条件①、四个样本 ：DLBCL002 、DLBCL007、 DLBCL008、DLBCL111
#筛选条件②、CD8 细胞
sample = c("DLBCL002","DLBCL007","DLBCL008","DLBCL111")
CellType = "T cells CD8"
annotation = annotation[annotation$Patient %in% sample,]
annotation = annotation[annotation$CellType %in% CellType,]
dim(annotation)  #4个样本，一共5千个CD8细胞
dat = dat[,annotation$ID]
dat[1:5,1:5]
dim(dat)
View(dat)
View(annotation)

#创建Seurat对象

代码语言：javascript复制
library(Seurat)
sce <- CreateSeuratObject(counts = dat, 
                          project = "sce", #项目名称
                          min.cells = 3, 
                          min.features = 200)
sce 
View(sce$nCount_RNA)
View(sce$nFeature_RNA)
View(sce$orig.ident)

sce$orig.ident<-annotation$Patient #修改sce 矩阵样本的分组信息
View(sce$orig.ident)

#识别表达量高变的基因

代码语言：javascript复制
sce = FindVariableFeatures(sce,selection.method = "vst", nfeatures = 2000)#它们在某些细胞中高表达，而在其他细胞中低表达，在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。
#在 Seurat 中的过程，通过直接建模单细胞数据中固有的均值-方差关系来改进以前的版本，并在FindVariableFeatures()函数中实现。默认情况下，为每个数据集返回 2,000 个特征基因用于下游分析。
鉴别前10的高变基因。
top10 <- head(VariableFeatures(sce), 10)
画出没有标签的高变基因图
plot1 <- VariableFeaturePlot(sce)
加入前十个标签
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1   plot2

#简单画一个DimPlot，判断是否要harmony

代码语言：javascript复制
sce = NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)
sce = FindVariableFeatures(sce)
sce = ScaleData(sce, vars.to.regress = c("nFeature_RNA", "percent_mito"))
sce = RunPCA(sce, npcs = 20)#PCA降维，对高变基因降维
print(sce [["pca"]], dims = 1:5, nfeatures = 5)#展示一部分结果
VizDimLoadings(sce , dims = 1:2, reduction = "pca")##点图形式展示
DimPlot(sce, reduction = "pca")#投影的降维图
DimHeatmap(sce, dims = 1:15, cells = 500, balanced = TRUE)#DimHeatmap()允许轻松探索数据集中异质性的主要来源，并且在尝试决定要包括哪些 PC 以进行进一步的下游分析时非常有用。单元格和特征都根据其 PCA 分数排序。

##非线性降维

代码语言：javascript复制
#Seurat 提供了多种非线性降维技术，例如 tSNE 和 UMAP，以可视化和探索这些数据集。这些算法的目标是学习数据的底层流形，以便将相似的单元格放在低维空间中。上面确定的基于图形的集群中的单元格应该在这些降维图上共同定位。作为 UMAP 和 tSNE 的输入，我们建议使用相同的 PC 作为聚类分析的输入。
sce = RunTSNE(sce, npcs = 20)
sce = Runtsne(sce, dims = 1:10)
sce = RunUMAP(sce, dims = 1:10)
sce = FindNeighbors(sce, dims = 1:20, k.param = 60, prune.SNN = 1/15)  #这个dims要和PCA的npcs对应起来。
DimPlot(sce,reduction = "umap",label=T,group.by = "orig.ident")
ggsave(filename="no_harmony_DimPlot.pdf")

#需要 harmony 去批次 因为是主动筛选4个样本出来的

代码语言：javascript复制
library(harmony)
sce <- RunHarmony(sce, "orig.ident")
names(sce@reductions)
sce = RunUMAP(sce,  dims = 1:15, reduction = "harmony")
sce = RunTSNE(sce, npcs = 20,reduction = "harmony")
sce = FindNeighbors(sce, reduction = "harmony", dims = 1:15)
DimPlot(sce,reduction = "umap",label=T,group.by = "orig.ident")
ggsave(filename="harmony_DimPlot.pdf")

##(1)设置分辨率

代码语言：javascript复制
#分辨率分类
#文中是分了6群
#设置不同的分辨率，观察分群效果(选择哪一个？)
for (res in c(0.01, 0.05, 0.1,0.15, 0.2, 0.3, 0.4,0.45,0.5,0.8,0.9,1)) {
  sce=FindClusters(sce, resolution = res, algorithm = 1)
}#细胞分类
colnames(sce@meta.data)
apply(sce@meta.data[,grep("RNA_snn_res",colnames(sce@meta.data))],2,table)
#0.8

(2)可视化高低分辨率的分群情况

代码语言：javascript复制
p1_dim = plot_grid(ncol = 4, DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.01")   
                     ggtitle("louvain_0.01"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.05")   
                     ggtitle("louvain_0.05"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.1")   
                     ggtitle("louvain_0.1"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.2")   
                     ggtitle("louvain_0.2"))
p1_dim

p1_dim = plot_grid(ncol = 5, DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.3")   
                     ggtitle("louvain_0.3"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.5")   
                     ggtitle("louvain_0.5"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.8")   
                     ggtitle("louvain_0.8"), DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.0.9")   
                     ggtitle("louvain_0.9"),DimPlot(sce, reduction = "tsne", group.by = "RNA_snn_res.1")   
                     ggtitle("louvain_1"))
p1_dim

##(3)聚类树

代码语言：javascript复制
p2_tree = clustree(sce@meta.data, prefix = "RNA_snn_res.")
p2_tree
ggsave(plot = p2_tree, filename="Tree_diff_resolution.pdf",width = 10,height = 11)

##(4)分群

代码语言：javascript复制
#最后分群选择的resolution=0.4，细胞聚类为7群。
sce1 <- FindClusters(sce, resolution = 0.4) 
p1 <- DimPlot(sce1, reduction = "umap", label = TRUE, repel = TRUE)
p1 

p2 <- DimPlot(sce1, reduction = "umap", group.by = "orig.ident")
p2

#也可以单独观察每个样本中细胞的聚类情况
p3=DimPlot(sce1, reduction = "umap", split.by = "orig.ident")
p3

#提取各个细胞类型的marker gene

代码语言：javascript复制
#利用 FindMarkers 命令，可以找到找到各个细胞类型中与其他类别的差异表达基因，作为该细胞类型的生物学标记基因。其中ident.1参数设置待分析的细胞类别，min.pct表示该基因表达数目占该类细胞总数的比例
#install.packages("metap")
library(metap)
#1_find all markers of cluster 1
cluster1.markers <- FindMarkers(sce, ident.1 = 6, min.pct = 0.25)
head(cluster1.markers, n = 5)
#利用 DoHeatmap 命令可以可视化marker基因的表达
pbmc.markers <- FindAllMarkers(sce, only.pos = TRUE, min.pct = 0.25)
#2_以cluster 6为例寻找conserved marker
DefaultAssay(sce1) <- "RNA"
nk.markers <- FindConservedMarkers(sce1, ident.1 = 6, grouping.var = "orig.ident", verbose = FALSE)
nk.markers #对每个cluster执行上述操作，就能找出所有细胞类的conserved marker
View(nk.markers)
head(nk.markers, n = 9)
#对marker在不同细胞中的丰度进行可视化
p4=FeaturePlot(sce1, features = c("RPL41", "RPL30", "RPS15A", "RPS27", "RPS12", "RPS28", "RPLP1",
"RPLP10", "RPS26"), min.cutoff = "q9")
ggsave(plot = p4, filename="FeaturePlot_conserved_marker.pdf",width = 10,height = 11)

#探索感兴趣的基因

代码语言：javascript复制
#Seurat提供了许多方法使我们能够方便的探索感兴趣的基因在各个细胞类型中的表达情况
VlnPlot(sce1, features = c("RPL41", "RPL30")) #这里，我们随便选了两个基因
文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

  • 生物信息学马拉松授课（买一得五） ，你的生物信息学入门课
  • 时隔5年，我们的生信技能树VIP学徒继续招生啦
  • 144线程640Gb内存服务器共享一年仍然是仅需800
  • 千呼万唤始出来的独享生物信息学云服务器
  • 生信技能树知识整理实习生又又又开放申请啦
  • 生信共享办公室出租
marker 可视化 数据 算法 线程

0 人点赞