今天要做的是seurat 的3000 PBMC 的分析。(https://www.bilibili.com/video/BV14q4y1J7MS?spm_id_from=333.337.search-card.all.click)基本流程如下：
2022-01-12 Seurat官网单个样本标准流程

image.png

标准分析 / 鉴定细胞类型的话，有两种方法： 1、Maeker 基因鉴定； 2、SingleR 鉴定

Seurat 核心流程：

data_dir <- "data/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/"
pbmc.counts <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc.counts)
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
pbmc <- FindNeighbors(object = pbmc)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunTSNE(object = pbmc)
pbmc <- RunUMAP(object = pbmc)
DimPlot(object = pbmc, reduction = "tsne")

1、读取数据

这次的数据还是经典的三个文件

image.png

代码：

rm(list=ls())

library(dplyr)
library(Seurat)
library(patchwork)

## =============1.Load the PBMC dataset
data_dir <- "D:/2022-04-01-单细胞/code/data/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19"
pbmc.data <- Read10X(data.dir = data_dir)
pbmc.data

image.png

可以看到在稀疏矩阵里面读取到了32738个基因和2700个细胞。

2、创建 Seurat 对象

创建 Seurat 对象进行一个简单的过滤：
min.cell：每个feature至少在多少个细胞中表达
min.features：每个细胞中至少有多少个feature被检测到

pbmc <- CreateSeuratObject(counts = pbmc.data, 
                           project = "pbmc3k", 
                           min.cells = 3,
                           min.features = 200)
pbmc

image.png

过滤之后可以看到只有 13714 个基因了。
稀疏矩阵就是用 . 代替 0 值，这样更节省运行空间。

3、质控

质控的原因：
1.单个细胞里面
检测到的基因太少，可能是低质量的细胞或者是空的
检测到的基因过多，可能是一个液滴里面包含两个以上的细胞
2.线粒体基因
低质量或者快死的细胞线粒体基因表达会升高

## =============3.QC
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
grep(pattern = "^MT-", rownames(pbmc), value = T)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

grep 查找 pbmc的行名也就是基因中，人线粒体（是以MT-了开头的）基因，查到了13个基因，如下：

image.png

# QC指标
head(pbmc@meta.data, 5)
summary(pbmc@meta.data$nCount_RNA)

# QC指标使用小提琴图可视化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# 密度曲线图
data <- pbmc@meta.data
library(ggplot2)
p <- ggplot(data = data, aes(x=nFeature_RNA)) + geom_density()
p

# 两个指标之间的相关性
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

# 过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

subset 取子集 基因大于200个 同时小于 2500 个 线粒体基因小于 5%

4、数据标准化

## =============4.标准化
# LogNormalize
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

# 标准化后的值保存在：pbmc[["RNA"]]@data
normalized.data <- pbmc[["RNA"]]@data
normalized.data[1:20,1:4]
dim(normalized.data)

5、选择高变基因

高变基因意思是：该基因在有些细胞高表达，但是在其他细胞里面低表达。用到的是均值和方差之间的关系。

## =============5.鉴定高变基因
# 高变基因：在一些细胞中表达高，另一些细胞中表达低的基因
# 变异指标： mean-variance relationship
# 默认返回两千个高变基因，用于下游如PCA降维分析。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 提取前10的高变基因
top10 <- head(VariableFeatures(pbmc), 10)
top10

# 展示高变基因
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

视频里这里展示高变基因可视化的时候是可以的，我这里报错，不知道原因是什么。
横轴是均值，每一个点代表一个基因。纵坐标是方差，表示离散值。红色的就是高变基因。

image.png

6、Scaling the data 归一化数据

这里说一下归一化和标准化的区别：

1、区别介绍
（1）Scaling
归一化用于将数据归一到某一个范围，如[0,1]或者[0,10]，主要是范围上的变化
（2）Normalization
标准化用于改变数据的分布情况，通过将数据分布转变成正态分布
2、适用场景
（1）Scaling
当应用与距离与相似度度量的时候，归一化起到弥足关键的作用，如支持向量机（SVM）和K近邻算法（KNN）
（2）Normalization
当应用要求数据满足正态分布时，对原始数据进行标准化操作，将有利于算法的进行，如线性判别器 (LDA)和高斯朴素贝叶斯（Gaussian naive Bayes）
参考：https://blog.csdn.net/weixin_44704985/article/details/108256543
————————————————

