R中单细胞RNA-seq数据分析教程 (3)

数据科学工厂

引言

本系列开启R中单细胞RNA-seq数据分析教程，持续更新，欢迎关注，转发！

5. 数据缩放

由于各个基因的基础表达水平和分布情况各异，如果不对数据进行转换，那么在分析中每个基因所起的作用就会有所不同。不希望出现这种情况，因为不希望分析结果只受到那些表达量高的基因的影响。因此，会对数据进行标准化处理，利用所选择的特征来进行，这在数据科学领域是一种常见的做法。

seurat <- ScaleData(seurat)

在这个阶段，可以通过设定 var.to.regress 参数，来排除数据集中不必要的变异因素。比如，

seurat <- ScaleData(seurat, vars.to.regress = c("nFeature_RNA", "percent.mt"))

在数据分析中，通常会考虑排除一些变量，以减少数据中的不必要变异。这些变量包括检测到的基因/转录本的数量（nFeature_RNA / nCount_RNA）、线粒体转录本的百分比（percent.mt）以及与细胞周期相关的变量（具体见下文）。做法是，首先建立一个线性回归模型，将基因的标准化表达水平设为因变量，而将那些需要排除的变量设为自变量。然后，取这个线性模型的残差，这些残差就代表了在去除这些变量线性影响后的信号。需要注意的是，这种排除变量的过程会显著减慢整个分析的速度，并且不能保证结果一定会令人满意，因为不需要的变异可能并非线性关系。因此，一个普遍的建议是，在初次探索数据时不要进行任何变量排除，而是先检查结果。如果发现有某种不需要的变异主导了细胞的异质性，再尝试排除相应的变量，并观察结果是否有所改善。

6. PCA降维

在确定了变异性大的基因并缩放数据之后，原则上可以开始分析细胞的异质性。但是，强烈建议在进行更深入分析之前先进行线性降维，有时这甚至被认为是必要的步骤。这样做的好处包括但不限于：

1. 数据变得更加简洁，计算速度因此得到显著提升。
2. 由于单细胞RNA测序（scRNA-seq）数据天生稀疏，对相关特征的测量进行汇总可以大幅提高信号的稳定性。

有什么不足之处吗？基本上没有。当然，你需要在脚本中添加一些额外的代码，并决定在后续分析中使用多少降维后的维度，仅此而已。

对于单细胞RNA测序数据，所谓的线性降维通常是指主成分分析，简称PCA。

seurat <- RunPCA(seurat, npcs = 50)

对于一个数据集，可以计算的主成分（PCs）的数量，取决于高变异基因的数量或细胞的数量，取两者中的较小值。但是，这些主成分中大多数并不提供有用的信息，它们仅仅代表了随机噪声。只有最重要的那些主成分才是有信息量的，它们揭示了不同细胞群体间的差异。因此，Seurat采用截断的主成分分析（PCA），默认情况下只计算前50个最重要的主成分，而不是计算所有可能的主成分。这个数量可以通过设置 npcs参数来调整。

即便如此，也不一定需要使用所有计算出的主成分。决定使用多少个最重要的主成分是一种技巧，没有统一的标准，每个人都有自己的见解。所谓的肘部图方法可以帮助做出选择。这种方法涉及绘制每个主成分解释的变异量，并在曲线的拐点处确定使用多少个主成分。

ElbowPlot(seurat, ndims = ncol(Embeddings(seurat, "pca")))

根据定义，高等级的主成分（PCs）比低等级的主成分在数据中解释了更多的变异（具有更高的标准差）。但是，标准差的减少并不是线性的。肘部图的曲线在最初的几个PCs时急剧下降，随后减缓并趋于平缓。可能会认为曲线的初期阶段代表了与细胞群体间生物学差异相关的‘真实’信号，而后期阶段主要反映了技术变异或单个细胞的随机性。从这个观点出发，选择前15个PCs可能是合适的，而排名低于20的PCs看起来不太必要。然而，尽管这是一个不错的参考，但并不完美：

• 精确确定曲线的肘点或转折点非常困难，因为通常它并不是一个完美的肘部形状。
• 高等级的PCs确实比低等级的PCs解释了更多的变异，但更多的解释变异并不一定意味着更高的信息量。有时候，有趣的但微弱的信号可能隐藏在噪声中，因此成为低等级PCs的一部分。

