单细胞 Seurat 聚类标准流程。

参考链接：

  1. Seurat - Guided Clustering Tutorial • Seurat (satijalab.org)
  2. vignettes/pbmc3k_tutorial.Rmd
  3. 单细胞分析 | Seurat基础流程 | 保姆级教程
  4. Introduction to scRNA-seq integration • Seurat (satijalab.org)
  5. Integrative analysis in Seurat v5 • Seurat (satijalab.org)

Seurat - Guided Clustering Tutorial

1 Setup the Seurat Object

从单细胞计数矩阵开始: 由 ==cellranger 10X== 流程对下机数据进行处理而产生的，里面存储的是UMI计数矩阵，矩阵中的值表示的是在每个细胞（列）中检测到的每个特征（即基因，行）的分子数。但要注意，最新版本的 cellranger 也使用 h5 file format输出，可以使用 Seurat 中的 Read10X_h5() 函数读取。

[!done] counts"矩阵"比较：

  1. UMI计数矩阵： Value是每个单细胞的每个基因的分子数(表达量，counts数) 是在单细胞水平的分析。比较的是每个细胞。
  2. RNAseq的counts矩阵： Value是每个样本的每个基因的表达量(counts数) 是在bulk样本水平的分析。比较的是每个样本。
  3. scRNAseq的"counts矩阵"：Value是每个Cluster的每个基因的表达量(counts数) 是在单细胞水平的分析。比较的是每个细胞簇。

Seurat流程是输入的单细胞的基因表达量(1)，而普通的DEG是在bulk水平的(2)，不过差异分析显然可以在单细胞分辨率下进行(3)，即Seurat聚类以后，按Cluster来选取DEG即可。

这里有关对于单细胞的差异分析的理解不够准确，这一点可使用dreamlet，见文章结尾

[!todo]

  • cellranger流程是什么？
    • 是一种用于单细胞的上游分析流程的软件
  • 单细胞下机数据的格式是什么？
    • 常用的是FASTQ
  • scRNA下机数据与DE上游分析用的bulk下机数据有何不同？
    • scRNA-seq除了FASTQ格式的原始测序数据外，包括细胞barcode和UMI（Unique Molecular Identifier）信息，用于追踪每个细胞的转录本。具体可以再看
1
2
3
4
5
6
7
8
9
10
library(dplyr)
library(Seurat)
library(patchwork)

# Set the work directory
setwd("/hanlitian/disk2/gmq/240312P/scRNA_raw_data/")

# Load the dataset
myData <- Read10X(data.dir = "2W")

读取的dataset需要包括以下3个文件：

  1. barcodes.tsv: 标记不同细胞的条形码序列，行数对应细胞数量。
  2. genes.tsv: 基因信息，分gene id（第一列）和gene symbol（第二列），列数不固定，主要作为矩阵的行名。行数对应的是基因数量
  3. matrix.mtx: 有表达量的基因的counts数量（第三列信息，即表达量信息），前三行为注释信息，3个数字依次分别为：基因位置(genes)、细胞位置(barcodes)、counts量
1
2
3
4
5
6
7
# Initialize the Seurat object with the raw (non-normalized data).
mySeurat <- CreateSeuratObject(
	counts = myData, 
	project = "2W",     # 示例pbmc数据中 project = pbmc3k
	min.cells = 3,      # pre-filter when initializing
	min.features = 200  # pre-filter when initializing
)

What does data in a count matrix look like? 是一种以 . 代替 0 存储单细胞各基因的counts数据的稀疏矩阵，可以极大的减少数据的内存占用，提升运行速度。

以上仅是读取1个样本数据并创建了1个Seurat对象，当有多个测序样本时，需要合并，同时，每个样本都需要做预处理的流程，因此先后就是一个问题：

  • 解决思路是：先对样本合并(merge)，（自动对每个样本进行的）预处理(pre-processing)，再将处理后的样本整合(integration)
  • 这样的好处是，merge以后的预处理流程避免了scale对于integration可能的潜在影响，因此之前认为的先预处理的数据量大小相同，而合并步骤的数据量更小（因为处理过了），即先预处理再merge是不对的。

2 Standard pre-processing workflow (before PCA)

以下步骤涵盖了 Seurat 中 scRNA-seq 数据的标准预处理工作流程：

  • 据质量控制指标选择和过滤细胞(QC and selecting cells)
  • 数据归一化(Normalization)
  • 缩放(Scaling)
  • 检测高变异特征(Identification of highly variable features)

2.1 QC and selecting cells for further analysis

2.1.1 Quality metrics index

Seurat中常用的质控（QC）指标包括以下3个：

  • nFeature_RNA：每个细胞中检测到的基因的数量。
    • 低质量的细胞或空液滴通常含有很少的基因；
    • 细胞双胞体或多胞体可能表现出异常高的基因计数(gene count)。

