数据降维及细胞亚群分类

生信喵

一、数据降维
        单细胞数据中包含很多细胞以及很多基因，是一个较大的数据集，维度较大，需要对数据进行降维。降维就是对原始数据进行特征提取，经常会得到高维度的特征向量。通过降维的方式来寻找数据内部的特性，提升特征表达能力，降低模型的训练成本。
        数据降维的目的：数据降维，直观地好处是维度降低了，便于计算和可视化，其更深层次的意义在于有效信息的提取综合及无用信息的摈弃；
        数据降维的方法：主要的方法是线性映射和非线性映射方法两大类。
        线性降维方法主要包括主成分分析 PCA（Principal component analysis），LDA（Linear Discriminant Analysis），非线性包括 MDS（多维缩放）和 ISOMAP（等度量映射）等。
        单细胞数据分析中常用的降维方法包括 PCA，以及 UMAP，tSNE。对多个细胞进行聚类分群。细胞亚群分类是 10X ScRNA-seq 数据分析的核心步骤，不同软件有不同的算法。因为10X ScRNA-seq 数据量大，数据结构复杂，对应的分群算法也比较复杂。
        分群的基本原理就是利用基因表达量的信息，计算各个细胞间表达模式的差异度，然后基于一定的标准将所有细胞归为多个亚群（将差异度小于值的细胞归为一个亚群）。

二、PCA 与 tSNE 以及 UMAP 之间比较
PCA：
        PCA 是线性降维的方法。以转录组数据为例，PC 的主要目的是从大量基因的表达量信息中，提取对整体基因表达量影响最大的效应（例如，实验处理效应，实验批次效应，异常样本效应等）。这些效应被称为主成分（PC），然后利用主成分绘图。PCA 分析就是将数据中大量基因的信息浓缩到少数几个变量(代表样本中的主要效应)中。而后，只要 2~3 个变量（命名为 PC1，PC2，PC3）就可以代表原来几万个基因含有的大部分信息。那么细胞之间表达量差异，就体现在 PC1、PC2 这些变量数值上的差异。
tSNE：
        t-SNE (t-distributed Stochastic Neighbor Embedding)算法是用于降维的一种机器学习算法，由 Laurens van der Maaten 和 Geoffrey Hinton 在 2008 年提出来。t-SNE 是一种用于探索高维数据的非线性降维算法，非常适用于将高维数据降维到二维或者三维，再使用散点图等基
本图表进行可视化。
UMAP：
        UMAP ，全称 Uniform Manifold Approximation and Projection，统一流形逼近与投影，是基于曼和几何代数据的理论框架结构构建的。在处理大数据集时，UMAP 优势明显，运行速度块，占用内存小。Etienne Becht 等人 2019 年在 Nature Biotechnology 上发表的一篇文章将其应用在生物学数据上并分析了 UMAP 在处理单细胞数据方面的应用和优势。UMAP 应该说是目前最好的降维算法了，现在的 10X 单细胞的降维图都选择了 UMAP，因为其能最大程度的保留原始数据的特征同时降低特征维数。
        tSNE 与 UMAP 三维展示：https://pair-code.github.io/understanding-umap/
        
        tSNE 与 UMAP 三维展示

三、为什么 tSNE 与 UMAP 更适合单细胞数据降维？
        PCA 的方法侧重于去抓住样本中隐含的主要效应，从而让差异大的样本在图中呈现较远的距离。在常规 RNA-seq 项目中，一般样本不多，实验处理效应组合数通常不会超过 10 种（例如，2 类病人× 3 个时间点取样 = 6 种处理组合），因此每个实验处理效应在所有因素的总体效应中占比都比较大，属于效应比较大的因素。
        另外，实验批次效应，离群样本等也属于比较大的效应。以上的效应都易于被 PCA 获取，因此 PCA 的方法可以良好地去展示常规 RNA-seq 项目中的处理效应、批次效应和离群样本。而单细胞测序，一次测序成千上万个细胞，可以区分出几十种细胞亚群，包括一些稀有的细胞，PCA 则无法对这些样本进行准确区分。
        tSNE 算法就属于这种可以同时兼顾局部结构和全局结构的非线性降维可视化算法。
