单细胞测序—基础分析流程

单细胞测序—基础分析流程

注：基于官方文档

参考官方文档：https://satijalab.org/seurat/articles/pbmc3k_tutorial

数据：https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz

1 读取数据与创建Seurat对象

10X的输入数据是固定的三个文件，在工作目录下新建01_data/，把三个文件放进去。

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

pbmc.data <- Read10X(data.dir = "01_data/")
#即32738个基因，2700个细胞
dim(pbmc.data)
[1] 32738  2700

pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2  3 . . . . . 3 4 1 5
TCL1A . .  . . . . . . 1 . . . . . . . . . . .  . 1 . . . . . . . .
MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .

range(pbmc.data)
[1]   0 419
#一个基因至少要在3个细胞里面有表达，才被保留；
#一个细胞里面至少要表达两百个基因，才被保留。
pbmc <- CreateSeuratObject(counts = pbmc.data, 
                           project = "pbmc3k", 
                           min.cells = 3, 
                           min.features = 200)
pbmc

注：

10X的输入数据是固定的三个文件，这三个文件分别存储了什么信息？

• barcode.tsv

这个文件包含了细胞条形码（cell barcodes），每一行对应一个细胞。条形码是用于区分不同细胞的唯一标识符。通常，每个条形码由一串字母和数字组成。

• genes.tsv（新版数据格式中为features.tsv））

这个文件包含了基因的信息，每一行对应一个基因。通常包含两列数据：

第一列是基因的唯一标识符（如Ensembl ID）。

第二列是基因的常用名称（如“CD3D”）。在新版数据格式中，可能还有第三列标识特征的类型（如Gene Expression, Antibody Capture等）。

• matrix.mtx

这个文件是一个稀疏矩阵文件，存储了每个基因在每个细胞中的计数数据。矩阵的行对应genes.tsv中的基因，列对应barcodes.tsv中的细胞。文件内容包括：

行数（基因数）。

列数（细胞数）。

非零元素的数量。

具体的计数值（基因在细胞中的表达量），以三元组形式存储：行索引、列索引和计数值。

这些文件结合起来，提供了每个细胞的基因表达信息，通常用于后续的单细胞RNA测序数据分析。

稀疏矩阵

矩阵中的 . 值表示 0（未检测到分子）。由于 scRNA-seq 矩阵中的大多数值为 0，因此 Seurat 尽可能使用稀疏矩阵表示。这为 Drop-seq/inDrop/10x 数据节省了大量的内存和速度。

>pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2  3 . . . . . 3 4 1 5
TCL1A . .  . . . . . . 1 . . . . . . . . . . .  . 1 . . . . . . . .
MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .

CreateSeuratObject：这个函数从原始计数数据创建一个Seurat对象。参数包括：

counts：计数矩阵。

project：表示项目名称的字符串。

min.cells：基因在至少多少个细胞中表达才被包括。

min.features：每个细胞中检测到的最少基因数

pbmc：这个变量存储创建的Seurat对象，其中包括元数据和标准化数据等。

2 质控

这里是对细胞进行的质控，指标是：

• 线粒体基因含量不能过高
• nFeature_RNA 不能过高或过低

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
head(pbmc@meta.data)

VlnPlot(pbmc, 
        features = c("nFeature_RNA",
                     "nCount_RNA", 
                     "percent.mt"), 
        ncol = 3,pt.size = 0.5)

根据这个三个图，确定了这个数据的过滤标准：

nFeature_RNA在200~2500之间；线粒体基因占比在5%以下。

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

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

注：

相关代码解释

head(pbmc@meta.data)：

• 这个命令显示了Seurat对象pbmc中meta.data的前几行。meta.data包含了每个细胞的元数据，例如每个细胞的基因和UMI计数等。

percent.mt

计算线粒体基因的比例，并将其存储在percent.mt中。PercentageFeatureSet函数的pattern参数用于匹配基因的名字，这里使用正则表达式^MT-来匹配所有以“MT-”开头的基因，这些基因通常代表线粒体基因。再次查看meta.data，现在可以看到多了一个percent.mt的列。

ncount与nfeature辩析

nFeature_RNA是每个细胞中检测到的基因数量。

nCount_RNA是细胞内检测到的分子总数(UMI)。