[!done] "doublets or multiplets"的nFeature和nCount

  1. 疑问：这很奇怪，nFeature应该是基因的数量(number of genes)，而不是基因表达的数量(read counts, 这其实是nCount)，那这里下面的这句"细胞双胞体或多胞体可能表现出异常高的基因计数(gene count)"和nFeature有什么关系？
  2. 误解："doublets or multiplets" 听起来像是一些nFeature值正常(不少也不多)但是nCount很大的情况，这是不是放在nCount里比较合适，在这里弄的我有点摸不清和感到混乱。
  3. 解决：通过后面的nFeature与nCount的散点关系图可以看出，在nFeature很大的时候，nCount也很大，这与下面nCount的解释是一致的，这其实提示我们，融合的"单细胞"其融合前的每个细胞的nFeature其实也是不同的，因此融合以后nFeature其实不会是正常的大小了，而是比较高的值。
  4. 对第3点所观察到的现象的补充：之前认为Feature正常而Count异常，实际上是假定了不同细胞的Feature是一样的，因此觉得融合之后，基因的种类的数量，也就是feature还是那么多，但这其实是不正确的，因为这里做的是RNAseq，每个细胞所发现的基因是这个细胞表达的基因，而非是共同的基因组。显然，不同的细胞表达的基因是不同的，所以融合的细胞不仅仅是Count值会异常偏高，其Feature值常常也是异常偏高的。
  5. 讨论：原文的 "Cell doublets or multiplets may exhibit an aberrantly high gene count" 其实说的是没有错的，但我想如果表达为 "Cell doublets or multiplets may exhibit an aberrantly high gene count and feature"更合适一些，更利于理解。
  • nCount_RNA：细胞内检测到的分子总数。
    • 即检测到的基因表达的总数量（一个基因可能会因为有多个转录本而被检测到多次），因此分子总数与基因表达数量密切相关。
    • 另外，随着测序深度的增加，单细胞所检测到的基因数量也在增加，两者是有一定关联性的。这可以在下文的FeatureScatter中直观地看到。
    • 对于RNAseq，"基因表达量"是通过RNA测定的，nCount就是read counts。
  • percent.mt：对应线粒体基因组的序列的百分比。
    • 低质量/濒死的细胞通常表现出广泛的线粒体污染。
    • 我们通常使用PercentageFeatureSet()函数来计算线粒体基因相关的reads的百分比。代码如下：
1
2
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
mySeurat[["percent.mt"]] <- PercentageFeatureSet(mySeurat, pattern = "^MT-")

关于这行代码的解释：

  • "[[]]"运算符可以向操作的Seurat对象的元数据meta.data中添加列。
  • meta.data是一个存储QC统计数据的好地方
1
2
3
4
5
head(pbmc@meta.data)  
                   orig.ident nCount_RNA nFeature_RNA
# AAACATACAACCAC-1     pbmc3k       2419          779  
# AAACATTGAGCTAC-1     pbmc3k       4903         1352  
# AAACATTGATCAGC-1     pbmc3k       3147         1129
  • 这里的第一列是细胞的barcode，代表每一个单细胞，也就是说，对于这个数据框(maybe)，行是细胞(barcode)
    • orig.ident是在创建Seurat对象时 (CreateSeuratObject())定义的project
    • 创建Seurat对象过程中会自动计算基因的数量nFeature_RNA和分子总数nCount_RNA来存储在对象的meta.data中。
  • 注意到这里选取线粒体基因的方法：pattern = "^MT-"
    • 我们可以通过查找rownames(mySeurat)中以MT开头的基因(即^MT-)的方法查看标准数据中的线粒体基因，
      • ^: 表示匹配字符串的开头。
      • MT-: 表示要匹配以"MT-"开头的字符串。
      • rownames(mySeurat)[grep("^MT-", rownames(mySeurat))]
  • 这里的MT的方法是对人的gene name使用的，小鼠的数据使用小写mt-
  • 当我们使用常规数据，\只知道线粒体基因的id，他们没有共同特征时该怎么办呢？方法就是把pattern展开来写，只有13个基因(13个仅针对官方示例的pbmc数据):
1
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "ATP6|ATP8|COX1|COX2|COX3|CYTB|ND1|ND2|ND3|ND4|ND4L|ND5|ND6")

[!todo] 获取数据的线粒体基因信息 如果不知道线粒体基因的信息怎么办呢？

  • scRNA-seq/lessons/mitoRatio.md at master · hbctraining/scRNA-seq · GitHub
  • 这个线粒体基因信息在gtf或者gff文件中都会有的，去找一下基本都会找到的，找到后再像上面一样匹配即可。
  • 我觉得这个线粒体基因如何在 gtf / gff 文件中如何寻找：
    • 有一列代表基因名(features)，可能是 Symbol or Ensembl ID 之类的
    • 有一列代表基因位置，可能是 chromosomal / mitochondrial 等等

2.1.2 Quality metrics threshold

