2021-01-20-单细胞转录组分析实战

参考来源：冻春卷的简书、基迪奥生物

1.准备工作

1.1 项目介绍

A basic overview of Seurat that includes an introduction to:

(1) QC and pre-processing，QC和数据预处理
(2) Dimension reduction，降维
(3) Clustering，聚类
(4) Differential expression，差异分析

1.2 数据来源及下载

本文的范例数据为seurat官网的pbmc-3k数据，可以从官网直接下载：raw data

文件1：barcodes
文件2：ensemble ID to gene name
文件3：真正的测序数据

1.3 R包安装

# 安装Seurat包
install.packages('BiocManager') # 没有的话就安装一个
BiocManager::install('multtest')
install.packages('Seurat')
install.packages("patchwork")
library(dplyr)
library(Seurat)
library(patchwork)

# 加载Seurat包
library(Seurat)

1.4 读入数据

# 读入pbmc数据
pbmc.data <- Read10X(data.dir = "./Seurat_practise/hg19/")

# 查看稀疏矩阵的维度，即基因数和细胞数
dim(pbmc.data)

# 预览稀疏矩阵（1~10行，1~6列）
pbmc.data[1:10,1:6]

2. 创建Seurat对象与数据过滤

2.1 创建Seurat对象并并设置条件筛选细胞

（1）能检测到某个基因的细胞数，即unique基因的分布情况，对应上面的min.cells = 3；
（2）每个细胞能检测到的基因数，对应上面的min.features = 200；
（3）线粒体gene的比例要足够小，使用PercentageFeatureSet函数计算，以MT-开头的则是线粒体gene。

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

# 举例查看数据
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]

在使用CreateSeuratObject()创建对象的同时，过滤数据质量差的细胞。保留在>=3个细胞中表达的基因；保留能够检测到>=200个基因的细胞。

2.2 添加线粒体百分比列

# 计算每个细胞的线粒体基因转录本数的百分比（%），使用[[]]操作符放到metadata中
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc,pattern = "^MT-")

2.3 可视化QC矩阵

把nFeature，nCount，还有MT percentage画在同一个图上，检查数据情况。

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

可以看到，细胞的unique gene（nFeature）和total molecular（nCount）都挺集中的，分别在1000和2000左右，而线粒体基因的比例在百分之2.5左右，合格标准大概设置为多少，需要查一下。

2.4 查看不同特征之间的关系

FeatureScatter顾名思义，就是做Feature的Scatter散点图，这个命令可以用于任何你感兴趣的feature之间的interaction探索。

plot1 <- FeatureScatter(pbmc,feature1="nCount_RNA",
                             feature2="percent.mt")+NoLegend()
plot2 <- FeatureScatter(pbmc,feature1="nCount_RNA",
                             feature2="nFeature_RNA")+NoLegend()
CombinePlots(plots=list(plot1,plot2))

Note:和”nCount_RNA”为细胞的总RNA分子数；nFeature_RNA为细胞的unique基因数。

2.5 设定标准，筛选数据

过滤细胞：保留gene数大于200小于2500的细胞；目的是去掉空GEMs和1个GEMs包含2个以上细胞的数据；而保留线粒体基因的转录本数低于5%的细胞，为了过滤死细胞等低质量的细胞数据

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

3.表达量数据标准化（Normalizing the data）

从数据集中删除不需要的单元格后，下一步是标准化数据。默认情况下，使用全局缩放规范化方法LogNormalize，该方法通过总表达式对每个单元格的特征表达式度量进行标准化，并将其乘以一个缩放因子(默认为10,000)，然后对结果进行log转换。标准化值存储在pbmc[[“RNA”]]@data中。

#表达量数据标准化：LogNormalize的算法：A=log(1+(UMI_A/UMI_Total)*10000)
pbmc <- NormalizeData(pbmc,normalization.method = "LogNormalize",
                           scale.factor = 10000)

# 在默认设定的情况下，上面的的代码可简写为：
pbmc <- NormalizeData(pbmc)

4.鉴定表达高变基因（Feature selection）

使用FindVariableFeatures完成差异分析，选择数据集中差异较高的特征基因（默认2000）并用于下游分析。

# 鉴定表达高变基因（2000个），用于下游分析，如PCA；
pbmc <- FindVariableFeatures(pbmc,selection.method = "vst",
                                  nfeatures = 2000)

# 提取表达量变化最高的10个基因；
top10 <- head(VariableFeatures(pbmc),10)
top10