## =============6.Scaling the data
# 归一化处理：每一个基因在所有细胞中的均值变为0，方差标为1，对于降维来说是必需步骤
# 归一化后的值保存在：pbmc[["RNA"]]@scale.data
pbmc <- ScaleData(pbmc)
scale.data <- pbmc[["RNA"]]@scale.data
dim(scale.data)
scale.data[1:10,1:4]


# 可以选择全部基因归一化
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

image.png

归一化处理：每一个基因在所有细胞中的均值变为0，方差标为1，对于降维来说是必需步骤。
也可以用全部基因做归一化，后面的热图会用到。如果你的目标基因不在那2000里面，热图里面就没有。

7、降维

PCA降维，默认使用前面2000个高变基因，可以使用features改变用于降维的基因集。这里还可以加入其他基因集，可以加在object=pbmc 后面。

## =============7.降维
# PCA降维，默认使用前面2000个高变基因，可以使用features改变用于降维的基因集
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# 可视化
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca") 画出来的是下图：表示主成分1和2 ，分别有哪些基因构成，或者是主成分的权重排名。

image.png

DimPlot(pbmc, reduction = "pca") 主成分散点图

image.png

DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)如下图：

image.png

DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE) 也可以画15个，如下图：

image.png

8、细胞聚类

在进行聚类之前，要先确定使用 PC 的个数

## =============8.确定使用PC个数
# each PC essentially representing a ‘metafeature’
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot(pbmc)

JackStrawPlot(pbmc, dims = 1:15) 画出的图如下：要选PC 小于 0.05 的，也就是从 PC13 开始。选择 13个 PC 进入分析。

image.png

ElbowPlot(pbmc) 画出的图如下： 纵坐标是标准差，横坐标是 PC 贡献的累积度。从图可以看出，第一个贡献是最大的，从 PC10 开始，贡献就趋于平滑。官方在这里选择了10 。这里选PC的时候尽可能宽松一些。所以一般按照 p<0.05 选就行。

image.png

## =============9.对细胞聚类
# 首先基于PCA空间构建一个基于欧氏距离的KNN图
pbmc <- FindNeighbors(pbmc, dims = 1:10)

# 聚类并最优化
# resolution参数：值越大，细胞分群数越多，
# 0.4-1.2 typically returns good results for single-cell datasets of around 3K cells
# Optimal resolution often increases for larger datasets. 
pbmc <- FindClusters(pbmc, resolution = 0.5)

# 查看聚类数ID
head(Idents(pbmc), 5)

# 查看每个类别多少个细胞
head(pbmc@meta.data)
table(pbmc@meta.data$seurat_clusters)

这里的 resolution参数：值越大，细胞分群数越多，官方建议：1、0.4-1.2 typically returns good results for single-cell datasets of around 3K cells； 2、Optimal resolution often increases for larger datasets.
3000细胞的时候，选择0.4-1.2。 大的数据集的时候，resolution相应加大。 看到共聚到9类细胞，下图。

image.png
image.png

聚类信息保存在 pbmc@meta 里面了。

9、将细胞在低维空间可视化UMAP/tSNE

## =============10.将细胞在低维空间可视化UMAP/tSNE
pbmc <- RunUMAP(pbmc, dims = 1:10)
pbmc <- RunTSNE(pbmc, dims = 1:10)

# 可视化
DimPlot(pbmc, reduction = "umap", label = T, label.size = 5)
DimPlot(pbmc, reduction = "tsne", label = T, label.size = 5)
saveRDS(pbmc, file = "data/pbmc_tutorial.rds")

DimPlot(pbmc, reduction = "umap", label = T, label.size = 5) 如下图： umap 图细胞和细胞之间隔的远近是有意义的。

image.png

DimPlot(pbmc, reduction = "tsne", label = T, label.size = 5) 如下图： tsne 图没有距离之间的关系。

image.png

10、差异表达分析

## =============11.差异表达分析
# 在cluster2 vs else中差异表达
cluster1.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster1.markers, n = 5)

ident.1 = 2, min.pct = 0.25 就是选择细胞群 2 为一组，其余的细胞群为另一组进行比较。 min.pct 是控制基因在每一类里面的表达比例。

image.png

如果想指定两个类cluster 5 from clusters 0 and 3

cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)