New_trial 7里包括一些实践中阈值确认的尝试。根据其链接的教程提供了更多的可视化方案，有助于更好地确认 QC metrics threshold，但实际上仅仅根据 Seurat 标准流程进行简单的筛选就可以了，没有必要搞这么复杂，因为数据比较多的时候很难这样面面俱到，此外，严格的质控不一定是一件好事。

我们需要将质量控制指标可视化，并使用这些指标来筛选细胞。可视化的方法包括做小提琴图( VlnPlot() )和特征散点图( FeatureScatter() )。小提琴图展示了任意位置数值的密度，可以知道哪些位置的数值的密度较高，这样，我们可以可视化数据的分布，从而筛选掉低质量的离群值。

1
2
3
4
5
6
7
8
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

官方教程给出的阈值筛选：

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

根据小提琴图和散点图的情况，完全可以选择如下的质控过滤指标，即：

  • nFeature_RNA > 200（最开始的min.features = 200）& nFeature_RNA < 2000
  • nCount_RNA < 7000
  • percent.mt < 7

2.2 Normalizing the data

从数据集中移除不需要的细胞后，下一步就是对数据进行归一化处理。

[!done] 为什么做Normalize？ 排除技术误差，让测序深度和文库大小不同的细胞的基因表达量具有可比性，进而更好地呈现基因水平的差异。

默认情况下，采用全局缩放归一化方法 "LogNormalize"（对数归一化）:

  • 将每个细胞的特征表达测量值按总表达量归一化，再乘以缩放因子（默认为 10,000），然后对结果进行对数转换。
    Valuelog⁡=ln⁡(GeneCountsAllGeneCounts×ScaleFactor+1) Value_{\log} = \ln(\frac{GeneCounts}{AllGeneCounts} × ScaleFactor +1)
  • 在 Seurat v5 中，归一化值存储在 pbmc"RNA"$data。
1
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

由于这些参数都是默认的，因此也可以省略，写成：

1
pbmc <- NormalizeData(pbmc)

虽然这种归一化方法是标准的，并广泛应用于 scRNA-seq 分析，但全局缩放依赖于一种假设，即每个细胞原本含有相同数量的 RNA 分子。我们和其他人已经为单细胞预处理开发了不需要这些假设的替代工作流程:

  • SCTransform() 归一化工作流程。
  • 我们的论文对该方法进行了介绍。
  • 这里还有一个使用 Seurat 的讲解。
  • 使用 SCTransform 会取代 NormalizeData、FindVariableFeatures 及 ScaleData。

[!todo] FPKM,TPM与LogNormalize, SCTransform的区别和联系？

  1. FPKM, TPM 是bulk RNA分析的常用归一化方法。
  2. LogNormalize, SCTransform是scRNA使用Seurat的方法。

2.3 Identification of highly variable features

我们接下来计算在数据集中表现出细胞间高差异性的基因子集（即它们在一些细胞中高度表达，在另一些细胞中低度表达）。研究发现，在下游分析中这些基因有助于突出单细胞数据集中的生物信号。

FindVariableFeatures：鉴定高可变基因的函数，默认返回每个数据集的前2000个基因用于像PCA这样的下游分析。后面分析就主要用这2000个基因了，这样可以减小内存和计算资源。

1
2
3
4
5
6
7
8
9
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

[!todo] Why VST here?

  • 注意到这里有一个选择方法是vst，或许这里需要更多地来了解一下这个FindVariableFeatures是如何工作的？
  • 我记得vst与log，rlog都是数据的归一化方法，那这为啥前面用了一次LogNormalize以后又在这里写vst？
  • 在 FindVariableFeatures 中使用的 vst（方差稳定转换）方法是一种用于识别单细胞 RNA-seq 数据中高变异基因的技术。vst 的目标是通过对基因表达数据进行方差稳定化，以减小高表达基因和低表达基因之间的差异，使得高变异基因更容易被识别。

2.4 Scaling the data

在计算完高可变基因后，需要对数据集进行一个线性转换（缩放）。

这个步骤是降维（如PCA）数据处理之前的一个标准的预处理过程，由ScaleData()这个函数执行。该函数的作用解析如下：

  • 调整每个基因的表达量，使得细胞间的平均表达量为0
  • 缩放每个基因的表达量，使得细胞间的方差为1

[!todo] Why choose scale?

  • 通常来说，PCA只有在方差变化比较小的数据中表现好，而我们的RNAseq数据的特点是不符合的，通常会使用一些log或vst方法来统一方差，这里即选用了相似的scale方法，使其均值为0，方差为1。
  • 但我其实有点搞不懂这些方法的不同了，正如前面在高变基因筛选的步骤选了vst方法，这让我感觉有点混乱且对方法的选择摸不清。
  • 以上两种理解都是错误的，对于RNAseq数据分布出现的均值-方差关系，我们使用VST方法对其处理，使得方差一致，以得到更好的PCA结果，这一步骤实际上是在Normalize的时候实现的，scale不会改变数据的形状，而只改变大小。尽管如此，改变大小也是有好处的，这样可以避免高表达的基因占据过大的权重而影响PCA。
  • 结果存储在pbmc"RNA"$scale.data中
  • 默认只有上一步筛选得到的可变基因会被缩放转换，可以通过feature参数来改变操作对象的范围。