nFeature_RNA过低，表示该细胞可能已死/将死或是空液滴。太高的nCount_RNA和/或nFeature_RNA表明“细胞”实际上可以是两个或多个细胞。结合线粒体基因count数除去异常值，即可除去大多数双峰/死细胞/空液滴，因此它们过滤是常见的预处理步骤。

FeatureScatter函数绘制两个散点图，plot1展示了nCount_RNA与percent.mt的关系，plot2展示了nCount_RNA与nFeature_RNA的关系，得到前者的相关系数为-0.13，后者的相关系数为0.95。这说明了什么，为什么要做这一步?

nCount_RNA与percent.mt的关系（相关系数为-0.13）：

• 这个相关系数接近于0，表明nCount_RNA和percent.mt之间几乎没有线性关系。这意味着线粒体基因的比例在不同细胞中与总的RNA计数之间没有明显的关联。
• 但如果观察到明显的负相关（相关系数为负且绝对值较大），可能意味着细胞存在线粒体基因的异常高表达（如细胞凋亡）。因此，这一步有助于识别和过滤掉可能有问题的细胞。

nCount_RNA与nFeature_RNA的关系（相关系数为0.95）：

• 这个相关系数非常接近1，表明nCount_RNA和nFeature_RNA之间有很强的正相关关系。这是符合预期的，因为一般来说，检测到的基因数量（nFeature_RNA）与总的RNA计数（nCount_RNA）应该成正比。
• 强正相关关系表明数据质量良好，数据没有受到过多的噪音或双重细胞的干扰。

3 标准化与寻找高变基因HVG（特征选择）

计算在数据集中表现出高细胞间差异的特征子集（即它们在某些细胞中表达高，而在其他细胞中表达低）。下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。

#消除测序深度的影响
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
top10 <- head(VariableFeatures(pbmc), 10);top10

#这里选了2000个，把前十个在图上标记出来
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, 
                     points = top10, 
                     repel = TRUE)
plot1 + plot2

注：

NormalizeData

这个步骤使用NormalizeData函数对数据进行标准化。标准化是为了消除不同细胞之间测序深度的差异，从而使不同细胞之间的表达水平可以进行比较。通常，标准化会将每个细胞中的基因表达值除以该细胞中的总表达量，然后乘以一个标量（如1e4），最后取对数转化。

pbmc <- FindVariableFeatures(pbmc)

top10 <- head(VariableFeatures(pbmc), 10);top10

使用FindVariableFeatures函数识别数据集中变异性较高的基因。这些基因在下游分析中（如聚类和降维）起到重要作用，因为它们能更好地区分不同的细胞类型或状态。

提取并显示了变异性最高的前10个基因。这些基因是根据变异度排序的，可以用于进一步的分析和注释。

plot1/plot2

VariableFeaturePlot函数绘制高变异基因的散点图，展示基因的平均表达水平（平均表达值）与变异程度（标准差）的关系。

LabelPoints函数用于在图中标注特定的基因，这里是标注前10个高变异基因。repel = TRUE参数表示避免标签重叠，使图更加清晰。

4 缩放数据与降维

4.1 PCA

pbmc <- ScaleData(pbmc, features = rownames(pbmc))
pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
#每个主成分对应基因的热图
DimHeatmap(pbmc, dims = 1:15, cells = 500)
# 应该选多少个主成分进行后续分析
ElbowPlot(pbmc)
# 限速步骤
f = "jc.Rdata"
if(!file.exists(f)){
  pbmc <- JackStraw(pbmc, num.replicate = 100)
  pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
  save(pbmc,file = f)
}
load(f)
JackStrawPlot(pbmc, dims = 1:20)

注：

ScaleData(pbmc, features = rownames(pbmc))

ScaleData函数对数据进行中心化和标准化。中心化是指减去平均值，标准化是将数据除以标准差。这一步使得每个基因在所有细胞中的表达值具有相同的量纲，防止高表达基因对下游分析的影响。这里features = rownames(pbmc)表示对所有基因进行缩放。

问：在之前的代码中，执行过pbmc <- FindVariableFeatures(pbmc)，我的理解是高变基因集已经被复制为pbmc,这种理解是错误的吗?