plot3 <- VariableFeaturePlot(pbmc)+NoLegend()
plot4 <- LabelPoints(plot=plot3,
                     points=top10,
                     repel=TRUE,
                     xnudge = 0,ynudge = 0)
plot4

5.缩放数据（Scaling the data）

PCA分析数据准备，使用ScaleData()进行数据归一化；默认只是标准化高变基因（2000个），速度更快，不影响PCA和分群，但影响热图的绘制。

5.1 缩放的标准

ScaleData函数:

（1）使每个基因在所有细胞间的表达量均值为0
（2）使每个基因在所有细胞间的表达量方差为1

缩放操作给予下游分析同等的权重，这样高表达基因就不会占主导地位，其结果存储在pbmc[[“RNA”]]@scale.data中**

5.2 对所有基因进行缩放

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

缩放是Seurat工作流程中的一个重要步骤，但仅限于将要进行PCA的基因。因此，ScaleData中的默认值只是对前面确定的变量特性(默认值为2,000)执行缩放。要做到这一点，可以忽略前面函数调用中的features参数，即pbmc <- ScaleData(pbmc)。

由于Seurat heatmap(与DoHeatmap一起产生，如下图所示)，需要对heatmap中所有基因进行缩放，以确保高表达的基因不会在heatmap中占主导地位。为了确保我们在后面的热图中不会遗漏任何基因，我们将缩放所有基因。

5.3 删除不需要的变量

pbmc <- ScaleData(pbmc,vars.to.regress = "percent.mt")

使用ScaleData函数从单细胞数据集中删除不需要的变量。例如，我们可以“regress out”与细胞周期阶段或线粒体污染相关的异质性（去除不需要的变量，即校正协变量，校正线粒体基因比例的影响）。

6.线性降维（PCA）

# 线性降维（PCA），默认用高变基因集，但也可通过features参数自己指定；
pbmc <- RunPCA(pbmc,features=VariableFeatures(object=pbmc))

# 检查PCA分群结果，这里只展示前12个PC，每个PC只显示3个基因；
print(pbmc[["pca"]],dims=1:12,nfeatures = 3)

Seurat支持多种方式来可视化PCA 数据：VizDimReduction，DimPlot，DimHeatmap。

# VizDimReduction，小范围试一下
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

# 绘制PCA散点图；
DimPlot(pbmc,reduction = "pca")+NoLegend()

DimHeatmap可以探索数据集的异质性来源，帮助确定用于下游分析的主成分。细胞和features都是根据它们的PCA分数排序。

DimHeatmap(pbmc, dims = 1, 
                 cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, 
                 cells = 500, balanced = TRUE)

可以看到PC1-6对于数据分离比较有帮助，而越往后面帮助越小，可以通过这种方式来确定，下游分析中用于分析的主成分。

7.确定数据集的分群个数

# 方法1：Jackstraw置换检验算法：重复取样（原数据的1%），重跑PCA，鉴定p-value较小的PC；计算’null distribution‘（即零假设成立时）时的基因scores；
pbmc <- JackStraw(pbmc,num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc,dims=1:20)
JackStrawPlot(pbmc,dims=1:15)

# 方法2：肘部图（碎石图），基于每个主成分对方差解释率的排名；
ElbowPlot(pbmc)

分群个数这里选择10，建议尝试选择多个主成分个数做下游分析，对整体影响不大；在选择此参数时，建议选择偏高的数字（“宁缺毋滥”，为了获取更多的稀有分群）；有些亚群很罕见，如果没有先验知识，很难将这种大小的数据集与背景噪声区分开来。

8.细胞聚类（Cluster the cells）

# 基于PCA空间中的欧式距离计算nearest neighbor graph，优化任意两个细胞间的距离权重（输入上一步得到的PC维数）；
pbmc <- FindNeighbors(pbmc, dims = 1:10) # 前10个PC

#接着优化模型，resolution参数决定下游聚类分析得到的分群数，对于3k左右的细胞，设为0.4-1.2能得到较好的结果（官方说明）；如果数据量增大，该参数也应该适当增大；
pbmc <- FindClusters(pbmc, resolution = 0.5) # 分辨率是0.5

# 结果
# Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

# Number of nodes: 2638
# Number of edges: 95930

# 使用Ident()函数可查看不同细胞的分群；
head(Idents(pbmc), 8)

9.非线性降维（UMAP/tSNE）

Seurat提供了几种非线性降维的方法进行数据可视化（在低维空间把相似的细胞聚在一起），比如UMAP和t-SNE，运行UMAP需要先安装‘map-learn’包。