1
2
3
4
5
6
7
8
# 默认是这样的
pbmc <- ScaleData(pbmc)  
# 也可以这样运行  
scale.genes <-  VariableFeatures(pbmc)  
pbmc <- ScaleData(pbmc, features = scale.genes)
# 再比如，如果想对所有的gene都scale，运行代码如下：
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
  • 但对all.genes都进行scale是不太划算的，因为这一步的运行时间相当长。就是不全都scale一下就不那么放心。

但这一步的目的是很简单的，即赋予每个基因在下游分析中同样的权重，不至于使高表达基因占据主导地位，使得PCA的表现更好。scale的作用是使数据具有一个统一的尺度，这并不会影响每个点的相对位置，只是使他们的表达量尺度统一起来。是PCA分析的必要步骤。

在 Seurat 中，我们还可以使用 ScaleData() 函数移除单细胞数据集中不需要的变异来源。例如，我们可以 "回归 (regress out)"与细胞周期阶段或线粒体污染有关的异质性：

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

但建议使用我们的新规范化工作流 SCTransform()完成这一任务，可以参见其论文与讲解。与 ScaleData() 一样，函数 SCTransform() 也包含一个 vars.to.regress 参数。

[! done] Scale应当在merge之后，在Integration之前进行 多样本或批次的数据整合分析时，是否需要按样本分别进行ScaleData处理？

3 Perform linear dimensional reduction

对于预处理结束后的数据，就可以对其进行PCA了。

1
2
3
pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))
# 简单查看pca的结果
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

接下来进行PCA结果的可视化。可视化的作用大概就是向别人直观地进行展示我们的pca结果，他对后续的处理没有影响。 Seurat提供了几种方法来对PCA结果的细胞和基因进行可视化，包括VizDimReduce()、 DimPlot()和DimHeatmap()。

第一种：VizDimLoadings()：VizDimLoadings是一个R包中的函数，用于绘制维度加载图。维度加载图是一种用于可视化主成分分析（PCA）或其他降维方法中的变量与主成分之间的关系的图形。它显示了每个变量在每个主成分上的贡献程度，可以帮助我们理解变量在不同主成分上的权重和相关性。这种图形通常以散点图或箭头图的形式呈现。

1
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

第二种：DimPlot()：DimPlot函数可以将单细胞数据降维后的结果可视化，以便于我们观察不同细胞之间的相似性和差异性。在散点图中，==每个点代表一个单细胞==，不同颜色或形状的点表示不同的细胞类型或状态。通过观察散点图，我们可以发现细胞之间的聚类模式、相似性和差异性，进而对细胞类型或状态进行分类和注释。此外，DimPlot函数还可以根据其他降维方法（如t-SNE和UMAP）进行可视化，以便于我们比较不同降维方法的效果。

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

第三种：DimHeatmap()：DimHeatmap函数可以将单细胞数据在降维空间中的表达模式可视化，以便于我们观察不同基因或基因集之间的表达模式和差异性。在热图中，每一列代表一个细胞，每一行代表一个基因或基因集，颜色的深浅表示该基因或基因集在该细胞中的表达水平。通过观察热图，我们可以发现不同基因或基因集之间的相关性、表达模式和差异性，进而对细胞类型或状态进行分类和注释。在DimHeatmap函数中，我们可以指定降维的主成分数量和用于绘制热图的细胞数量，以便于我们控制可视化结果的精度和速度。此外，如果存在批次效应或其他技术因素导致的表达偏差，可以使用平衡化方法来消除这些效应，从而提高可视化结果的准确性和可靠性。

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

4 Determine the ‘dimensionality’ of the dataset

为了克服scRNA-seq数据的任何单一基因的广泛技术噪音，Seurat基于PCA分数对细胞进行聚类，每个PC基本上代表一个“metafeature”(==metafeature或许可以理解为一组基因，这群基因具有相似的表达pattern==)，最高主成分代表了对数据集的稳健压缩，因此，如何准确地确定降维的“维度”对后续的分析非常重要！然而，我们应该选择包含多少个成分呢？10? 20? 100?

Robust：稳健的，音译为鲁棒性，是指在面对异常情况或噪声时，系统或方法仍能保持稳定和可靠的能力。在数据分析中，鲁棒性通常指的是某个统计模型或算法对于数据中的离群点、缺失值、错误值等异常情况的容忍度。具有较高鲁棒性的模型或算法能够有效地识别和处理异常值，从而提高分析结果的准确性和可靠性。