四、PCA 分析数据降维
        PCA 分析数据准备，使用 ScaleData()进行数据归一化。进行一种线性转换，对每个基因进行转换，最终所有基因均值为 0，方差为 1。默认只是标准化高变基因（2000 个），速度更快，不影响 PCA 和分群，但影响热图的绘制。如果相对所有基因进行处理，可以通过 features选项修改。

1. pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 22000)
   
2. all.genes <- rownames(pbmc)
   
3. all.genes
   
4. 
   
5. #默认处理2000个基因，添加features对所有基因进行处理
   
6. pbmc <- ScaleData(pbmc,features = all.genes)
   
7. pbmc
   
8. #同时过滤掉不想要的基因组
   
9. # pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")

复制代码

五、PCA 分析

1. #PCA分析
   
2. pbmc <- RunPCA(pbmc,features = VariableFeatures(object = pbmc))
   
3. #输出前五个PCA结果
   
4. print(pbmc[['pca']],dims = 1:5,nfeatures = 5)

复制代码

六、绘制 PCA

1. #PCA可视化，Seurat提供三个函数VizDimReduction(), DimPlot(), DimHeatmap()
   
2. VizDimLoadings(pbmc,dims = 1:2,reduction = "pca")

复制代码

      
        前两个主成分中的基因
        绘制 PCA 散点图

1. DimPlot(pbmc,reduction = "pca" )

复制代码

      
        第一主成分散点图
        绘制两个主成分的热图

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

复制代码

      
        绘制多个主成分热图

1. #绘制15个
   
2. DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

复制代码

      


七、非线性降维
7.1 确定数据的分群个数
       在进行分群之前，需要首先定义数据集分群个数，这里我们需要选择出主成分的数目，用于后续细胞分类。这里定义的“维度”并不代表细胞类型的数目，而是对细胞分类时需要用到的一个参数。
方法 1：Jackstraw 置换检验算法；重复取样（原数据的 1%），重跑 PCA,鉴定 p-value 较小的PC；计算‘null distribution’(即零假设成立时)时的基因 scores。

1. #确定数据的分群个数，比较耗时
   
2. pbmc <- JackStraw(pbmc,num.replicate = 100)
   
3. pbmc <- ScoreJackStraw(pbmc,dims = 1:20)

复制代码

      这里我们需要选择出主成分的数目，用于后续细胞分类。这里定义的“维度”并不代表细胞类型的数目，而是对细胞分类时需要用到的一个参数。

1. JackStrawPlot(pbmc,dims = 1:15)

复制代码

      
       JackStrawPlot

方法 2：肘部图（碎石图），基于每个主成分对方差解释率的排名。

       分群个数这里选择 10，建议尝试选择多个主成分个数做下游分析，对整体影响不大；在选择此参数时，建议选择偏高的数字，些亚群很罕见，如果没有先验知识，很难将这种大小的数据集与背景噪声区分开来。

      
       碎石图检验
       JackStraw 和 Elbow 都可以决定数据的“维度”。但是 Elbow 比较直观，我们选择 Elbow 结果进行解读。可以看到，主成分（PC）12 到 15 之间，数据的标准差基本不再下降。所以我们需要在 12 到 15 之间进行选择，（官网的建议10），我们选取 15，即前 15 个主成分用于细胞的分类。


7.2 非线性降维
       基于 PCA 空间中的欧氏距离计算 nearest neighbor graph，优化任意两个细胞间的距离权重（输入上一步得到的 PC 维数）。选择不同的 resolution 值可以获得不同的 cluster 数目，值越大 cluster 数目越多，默认值是 0.5。

1. #细胞聚类
   
2. pbmc <- FindNeighbors(pbmc,dims = 1:15)

复制代码

      接着优化模型，resolution 参数决定下游聚类分析得到的分群数，对于 3K 左右的细胞，设为 0.4-1.2 能得到较好的结果(官方说明)；如果数据量增大，该参数也应该适当增大。

1. pbmc <- FindClusters(pbmc,resolution = 0.5)

复制代码

      使用 Idents（）函数可查看不同细胞的分群；
       head(Idents(pbmc), 8)

生信喵