答：理解混淆，FindVariableFeatures函数的确是用来识别数据集中变异性较高的基因集，但它不会修改或复制原始数据集pbmc。相反，它会在pbmc对象的内部存储这些高变异基因的信息，以供后续分析使用。具体来说，FindVariableFeatures函数会计算每个基因的变异度，并将高变异的基因记录在pbmc对象的一个叫做VariableFeatures的属性中。这个属性包含了经过筛选后被认为在不同细胞中具有显著变异性的基因列表。

因此，执行pbmc <- FindVariableFeatures(pbmc)之后：

1. pbmc对象本身没有被复制或更改：它仍然包含所有的原始数据和元数据。
2. 高变异基因集（高变基因）被标记：这些基因的信息被存储在pbmc对象中，但这并不意味着其他基因被移除或对象被复制。

pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))

RunPCA函数对数据进行主成分分析（PCA），主要用于降维和特征提取。分析的特征是之前识别的高变异基因（VariableFeatures(pbmc)）。PCA帮助识别数据中变化最大的方向，并将这些方向作为新的坐标轴（主成分），减少数据的维度。

VizDimLoadings

VizDimLoadings函数可视化前两个主成分（PC1和PC2）上基因的加载值。加载值代表每个基因在主成分上的贡献大小，帮助识别哪些基因在特定主成分上有较大的影响。

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

DimHeatmap函数绘制前15个主成分的热图，展示每个主成分上的前500个细胞和基因。这可以帮助识别每个主成分上哪些基因对其有较大的贡献。

ElbowPlot(pbmc)

ElbowPlot函数绘制碎石图（Elbow plot），展示每个主成分的标准差。这帮助选择用于后续分析的主成分数目。图中通常会出现一个"肘部"，即标准差开始显著下降的点，选择这个点之前的主成分数目通常是合适的。

重要性：选取合适数量的主成分可以避免过拟合，同时保留足够的生物学信息用于下游分析。

JackStraw分析：

JackStraw函数使用置换试验（permutation test）来评估每个主成分的显著性。这里num.replicate = 100表示进行100次置换。

打分：ScoreJackStraw函数计算每个主成分的显著性分数。

保存和加载结果：如果文件jc.Rdata不存在，则运行JackStraw分析并保存结果；否则加载现有结果。

PCA聚类

#PCAPlot(pbmc) + NoLegend()
DimPlot(pbmc, reduction = "pca")+ NoLegend()
# 结合JackStrawPlot和ElbowPlot，挑选10个PC，所以这里dims定义为1:10
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5) #分辨率
# 结果聚成几类，用Idents查看
length(levels(Idents(pbmc)))
[1] 9

注：

DimPlot(pbmc, reduction = "pca") + NoLegend()

PCAPlot 和 DimPlot 函数用于可视化数据在前两个主成分上的分布。PCAPlot是Seurat v2版本的函数，而DimPlot是Seurat v3及更高版本的函数，后者功能更强大，可以选择不同的降维方法（如PCA、UMAP、t-SNE等）。NoLegend() 函数用于去除图例，使得图形更简洁。这些图显示了不同细胞在前两个主成分上的分布情况，有助于识别数据中是否存在明显的聚类。

pbmc <- FindNeighbors(pbmc, dims = 1:10)

结合JackStrawPlot和ElbowPlot，用户可以选择具有统计显著性的主成分数目。在本例中，用户选择了前10个主成分（dims = 1:10）用于后续分析。这意味着在接下来的步骤中，数据的主要变异性将由这10个主成分来表示。FindNeighbors 函数根据之前选择的主成分，构建每个细胞的K近邻图（K-nearest neighbor graph）。这个图表示每个细胞与其他细胞的邻近关系，是聚类分析的基础。

pbmc <- FindClusters(pbmc, resolution = 0.5)

FindClusters 函数基于K近邻图对细胞进行聚类。resolution参数控制聚类的分辨率，较高的分辨率会产生更多的簇（更小的聚类），较低的分辨率会产生更少的簇（更大的聚类）。在这里，resolution = 0.5 是一个常见的选择，用于生成适中的聚类数目。

4.2 UMAP/t-SNE

使用UMAP（Uniform Manifold Approximation and Projection）算法对单细胞RNA测序数据进行降维

pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")

注：

pbmc <- RunUMAP(pbmc, dims = 1:10)

RunUMAP 函数基于指定的主成分（这里是前10个主成分）运行UMAP降维算法。UMAP是一种非线性降维方法，旨在将高维数据映射到低维空间（通常是二维或三维）中，同时保留数据的全局和局部结构。