识别数据集的真实维度--对用户来说可能具有挑战性/不确定性。因此我们建议采用以下几种方法：

  • 第一种方法更具监督性，通过探索 PC 来确定相关的异质性来源，例如可与 GSEA 结合使用。
  • 第二种（ElbowPlot）
  • 第三种是常用的启发式方法，可以即时计算。在本例中，我们可能有理由选择 PC 7-12 之间的任何数值作为临界值。

[!todo] pcs 如何与 GSEA 结合使用？？？

通常情况下，我们使用最多的就是肘形图（ElbowPlot），因为其结果展现形式很直观（在V4官网上还有一种JackStraw图，肘形图实际上也来源于这种图，只不过JackStraw图并不直观）

1
ElbowPlot(pbmc)

基于每个元素(PC)解释的方差百分比对主要元素进行排序。在这个例子中，我们可以观察到在 PC9-10周围有一个“肘”，这表明大部分真实信号是在前10个 PC 中捕获的。

或许可以这样理解：==横坐标的PC是你接下来要选择分析的维度数，即多少组基因表达pattern==。当选择一组的时候，即所有基因都在一个表达pattern里面（我们从生物学角度看也肯定是不可能的），我们很容易想到，所有细胞都在一个组内，肯定是存在很大差异的，所以方差值（方差=标准差（Standard Deviation）的平方）就会很大。当PC增多，即分组增多后，每个组内的差异也就慢慢变小了，方差值自然也就变小了。再往后，随着选的PC的增多，基本上区分不出来组内的差异了，这个时候方差基本就稳定了下来，此时在选取更大的PC值也就没有更大的意义了，并且还会增加计算压力。因此，个人认为关于PC选取的最合适原则是在拐点右侧3-5个维度即可（本数据拐点在7-8左右，我们选10刚刚好），宽松一点的话，可以设置在拐点右侧10以内。不过最重要的一点还是尽可能地多尝试，选取最适合自己的PC值，这一点也是官网教程所建议的。

总之，我们其实目的在于选择方差稳定的pc数来进行后续的聚类。（这句话理解的可能还有点问题）

我们在此选择了 10，但鼓励用户考虑以下几点：对树突状细胞和 NK 有研究兴趣的研究员可能会认识到，与 PCs 12 和 PCs 13 密切相关的基因定义了罕见的免疫亚群（例如，MZB1 是质体 DC 的标记）。然而，这些亚群非常罕见，在没有事先了解的情况下，很难将其从如此规模的数据集的背景噪声中区分出来。

  • 我们鼓励用户使用不同数量的 PCs（10、15 甚至 50 个！）重复下游分析。
  • 正如您所观察到的，结果通常不会有太大差异。
  • 我们建议用户在选择该参数时偏高一些。例如，仅使用 5 个 PC 进行下游分析确实会对结果产生重大不利影响。

[!Todo] How to understand PCA? 参考：解释PCA结果

PCA的核心思想是将原始数据投影到一个新的坐标系中，使得数据在新的坐标系下具有最大的方差。这意味着我们可以通过保留最重要的方差来保留数据的关键特征，同时将维度降低到一个更可控的水平。

具体而言，PCA的步骤如下：

  1. 数据标准化： 首先，我们对原始数据进行标准化处理，以确保每个特征具有相同的重要性。
  2. 计算协方差矩阵： 接下来，我们计算数据的协方差矩阵，该矩阵描述了数据特征之间的相关性。
  3. 特征值分解： 然后，通过对协方差矩阵进行特征值分解，我们可以得到特征向量和特征值。
  4. 选择主成分： 最后，根据特征值的大小，我们选择最重要的特征向量作为主成分，从而构建新的特征空间。

主成分的含义：

  • 主成分是原始特征的线性组合，它们是以方差最大的方向来表示数据的。在这个例子中，PC1和PC2分别代表了数据中的最大和第二大的方差。

如何选择更重要的变量：

  • 可以通过查看主成分的贡献率来判断每个主成分对原始数据的解释程度。通常情况下，我们会选择贡献率较高的主成分作为数据的代表，从而实现降维。

以下是一个基本的方法来选择更重要的变量：

  1. 查看主成分的贡献率：首先，计算每个主成分对总方差的贡献率。这可以通过主成分的方差解释比例来获得，即每个主成分的方差除以所有主成分方差之和
  2. 确定阈值：设定一个阈值，例如80%或90%的方差解释比例。这个阈值表示我们想要保留的总方差的比例。
  3. 累积贡献率：计算累积贡献率，即逐步将主成分的贡献率相加，直到达到或超过设定的阈值。
  4. 选择变量：根据累积贡献率选择变量。通常情况下，我们会选择使累积贡献率超过阈值的所有主成分对应的原始特征作为重要变量。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# 加载必要的库
library(ggplot2)library(pracma)