# 安装UMAP
reticulate::py_install(packages ="umap-learn")

# 用umap的方法，并可视化
pbmc <- RunUMAP(pbmc, dims = 1:10) 
umapplot <- DimPlot(pbmc, reduction = "umap", pt.size=1) 
umapplot

用DimPlot()函数绘制散点图，reduction=”tsne”，指定绘制类型；如果不指定，默认先从搜索umap，然后tsne，再然后pca；也可以直接使用这3个函数PCAPlot()、TSNEPlot()、UMAPPlot()；cols、pt.size分别调整分组颜色和点的大小；

# 用tsne的方法，并可视化
pbmc <- RunTSNE(pbmc,dims=1:10)
tsneplot <- TSNEPlot(pbmc,label=TRUE,
                          pt.size=1.5)+NoLegend()
tsneplot

#绘制Marker基因的tsne图；
FeaturePlot(pbmc,features = c("MS4A1","CD14"))

10.寻找亚群的marker基因

10.1 解释标准或参数

Seurat可以找到通过差异表达式定义集群的标记。默认情况下，它识别单个簇的阳性和阴性标记(在ident1中指定)，与所有其他细胞相比较。findallmarker为所有集群自动化这个过程，但是您也可以测试集群组之间的相互关系，或者测试所有细胞。

(1) min.pct参数要求在两组单元中至少检测到一个特征；
(2) thresh.test参数要求一个特性在两组之间有一定的差异(平均)；若两个值都设置为0，会非常耗时——因为这将测试大量不太可能具有高度差异性的features。
(3) 为了加速，可以设置max.cells.per.ident参数，这将向下采样每个标识类，使其单元格数不超过设置的单元格数。虽然通常会有功率的损失，速度的增长可能是显著的，最高度差异表达的特征可能仍然会上升到顶部。

10.2 举例

# 找cluster 1的markers
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1,  # 至少检测到1个feature
                                      min.pct = 0.25) # 设置后会加速，默认是0.1
# 耗时23 s
head(cluster1.markers, n = 5) # 显示！

# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, # 找cluster5的marker
                                      ident.2 = c(0, 3), # 和cluster 0-3对比
                                      min.pct = 0.25)
head(cluster5.markers, n = 5)

# find markers for every cluster compared to all remaining cells, 
# report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, 
                                     min.pct = 0.25, 
                                     logfc.threshold = 0.25)

pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)

10.3 更加精细的寻找marker的方法

我们可以用 test.use 参数设置多种检验方式。例如，用 ROC 返回任何单个标记的“分类能力”。

# 用ROC方法
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, 
                                      logfc.threshold = 0.25, 
                                      test.use = "roc", 
                                      only.pos = TRUE)

10.4 小提琴图可视化marker

画法1

# 小提琴图将marker可视化，随便挑上部分结果中的3个基因
# LDHB, CCR7,LEF1可以将cluster0分别于其他cluster
VlnPlot(pbmc, features = c("LDHB", "CCR7","LEF1")) 

画法2

# raw counts也可以画
VlnPlot(pbmc, features = c("LDHB", "CCR7","LEF1"), 
              slot = "counts", 
              log = TRUE)

画法3

# cluster 图展示
FeaturePlot(pbmc,features = c("LDHB", "CCR7","LEF1"))

画法4

DoHeatmap为给定的细胞和特征生成一个表达式heatmap。在本例中，绘制每个集群的前20个标记(如果小于20，则绘制所有标记)。

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

11.为分群重新指定细胞类型

Fortunately in the case of this dataset, we can use canonical markers to easily match the unbiased clustering to known cell types:

new.cluster.ids <- c("Naive CD4 T","Memory CD4 T","CD14+Mono",
                     "B","CD8 T","FCGR3A+Mono","NK","DC","Platelet")

#将pbmc的水平属性赋值给new.cluster.ids的names属性；
names(new.cluster.ids)
levels(pbmc)
names(new.cluster.ids) <- levels(pbmc)
names(new.cluster.ids)
pbmc <- RenameIdents(pbmc,new.cluster.ids)

#绘制tsne图和UMAP图（修改标签后的）
tsneplot2 <- TSNEPlot(pbmc,label=TRUE,
                           pt.size=0.5)+NoLegend()
tsneplot2

uampplot2 <- DimPlot(pbmc, reduction = "umap", 
                     label = TRUE, 
                     pt.size = 1.2,
                     label.size=5) + NoLegend()
uampplot2

Seurat学习笔记