在scRNA-seq数据通过一些列的预处理,质控以及标准化之后,后续分析步骤需要通过识别细胞间的基因表达差异来进行聚类, 分析不同细胞群的差异性。这就涉及到了单细胞RNA数据处理中的特征选择,降维以及如何使数据可视化。
数据降维的必要性
如果将单个细胞看作一个数据点,那么检测的基因数就是其对应的变量数,也就是我们所说的维数。一般一个人类scrna-seq 能检测到~25k基因,于是基因表达数据也被称为高维数据。我们知道, 并不是所有的基因都会表达,并且由于单细胞测序技术的限制,不是所有的转录分子都能被成功捕获,再加上测序深度的差异, 每个细胞中约能检测到10%~50%的转录分子,这导致了许多基因计数为0. 虽然通过数据的预处理我们能过滤掉零计数基因,剩下的基因维数仍然可以高达~15k。并且基因之间的高度相关性也使得表达数据包含了许多冗余信息, 从而掩盖了真正的生物学差异。因此,我们需要通过降维抽取数据的概要(component),减少数据噪音,并提高下游分析速度。降维的另外一个好处就是可以在低维中实现数据的可视化。
特征选择 (feature selection)
在scRNA-seq中常见的降维第一步是特征选择 (feature selection)。在很多pipeline中,会用feature来指代基因。在不同细胞之间,由于大部分的细胞表达是类似的,只有少数基因表达产生差异,特征选择的目的是筛选出高度可变基因(highly variable genes),用来代表数据的主要差异性。一般会选择1k~5k个基因,具体数目根据数据的复杂性而有所不同。在实际操作中,如果已知样本包含多种细胞亚型,如免疫细胞,我们建议大家尽量提高特征数目。
一种简单的特征选择方法是对每个基因基于其在所有细胞中的平均表达值来分组,每组中具有最高variance-to-mean ratio的基因被选为高度可变基因。单细胞分析中常用的r
R包Seurat就是使用这种方法,Seurat:: FindVariableFeatures()。
降维 (Dimensionality reduction)
在特征选择之后,可以通过降维算法对高度可变基因的表达矩阵维数进一步压缩,常见的降维方法有PCA (principle component analysis), MDS (multi-dimensional scaling), Sammon mapping, Isomap ,t-SNE (t-distributed stochastic neighbor embedding)和UMAP(Uniform Approximation and Projection method)。其中PCA, t-SNE和UMAP在scRNA-seq中使用非常普遍 [2]。
大家对于PCA应该都不陌生,它通过基因线性组合来捕获数据方差, 在实现方差最大化的同时实现降维。根据定义,前几个top 成分捕获了数据的主要差异,我们可以通过将下游分析限制在top PC上来进行降维, 这些top PC在低维构建了一个原始数据的最佳近似值。PCA应用于scRNA-seq数据时,我们的前提假设是生物过程是通过影响多个基因,以协调的方式来实现的。这也意味着,通过考虑多个基因的相关行为,可以捕获更多的差异,从而使得top PC尽可能的代表生物学结构。当然,PCA采用的线性组合方式也带来一个缺点,就是主成分自身的生物学意义很难解释。 理想情况下,我们希望降维后每个每个维度都能对应一个生物学结构, 但大多数时候是很难对应上。尽管如此,PCA方法的简单高效使其在scRNA-seq中得到了广泛使用。通常我们会使用PCA来进行一般性总结, 特别是用来查看有没有一些outlier cells,可能是在预处理时漏掉的low quality cells。
常用的R包以及方法有:
1. Scater::runpca()
2. Seurat::RunPCA()
那么我们应该如何选择主成分的个数呢?选择的PC越多,越能避免漏掉生物学信号,而代价就是可能引入更多的噪音。很多时候,大家会选择一个“合理”的任意值,通常在10~50 之间。 我们可以通过多次尝试来选到一个满意的数值。 除此之外,还有一些data-driven的方法来帮助我们选择, 例如下图中的elbow plot, 通过比较列出不同pc对应的方差百分比来选择曲线中的拐点elbow point作为“最佳”pc值。这个策略基于的假设是,捕获生物学信号的每个top PC都应比剩余的pc解释更多的variance,因此当我们越过“最佳” PC点时,其包含的方差百分比会急剧下降,形成一个明显的拐点, 如下图中的PC=10。在实际的操作中,拐点有可能不是非常的明确,特别是如果对下游的聚类效果不太满意时,可以回来调节一下pc值,看看结果会不会有显著的变化。
常用的r 包以及方法有:
1. PCAtools::findElbowPoint()
2. Seurat:: ElbowPlot()
t-SNE和UMAP是另外两种非线性的降维方法,由于其漂亮的可视化效果,这两种方法在单细胞数据教程中非常受欢迎。与PCA不同,这两种方法不仅限于线性变换,也不受制于准确表示远距离种群之间的距离,因此使得它们在低位空间如何对细胞进行排列具有更大的自由度,从而在其二维的可视化效果图中能够将复杂的细胞群分成许多不同的簇,使得效果图比PCA更加直观和容易解释。
t-SNE是以牺牲整体结构(global structure)为代价,着重于捕获局部结构(local similarly),因此它可能会夸大细胞群体之间的差异而忽略群体之间的潜在联系。如果简单的通过t-SNE结果图来解读细胞簇之间的关系,产生的结论可能会被图中的集群大小和位置所误导。在效果图中,t-SNE倾向于将密集的簇膨胀,并且压缩稀疏的簇,因此我们不能简单的通过图上的大小来衡量细胞群的差异。并且由于t-sne无法保证能保留距离较远的簇的相对位置,我们也不能简单通过图中的位置来确定远距离细胞簇之间的关系。相比之下,UMAP与t-SNE类似,同时在低维空间保留了高维空间细胞间的关系,因此UMAP更好的保留并反映了细胞群潜在的拓扑结构,对于细胞轨迹推断(trajectory inference)分析来说更实用 [4]。
在下图中,我们可以看到UMAP的可视化更趋向于一种紧凑的视觉效果,群簇之间的空间更大,也保留了更多的global structure,因此在选择可视化图中大家可以根据具体的需要来选择。根据小编的实战经验,如果对轨迹推断没有要求,只是看细胞簇群,UMAP的效果会更干净,但是在细胞数目非常大情况下,由于UMAP最大限度的保留了全局结构,这也使得每个簇群的分辨率降低,可能会使得簇群重叠,从而遮盖一些小的簇群,而t-SNE通常能将所有簇群尽可能的“铺开”,所以这种情况建议大家两种都画,然后比较一下。从速度方面来说,同一个数据UMAP的速度要比t-SNE快,这也是UMAP变得更受欢迎的重要原因。
常用的r 包以及方法有:
1. Scater:: runTSNE(); Scater:: runUMAP()
2. Seurat:: RunTSNE(); Seurat:: RunUMAP()
小编总结
对于数据降维处理可以主要归纳为两方面, 总结(summarization )和可视化( visualization)[3]。pca可以作为一个总结目的,用来查看有没有一些离群的细胞,它们有可能是一些低质量细胞,需要做进一步筛除。而t-sne和umap作为可视化工具,可以用来探索细胞群簇的关系,V信搜索:作图丫,可获取更多精彩内容。
参考文献
1. Seurat pipeline, https://satijalab.org/seurat/vignettes.html .
2. Kulkarni A, Anderson AG, Merullo DP, Konopka G. Beyond bulk: a review of single cell transcriptomics methodologies and applications. Curr Opin Biotechnol. 2019;58:129-136. doi:10.1016/j.copbio.2019.03.001.
3. Luecken MD, Theis FJ. Current best practices in single-cell RNA-seq analysis: a tutorial. Mol Syst Biol. 2019;15(6):e8746. Published 2019 Jun 19. doi:10.15252/msb.20188746.
4. Amezquita, R.A., Lun, A.T.L., Becht, E. et al. Orchestrating single-cell analysis with Bioconductor. Nat Methods 17, 137–145 (2020). https://doi.org/10.1038/s41592-019-0654-x.