# 生成随机数据
set.seed(123)
x1 <- rnorm(100)
x2 <- 2*x1 + rnorm(100)
x3 <- 0.5*x1 - rnorm(100)
data <- data.frame(x1, x2, x3)

# 数据标准化
scaled_data <- scale(data)
# 应用PCA
pca_result <- prcomp(scaled_data)
# 计算主成分的方差解释比例
percent_variance <- round(100 * (pca_result$sdev^2) / sum(pca_result$sdev^2), 1)
# 计算累积贡献率
cumulative_variance <- cumsum(percent_variance)
# 设定阈值
threshold <- 80
# 选择重要变量
important_variables <- which(cumulative_variance <= threshold)

cat("根据", threshold, "% 的阈值选择的重要变量为：\n")
cat(names(data)[important_variables])

5 Cluster the cells

Seurat 在（Macosko 等人）的初始策略基础上，采用了基于图的聚类方法。重要的是，驱动聚类分析的距离度量（基于先前确定的 PCs）保持不变。不过，我们将细胞距离矩阵划分为聚类的方法有了显著改进。我们的方法在很大程度上受到了近期文稿的启发，这些文稿将基于图的聚类方法应用于 scRNA-seq 数据 (SNN-Cliq, Xu and Su, Bioinformatics, 2015) 和 CyTOF 数据 (PhenoGraph, Levine et al., Cell, 2015)。简而言之，这些方法将细胞嵌入图结构--例如 K 近邻（KNN）图，在具有相似特征表达模式的细胞之间划出边，然后尝试将该图划分为高度相互关联的 ‘quasi-cliques’ or ‘communities’。

也就是说，==Seurat 的聚类是基于图表的方法==，但聚类的==距离度量则是基于先前识别的 PCs==，因而保持不变，这些方法将细胞嵌入到图形结构中，如 K 近邻图(KNN)，在具有相似特征表达模式的细胞之间绘制边界，然后将该图划分为高度相互关联的不同簇。

  • FindNeighbors(): 与 PhenoGraph 一样，我们首先根据 PCA 空间中的欧氏距离构建一个 KNN 图，然后根据任何两个单元的局部邻域中的共享重叠（Jaccard 相似性）来细化它们之间的边权重。这一步骤使用 FindNeighbors() 函数执行，并将之前定义的数据集维度（前 10 个 PCs）作为输入。
  • FindClusters(): 为了对细胞进行聚类，我们接下来会应用模块化优化技术，如卢万算法 Louvain algorithm（默认）或 SLM (SLM, Blondel 等人，统计力学杂志)，对细胞进行迭代聚类，目的是优化标准模块化函数。FindClusters() 函数实现了这一过程，并包含一个分辨率参数，用于设置下游聚类的 "粒度"，数值越大，聚类数量越多。我们发现，将该参数设置在 0.4−1.20.4-1.2 之间，通常能为 3k3k 左右的单细胞数据集带来良好的结果。对于更大的数据集，最佳分辨率通常会提高。
  • Idents(): 可以使用 Idents() 函数获取到 clusters。
1
2
3
4
5
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)

[!todo] Relationship between dimensional reduction, find clusters, and find variable genes? 通常情况下，聚类是基于降维后的数据进行的，因为降维后的数据更容易处理和分析。但在某些情况下，也可以直接在原始高维数据上进行聚类分析，尤其是当数据规模较小或者已经进行了预处理和特征选择。聚类后在cluster之间找差异基因。

6 Run non-linear dimensional reduction (UMAP/tSNE)

尽管这一步写在了240312 Seurat_Guided Clustering Tutorial-STAR之后，但这两步实质上是没有绝对的先后，这是因为UMAP的绘制与FindCluster是独立的。

他们都是基于FindNeighbors得到的graph(s)：使用PCA数据先输入KNN算法来基于基因表达量描述细胞间的关系，得到一个 k-nearest neighbor (KNN) graph, 而后KNN图用于产生一个 undirected shared nearest neighbor (SNN) graph.

这些图输入 clustering algorithms 来使得细胞聚类在一起

也是这些图，利用 t 分布随机邻域嵌入法（t-SNE）或统一逼近和投影法（UMAP）进一步进行非线性降维，以图形方式描述这些邻域在两个维度上的结构。 CITE: The impact of package selection and versioning on single-cell RNA-seq analysis. 看 Introduction 部分 doi: 10.1101/2024.04.04.588111

把聚类放直接在降维之后，这样的顺序在整体的理解上可以更清晰一些。

Seurat 提供了几种非线性降维技术，如 tSNE 和 UMAP，用于可视化和探索这些数据集。这些算法的目标是学习数据集中的潜在结构，以便在低维空间中将相似的细胞放在一起。因此，在上面确定的基于图的聚类中组合在一起的应在这些维度缩减图上共同定位。

