引言
本系列开启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降维
在确定了变异性大的基因并缩放数据之后,原则上可以开始分析细胞的异质性。但是,强烈建议在进行更深入分析之前先进行线性降维,有时这甚至被认为是必要的步骤。这样做的好处包括但不限于:
- 数据变得更加简洁,计算速度因此得到显著提升。
- 由于单细胞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:中脑与后脑交界处的标志物
通过这些信息,可以大致了解这个数据集中包含了哪些类型的细胞。