在Seurat中还有一个名为 JackStraw的程序，它也可以帮助确定在后续分析中应该考虑多少个PCs。然而，根据经验，这个程序运行非常缓慢，因为它依赖于数据的置换，而且它基本上没有提供比肘部图更多的信息。它的作用是估计每个PC的统计显著性，但一个‘显著’的PC并不一定意味着它是有信息的。随着细胞数量的增加，越来越多的PCs在统计上变得‘显著’，即使它们解释的变异并不显著。

除了做出无偏见的决定外，还可以检查哪些基因主要对每个高等级的PCs有贡献。如果知道这些基因及其在分析样本中的生物学意义，这将是有帮助的。这为提供了理解每个高等级PCs生物学意义的机会，使能够选择那些包含有用信息的PCs。

PCHeatmap(seurat, dims = 1:20, cells = 500, balanced = TRUE, ncol = 4)

请留意，仅选择那些由“有趣”基因代表的主成分（PCs）并不被推荐。这样做很可能会忽略掉一些有趣但出人意料的发现。

以这个例子来说，将采用前20个主成分进行接下来的分析。同样，使用更多或更少的主成分是完全可以的，实际上，这通常需要经过几次尝试才能做出最终的选择。同时，对于大多数数据而言，主成分的数量在10到50之间是合理的，而且在很多情况下，这个数量不会对结论产生太大的影响（但有时也会产生影响，因此仍需小心谨慎）。

7. 非线性降维

线性降维技术有利有弊。其优势在于，每个主成分（PC）都是基因表达的线性组合，这使得对PCs的解释变得简单直接。此外，数据在压缩过程中不会失真，因此大部分信息得以保留。然而，其缺点在于，通常需要超过10个主成分来捕捉大部分信息，这对于大多数分析是可行的，但在可视化时就不太适用了，因为对于普通人来说，超过三维的空间是很难理解的。

为了克服这一难题，引入了非线性降维技术。在单细胞RNA测序数据分析中，最常使用的非线性降维方法包括t分布随机邻域嵌入（t-SNE）和均匀流形近似与投影（UMAP）。这两种方法都旨在将样本放置在低维空间（2D/3D）中，以便在低维空间中尽可能保留原始空间中样本（此处指细胞）之间的距离或邻域关系。

seurat <- RunTSNE(seurat, dims = 1:20)
seurat <- RunUMAP(seurat, dims = 1:20)

plot1 <- TSNEPlot(seurat)
plot2 <- UMAPPlot(seurat)
plot1 + plot2

tSNE和UMAP哪个更好一直是热议的话题。根据经验，这两种方法各有千秋，没有哪一种总是比另一种更胜一筹。当细胞分化成不同的群体时，tSNE能提供极佳的可视化效果；而当数据中存在连续变化，比如细胞在发育和分化过程中的连续状态变化时，UMAP能更好地保持这种连续性结构。因此，尝试这两种方法，并选择更适合你数据的那一种，是一个明智的选择。

创建了tSNE或UMAP的嵌入之后，就可以开始检查数据中是否包含特定的细胞类型或状态，方法是绘制一些已知的、与感兴趣细胞类型相关的典型标记物的特征图。

plot1 <- FeaturePlot(seurat, c("MKI67","NES","DCX","FOXG1","DLX2","EMX1","OTX2","LHX9","TFAP2A"),
                     ncol=3, reduction = "tsne")
plot2 <- FeaturePlot(seurat, c("MKI67","NES","DCX","FOXG1","DLX2","EMX1","OTX2","LHX9","TFAP2A"),
                     ncol=3, reduction = "umap")
plot1 / plot2

对于不太了解这些基因的朋友，以下是它们的含义：

• MKI67：标志着细胞周期中的G2M阶段
• NES：神经前体细胞的标志物
• DCX：广泛用于识别未成熟神经元的标志物
• FOXG1：大脑皮层发育的标志物
• DLX2：腹侧大脑皮层的标志物
• EMX1：背侧大脑皮层（即大脑皮层）的标志物
• OTX2：中脑和间脑神经元的标志物
• LHX9：间脑和中脑神经元的标志物
• TFAP2A：中脑与后脑交界处的标志物

通过这些信息，可以大致了解这个数据集中包含了哪些类型的细胞。

微信扫一扫分享文章