UMAP常用于单细胞RNA测序数据的可视化，因为它能够有效地展示数据中的簇状结构（即不同的细胞群体）。

DimPlot(pbmc, reduction = "umap")

DimPlot 函数用于可视化降维结果。这里指定reduction = "umap"，表示使用UMAP降维结果进行绘图。这个图展示了每个细胞在UMAP空间中的位置，不同的颜色通常代表不同的聚类结果（即不同的细胞群体）。UMAP图可以帮助研究者直观地观察数据中的细胞群体，并识别不同细胞类型或状态。

UMAP图的目的是以一种易于理解和解释的方式展示数据中的复杂结构。相比于PCA，UMAP更适合用于展示数据中的非线性关系和复杂结构，尤其是在高维数据中。在单细胞RNA测序数据分析中，UMAP和t-SNE（t-distributed Stochastic Neighbor Embedding）是常用的降维和可视化方法。它们的目的是将数据中的高维特征压缩到2D或3D空间中，以便识别和解释数据中的簇或模式。

问：执行UMAP是否还有执行PCA的必要呢？单细胞测序的后续分析流程，是否是主要基于UMAP的分析结果呢？

答：执行UMAP之前仍然有必要先执行PCA。

原因如下：

1. PCA作为初步降维步骤

• 降噪和加速计算：PCA是线性降维方法，可以将高维数据投射到一个较低维度的空间，通常选取具有最大变异性的前几百个主成分。这有助于减少数据的噪声，并加速后续的非线性降维算法如UMAP和t-SNE的计算。
• 降维和数据压缩：PCA可以将大部分信息浓缩到少数几个主成分中，有效降低数据的复杂度。因此，使用PCA后提取的主成分作为UMAP输入，有助于减少计算负担，同时保留数据的主要结构。

1. UMAP的独特功能和优势

• 非线性降维：UMAP是一种非线性降维技术，能够更好地保留数据中的复杂和非线性关系。这使得它非常适合用于高维数据的可视化。
• 局部和全局结构：UMAP比PCA更适合展示数据的局部和全局结构，能够清晰地展示数据中的簇状结构（即不同的细胞群体）。

1. 单细胞测序数据分析流程中的UMAP和PCA

• PCA作为预处理步骤：尽管UMAP可以直接应用于原始数据，但通常先进行PCA以减少数据的维度和噪声，选择PCA提取的主成分作为UMAP的输入。
• UMAP用于可视化和聚类：在单细胞RNA测序数据分析中，UMAP图常用于展示细胞群体的分布和聚类结果。UMAP可以帮助识别不同的细胞类型或状态，因此经常用于数据的最终可视化步骤。

问：umap是基于PCA的结果执行，为什么在代码中没有看出来？

答：UMAP并不一定是必须基于PCA的结果执行的，但在实践中，常常会先进行PCA降维，然后再进行UMAP。这是因为PCA能够有效地降低数据的维度和噪声，从而加速UMAP的计算并提高结果的稳定性。

pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))这一步执行了PCA，并将结果存储在pbmc对象中;pbmc <- RunUMAP(pbmc, dims = 1:10)这一步运行UMAP降维算法，其中dims = 1:10表示使用PCA结果的前10个主成分作为输入。

尽管代码中没有显式地将PCA结果作为UMAP的输入参数传递，Seurat包的RunUMAP函数默认会使用之前通过RunPCA生成的主成分。具体来说：dims 参数：RunUMAP中的dims = 1:10指定使用PCA的前10个主成分。Seurat内部会自动从PCA结果中提取这些主成分用于UMAP降维。

5 Marker基因

差异表达基因分析，目的是找到在不同细胞群体（簇）中高表达的特征基因。

pbmc.markers <- FindAllMarkers(pbmc, 
                               only.pos = TRUE,
                               min.pct = 0.25)
                               
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
# A tibble: 18 × 7
# Groups:   cluster [9]
       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene         
       <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>        
 1 9.57e- 88       2.40 0.447 0.108 1.31e- 83 0       CCR7         
 2 1.35e- 51       2.14 0.342 0.103 1.86e- 47 0       LEF1         
 3 7.07e-139       7.28 0.299 0.004 9.70e-135 1       FOLR3        
 4 3.38e-121       6.74 0.277 0.006 4.64e-117 1       S100A12      