所有细胞群与其他群的差异表达基因
only.pos：只保留上调差异表达的基因

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
write.csv(pbmc.markers,file = "data/pbmc.markers.csv")
head(pbmc.markers)

image.png

pct.1意思是这个基因在 cluster 0 里面的表达比例 pct.2 就是指这个基因在剩余细胞里面的表达比例。

#  筛选
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

每一类筛选 top 2 得到下图：

image.png

# 选择一些基因进行可视化，作者这里根据自己的知识背景选择的相关基因
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features =   c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

可以看出，MS4A1基因基本只在细胞群3中表达。

image.png image.png image.png

11、细胞类型鉴定/注释

这里官网是用经典的Marker基因鉴定的细胞。Marker基因可以来自文献/知识背景，也可以用数据库。例如CellMarker 哈医大做的，还有GESA的C8。

new.cluster.ids <- c("Naive CD4 T",   # IL7R, CCR7
                     "CD14+ Mono",    # CD14, LYZ
                     "Memory CD4 T",  # IL7R, S100A4
                     "B",             # MS4A1
                     "CD8 T",         # CD8A
                     "FCGR3A+ Mono",  # FCGR3A, MS4A7
                     "NK",            # GNLY, NKG7
                     "DC",            # FCER1A, CST3
                     "Platelet")      # PPBP

names(new.cluster.ids) <- levels(pbmc)
new.cluster.ids

# 原来的细胞聚类名称
Idents(pbmc)

# 更改细胞聚类的名字
pbmc <- RenameIdents(pbmc, new.cluster.ids)
Idents(pbmc)
head(pbmc@meta.data)

# 可视化
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 1.2) + NoLegend()
DimPlot(pbmc, reduction = "tsne", label = TRUE, pt.size = 1.2) + NoLegend()

# 保存结果
pbmc@meta.data$cell_anno <- Idents(pbmc)
write.csv(pbmc@meta.data,file = "data/metadata.csv")
saveRDS(pbmc, file = "data/pbmc3k_final.rds")

image.png

接下来介绍其他方法注释细胞 scCATCH包 和 SingleR包

scCATCH包细胞类型比较少

SingleR包

# 清空当前环境变量
rm(list=ls())
options(stringsAsFactors = F)

# 加载包
library(Seurat)
library(SingleR)
#refdata <- HumanPrimaryCellAtlasData()
#save(refdata,file = "data/refdata.RData")
load("data/refdata.RData")
refdata
head(colnames(refdata))
head(rownames(refdata))

# 查看共有多少种细胞类型
unique(refdata@colData@listData[["label.main"]])

# 加载pbmc结果
pbmc <- readRDS("data/pbmc_tutorial.rds")

# 使用的数据为标化后的数据
testdata <- GetAssayData(pbmc, slot="data")
dim(testdata)
testdata[1:30,1:4]

clusters <- pbmc@meta.data$seurat_clusters
table(clusters)

cellpred <- SingleR(test = testdata,  
                    ref = refdata, 
                    labels = refdata$label.main,
                    method = "cluster", 
                    clusters = clusters,
                    assay.type.test = "logcounts", 
                    assay.type.ref = "logcounts")

str(cellpred,max.level = 3)
metadata <- cellpred@metadata
head(metadata)

celltype = data.frame(ClusterID = rownames(cellpred), 
                      celltype = cellpred$labels, 
                      stringsAsFactors = F)
celltype
write.csv(celltype, "data/celltype_anno_SingleR.csv")


# 打分热图上面的注释结果需要校正
p = plotScoreHeatmap(cellpred, clusters = rownames(cellpred), order.by = "cluster")
p

补充 Seurat 对象剖析

疑问1

新给的值会覆盖原来的值 经过这么多赋值过程，原来的pbmc是否被覆盖？

疑问2

细胞聚类基因为什么越来越少？

解答

Seurat 对象相当于是一艘船，里面的各个数据存放在数据槽里面，右图所示。
每个步骤新的赋值，不会保存在原来的地方，而是保存在新的箱子里。所以，原来的值不会被覆盖掉。

image.png

如果要取值（类比取船舱里面的某一个箱子），用 pmbc@ 符号。

取值用@

同时可以参考另一个Up主的视频，也是讲官方流程的。

https://www.jianshu.com/p/98fc6c80c216
https://www.bilibili.com/video/BV1fL411A7JP/?spm_id_from=333.788.recommend_more_video.0