虽然我们和其他人经常发现，tSNE 和 UMAP 等二维可视化技术是探索数据集的重要工具，但所有可视化技术都有局限性，不能完全体现基础数据的复杂性。特别是，这些方法旨在保留数据集中的局部距离（即确保基因表达谱非常相似的细胞共定位），但往往不能保留更多的全局关系。 我们可以利用 UMAP 等技术实现可视化，但要避免仅根据可视化技术得出生物学结论。

在这一点上，韩立天师兄给出了不错的解释，可以更多地参见240302 UMAP_hlt

[!todo] 降维方法的选择？

  • 按照这句话，聚类用PCA，是否意味着PCA在全局结构的保留更好？
  • 但聚类时，选择PCA降维，使用UMAP/t-SNE等来做可视化，这是否是由于UMAP在距离上是扭曲的，而聚类的生物学假设则认为，相似的细胞在低维的空间中有更近的距离？
  • 做轨迹是用局部结构，是否意味着UMAP的表现会优于PCA，slingshot的文档中演示的这2个都有，但他们有意避免了哪种方法最优的问题（“as this likely depends on the type of data, method of collection, upstream computational choices, and many other factors.”）那如何判断、比较、取舍？
  • 有人认为：UMAP展现的距离更真实，TSNE图展示的结果似乎更美观（TSNE图呈现出来的结果中，各个细胞团都很紧促，最终呈现出来细胞都抱团分布，出来的图个人感觉更美观一丢丢）单细胞分析 | Seurat基础流程 | 保姆级教程。在这一问题上的分歧有待进一步思考，可以问下师兄当时是读了哪些文献，咱也学习一下。
1
2
3
4
5
pbmc <- RunUMAP(pbmc, dims = 1:10)

# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap")

保存数据：

1
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

[!todo] .rds or .Rdata?

  • 似乎有这样的两种保存策略，但存的区别在哪我还不知道，不过我都是直接点击保存键(x
  • .rds 更加常用，用来保存某个关键变量的结果，.Rdata几乎不用，它会保存当前的所有变量到一个环境中，慢且冗余。

7 Finding differentially expressed features (cluster biomarkers)

在 Seurat v5 中，我们使用了 PACK-presto package 以显著提高 DE 分析的速度，尤其是大型数据集的分析速度。使用方法参见：240314 Seurat_Differential expression testing

  • GitHub - immunogenomics/presto: Fast Wilcoxon and auROC 可以通过 FindMarkers() 函数访问 Seurat 的大部分差异表达特征。默认情况下，Seurat 根据非参数 Wilcoxon 秩和检验执行差异表达（DE）测试。要测试两组特定细胞间的差异表达基因，请指定 ident.1 和 ident.2 参数。

对于没有使用 presto 的用户，Seurat可以通过差异表达(DE)来获取用于定义簇的标记。默认情况下，与所有其他细胞簇相比，它标识单个簇(在 ident.1中指定)的正标记和负标记。

  • FindAllMarkers()为所有簇自动完成这个过程，但是您也可以对任意簇之间进行比较，或者是某一簇对其他所有簇进行比较。了解 min.pct 和 logfc.threshold 参数，可以提高 DE 检验的速度。

[!todo] 如何选择差异分析的统计学检验方法？

这就又一次涉及到了有关RNA-seq数据的统计学建模的问题，在bulk数据的中，我们使用的是Deseq2或edgeR或limma等，以我现在唯一用过的Deseq2，其建模假设是数据的分布服从负二项分布，然后我记得似乎是通过了log变换等使得数据服从正态分布，然后使用了Wald检验（一种参数检验）。

事实上，这三个包都是使用参数检验来做差异分析的，但这些数据使用参数检验的主要依据是RNAseq的小样本量，第3个链接即是说在大样本（n>8）时应当使用wilcoxon秩和检验，或许这与Seurat选择这个方法有关，因为scRNA数据做伪差异分析的样本量(通常是指cluster的数量)是更大的，通常都是大于8的。

上面这行的解释实在是有点搞笑，我怎么会想这么解释的。。。关于这个问题我尚未理解，留作之后来解决吧

  • 此外，用 Wilcoxon 速度是真快。
  • 只计算某个簇的marker：
1
2
3
4
5
6
7
8
# find all markers of cluster 2
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2)
# list the former 5 markers
head(cluster2.markers, n = 5)

# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3))
head(cluster5.markers, n = 5)
  • 计算所有簇的marker，并选择positive的按cluster进行展示：
1
2
3
4
5
6
# find markers for every cluster compared to all remaining cells, report only the positive
# ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)
pbmc.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
1
2
3
4
5
cluster0.markers <- FindMarkers(pbmc, 
                                ident.1 = 0, 
                                logfc.threshold = 0.25, 
                                test.use = "wilcox", 
                                only.pos = TRUE)

logfc.threshold：表示基因在两组细胞间的平均表达量的差异倍数，该参数的默认值设置在0.1，增加该值可以加速运算速度，但可能会错过较弱的信号。通常应用中我使用的是0.25。 test.use：对数据进行差异表达测试，可以理解为参数检验，一般笔者采用的方法是wilcox，筛选阈值就是大家常用的p值，也就是0.05。

  • 保存marker：