...

注：

FindAllMarkers()

Seurat 包中的一个函数，用于识别不同细胞群体（簇）之间的差异表达基因（markers）。

only.pos = TRUE：仅返回在给定簇中表达上调（正向表达）的基因，而不包括在其他簇中下调的基因。这通常用于识别某个簇中特异性表达的基因。

min.pct = 0.25：一个基因必须在至少25%的细胞中表达，才能被认为是潜在的标记基因。这有助于排除在大多数细胞中都不表达的基因，从而减少噪声。

group_by(cluster)

group_by(cluster)：将数据按cluster（即细胞群体）分组。这里的cluster表示每个细胞群体的标识。

top_n(n = 2, wt = avg_log2FC)：从每个簇中选出前2个基因，这些基因在该簇中的平均log2FC最高。

问：Marker基因一定是高变基因吗？

Marker基因不一定是高变基因。虽然在单细胞RNA测序数据分析中，高变基因和Marker基因经常被研究者特别关注，但它们的定义和用途是不同的。

高变基因（Highly Variable Genes, HVGs）

• 定义：高变基因是指在整个数据集中，表达水平变化较大的基因。这些基因的变异性大，通常是用于描述不同细胞类型之间的异质性。
• 用途：高变基因常用于初步的降维和聚类分析，例如PCA，因为它们能够捕捉到数据集中不同细胞群体的主要变异性。然而，高变基因的选择标准通常是全局性的，并不一定与特定细胞类型相关。

Marker基因（标志基因）

• 定义：Marker基因是指在特定细胞群体中显著高表达，且能够区分该群体与其他群体的基因。它们通常是特异性表达在某个细胞类型或状态中的基因。
• 用途：Marker基因用于细胞类型注释和生物学功能分析。通过识别Marker基因，可以帮助研究者确定特定细胞群体的生物学特性，或将其标识为已知的细胞类型。

相关可视化：比较某个基因在几个cluster之间的表达量

# 小提琴图
VlnPlot(pbmc, features = c("PPBP", "S100A9"))
# 在umap图上标记
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)
# marker基因的热图
library(ggplot2)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
DotPlot(pbmc,features = unique(top10$gene))+RotatedAxis()
RidgePlot(pbmc,features = "RPS12")

6 细胞注释

注释文件：markers.txt

a = read.delim("../supp/markers.txt",header = F)
gt = split(a[,2],a[,1]) #unstack(a[,c(2,1)])

DotPlot(pbmc, features = gt,cols = "RdYlBu") +
  RotatedAxis()

注：

split(a[,2], a[,1])：

根据第一列（簇标识）将第二列（基因名）分组。split 函数返回一个列表，每个元素包含一个簇中的所有Marker基因。

将聚类得到的细胞群体重新命名，并在UMAP图上标注这些群体的新名称。

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

names(new.cluster.ids) <- levels(pbmc)
seu.obj <- RenameIdents(pbmc, new.cluster.ids)
p1 <- DimPlot(seu.obj, 
        reduction = "umap", 
        label = TRUE, 
        pt.size = 0.5) + NoLegend()
p1

注：

new.cluster.ids

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

new.cluster.ids 是一个字符向量，其中包含了新的群体名称。这些名称是基于对每个聚类结果的生物学特征和已知Marker基因的分析得出的，反映了每个群体可能对应的细胞类型。这些名称依次对应于原始聚类的顺序。

将 new.cluster.ids 的名称与 pbmc 对象的聚类级别（即原始聚类编号）进行关联。levels(pbmc) 返回原始的聚类级别，names(new.cluster.ids) <- 将 new.cluster.ids 向量的名称设定为这些聚类编号。

即：

> new.cluster.ids
             0              1              2              3              4 
 "Naive CD4 T"   "CD14+ Mono" "Memory CD4 T"            "B"        "CD8 T" 
             5              6              7              8 
"FCGR3A+ Mono"           "NK"           "DC"     "Platelet" 

seu.obj <- RenameIdents(pbmc, new.cluster.ids)

RenameIdents 函数使用 new.cluster.ids 中的新名称替换原始的聚类编号，这会将 pbmc 对象中的每个细胞重新标记为其对应的新群体名称。seu.obj 是重新命名后的 Seurat 对象，其Idents属性现在包含新的群体名称。