7.3 结果可视化
       Seurat 提供了几种非线性降维的方法进行数据可视化（在低维空间把相似的细胞聚在一起），比如 UMAP 和 t-SNE。
       t-SNE (t-distributed Stochastic Neighbor Embedding)算法是用于降维的一种机器学习算法，是由 Laurens van der Maaten 和 Geoffrey Hintn 在 2008 年提出来的。t-SNE 是一种用于探索高维数据的非线性降维算法。非常适用于将高维数据降维到二维或者三维，再使用散点图等基本图形进行可视化。PCA 是一种线性算法，它不能解释征之间的复杂多项式关系；而t-SNE 是基于在邻域图上随机游走的概率分布来找到数据内的结构。SNE 通过仿射（affinitie）变换将数据点映射到概率分布上，主要包括两个步骤。
       （1）SNE 构建一个高维对象之间的概率分布，使得相似的对象有更高的概率被选择，而不相似的对象有较低的概率被选择。
       （2）SNE 在低维空间里构建这些点的概率分布，使得这两个概率分布之间尽可能地相似。t-SNE 作为新兴的降维算法，也并非万能。
其中，t-SNE 主要不足有如下几点。
       （1）t-SNE 倾向于保存局部特征，对于本征维数（intrinsic dimersionality）本身就很高的数据集，是不可能完整的映射到二到三维的空间。
       （2) t-SNE 没有唯一最优解，且没有预估部分。如果想要做预估，则可以考虑在降维之后构建一个回归方程之类的模型。但是要注意，在 t-SNE 中，距离本身是没有意义的，都是概率分布问题。

1. pbmc <- RunUMAP(pbmc,dims = 1:10,label=T)
   
2. pbmc@reductions
   
3. pbmc@reductions$umap
   
4. pbmc@reductions$umap@cell.embeddings
   
5. head(pbmc@reductions$umap@cell.embeddings)
   
6. write.csv(pbmc@reductions$umap@cell.embeddings,file = 'umap.csv') #可以导出到loupe browser中查看
   
7. 
   
8. #DimPlot绘图,可绘制umap，tsne以及pca。
   
9. #也可以使用单独的函数PCAPlot，TSNEPlot，UMAPPlot
   
10. p1 <- DimPlot(pbmc,reduction = "umap")
    
11. p1
    
12. 
    
13. #调整颜色DiscretePalette(),"alphabet", "alphabet2", "glasbey", "polychrome", and "stepped"
    
14. DimPlot(pbmc,reduction = "umap",cols = DiscretePalette(n = 8,palette = "polychrome"))
    
15. library(RColorBrewer)
    
16. brewer.pal.info# 查看不同调色板颜色数，有的调色板不够上面我们需要的8种
    
17. 
    
18. #绘制tSNE
    
19. pbmc <- RunTSNE(pbmc,dims = 1:10)
    
20. head(pbmc@reductions$tsne@cell.embeddings)
    
21. write.csv(pbmc@reductions$tsne@cell.embeddings,file = 'tsne.csv',quote = FALSE) #可以导出到loupe browser中查看
    
22. 
    
23. p2 <- DimPlot(pbmc,reduction = "tsne")
    
24. p2
    
25. DimPlot(pbmc,reduction = "tsne",cols = DiscretePalette(n = 8,palette = "polychrome"))
    
26. p1+p2
    
27. 
    
28. p3 <- DimPlot(pbmc,reduction = "pca")
    
29. p1+p2+p3
    
30. # 美化
    
31. DimPlot(pbmc,reduction = "umap") + guides(color = 'none')
    
32. library(ggsci)
    
33. DimPlot(pbmc,reduction = "umap") + guides(color = 'none') +
    
34.   scale_color_npg() #nature主题
    
35. # lancet #柳叶刀
    
36. # nejm #新英格兰
    
37. 
    
38. p1 <- DimPlot(pbmc,reduction = "umap") + guides(color = 'none')
    
39. p2 <- DimPlot(pbmc,reduction = "tsne") + guides(color = 'none')
    
40. p1+p2+p3
    
41. #保存数据集为 RDS 文件
    
42. saveRDS(pbmc,file = "pbmc_tutorial.rds")

复制代码

     
       umap 绘图
      
       t-SNE 绘图
      
       3种降维的比较
       单细胞测序比传统bulk分辨率更高，比较不同细胞亚群的差异。