1
2
3
4
5
6
7
8
9
10
# 先把这个markers存入变量，写入变量为了后面好用  
All_markers <- pbmc.markers %>% group_by(cluster)  
  
# 保存  
write.table(all.markers[all.markers$p_val_adj <= 0.05,], 
			file = "./all_markers.xls", 
			sep = "\t", 
			col.names = TRUE, 
			row.names = F, 
			quote= F)

实际上，到了这里单细胞Seurat流程的分析基本上走完了，在这里保存RDS是比较合适的，因为里面也存好了计算的差异基因的情况，到时候重新加载RDS文件时，不需要运算即可获取。同时与注释获得完整的细胞类型还有一定的时间要求。所以，我们选择在这里保存第一版分析数据的RDS。

1
2
# 保存结果为rds文件  
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

下面的可视化也是基于以上分析结果进行的，所以这是在此保存.rds的原因。

8 Visulization

可视化可以直观地展示基因的表达情况，也是后续细胞簇注释的依据。

基因表达量的图是用来显示基因在各簇中的表达情况的，这个为该基因是否可以作为本簇的marker提供了直接的证据来源，如和已有报道相互印证（如XXX细胞特异性表达该基因），则可通过该基因对某一簇细胞进行细胞类型注释。细胞类型注释是一个人工主观的过程（自动注释的应用范围很少，只聚焦在极少数物种中，大多数的细胞类型注释还是依赖人工进行的），

以下是几种可用方法：

  • VlnPlot()
1
2
3
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
  • FeaturePlot()
1
2
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
    "CD8A"))
  • DoHeatmap()
1
2
3
4
5
6
pbmc.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
  • DotPlot()
1
2
DotPlot(pbmc, features = features) + RotatedAxis()   
# 在文章中标注markers的时候很好用，一个图能放下很多基因。
  • RidgePlot()
1
2
features <- c("LYZ", "CCL5", "IL32", "PTPRCAP")  
RidgePlot(pbmc, features = features, ncol = 2)

9 Assigning cell type identity to clusters

以下是官方教程的代码：

1
2
3
4
5
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)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
1
2
3
4
library(ggplot2)
plot <- DimPlot(pbmc, reduction = "umap", label = TRUE, label.size = 4.5) + xlab("UMAP 1") + ylab("UMAP 2") +
    theme(axis.title = element_text(size = 18), legend.text = element_text(size = 18)) + guides(colour = guide_legend(override.aes = list(size = 10)))
ggsave(filename = "../output/images/pbmc3k_umap.jpg", height = 7, width = 12, plot = plot, quality = 50)

它只是把每一个群都标好了，然后展示了如何在UMAP里面可视化

实际上，考虑分群的标注方法主要是根据已有的文献进行。并且分辨率不要太高，粗放地分一下就好。

主要的策略就是按照文献给出的Marker genes，用dotplot画一个 clusters -- marker genes 的基因表达量的图，然后根据 cluster 高表达marker的情况，给这些 clusters 注释为相应的细胞类群。

但是，总有一些 clusters 并不是很单纯的按照这些marker表达的，解决这个的方法之一是直接去看这个cluster的marker是哪些，然后看这些marker的功能是什么，从而来判断这个cluster是什么细胞亚群。应当牢记的一点是，我们关注需要的 clusters 就可以了。

对于应当考虑的细胞周期效应，和相应的计算分数，scale参见：Cell-Cycle Scoring and Regression • Seurat (satijalab.org)

Dreamlet

About Pseudobulk method

Overview

  • Quality control (seperated?)
  • Integration: harmony?
  • Get the myeloid cluster: myeloid marker?
  • Get the subclusters of myeloid cluster (monocyte, macrophage, DC, neutrophil?, NKT?, mito?)

Prepare the manipulating object

1
2
3
4
5
6
7
8
# 从 SeuratObject 转换为 SingleCellExperiment
library(Seurat)
library(SingleCellExperiment)

sce <- as.SingleCellExperiment(seurat_object)

# 从 SingleCellExperiment 转换为 SeuratObject
# seurat_object <- as.Seurat(sce)

Voom for pseudobulk

Variance Partitioning

  • For each gene it fits a linear (mixed) model and evalutes the fraction of expression variation explained by each variable. 对于每个基因，它都会拟合一个线性（混合）模型，并评估每个变量所解释的 expression variation 的比例。
  • Here we only included the stimulus status, but analyses of larger datasets can include covariates and random effects. With formula ~ StimStatus, an intercept is fit and coefficient StimStatusstim log fold change between simulated and controls. 这里我们只包括刺激状态，但对更大数据集的分析可以包括协变量和随机效应。用公式 ~ StimStatus 拟合截距和系数, StimStatusstim 模拟与对照之间的 lfc