来源

邓宏魁2022.4.13

nature medicine.圈图

0. 降维

• scRNA分析就是在于数据打交道。本来看不见摸不着的基因，现在体现在数据上，虽然还是摸不着，但能看得到。每个基因现在都作为数据的一个维度存在。假如我们有一个scRNA数据集，其中只有两个基因，其中两个坐标轴表示两个基因的表达量，图中的点就表示细胞，每个细胞在图上都由xy轴确定。也就是说，基因的表达决定了细胞的分布。如果有两个基因，那么细胞的位置就由两个维度决定;三个基因，就由三维空间决定;而我们有2万多个基因，那么就存在2万多个维度（高维空间)。
• 单细胞分析的一个重点就是: 用尽可能少的维度来展示数据真实的结构。降维的过程其实就是去繁（基因的变化小）存简，每个基因的变化都对整体有影响，对细胞来讲都是一个变化维度。但这么多维度我们看不了也处理不了，需要尽可能保证真实差异的前提下减少维度的数量，因此需要挑出那些更能代表整体差异的基因。
• 如何挑选有价值的基因(feature)
  一般来说，如果一个基因在细胞群体中变化幅度很大，就是受关注对象，我们也会认为是生物因素导致了这么大的差异。这样的基因叫做高度变化基因(highly variable genes ,HyGs)，y降维和聚类分群也是常用的操作，它们都是依赖于基因表达量（例如综合细胞中各个基因的表达量，然后把表达模式相近的细胞聚在一起)。因此挑选哪些基因进行聚类分群和降维分析，这是非常关键的，挑选的基因一定要有代表性，尽可能多的包含生物信息而剔除其他技术噪音。
• 官方推荐使用2000个高变基因做后续的降维

1. 挑选高变基因

FindVariableFeatures(function)

scRNA1 <- FindVariableFeatures(scRNA1, selection.method = "vst", nfeatures = 3000) 
#nfeature的小提琴图中，大部分细胞的基因表达在2000左右，高变基因阈值越高，丢失的信息会越少（维度越多，信息量越多）
#但是选的太多，降维幅度不大，挑选高变基因意义也不会太大
top10 <- head(VariableFeatures(scRNA1), 10) 
#把top10高变基因挑选出来
plot1 <- VariableFeaturePlot(scRNA1) 
#画出不带标签的高变基因图
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
#把top10的基因加到图中
plot <- CombinePlots(plots = list(plot1, plot2),legend="bottom") 

plot 

S4对象里的高变基因：var.features即高变基因

• 结果

  
  
  横坐标：基因的平均表达量；纵坐标：每个基因的变异程度

2.scaledata

#如果内存足够最好对所有基因进行中心化
scale.genes <-  rownames(scRNA1)
scRNA1 <- ScaleData(scRNA1, features = scale.genes)
##如果内存不够,可以只对高变基因进行标准化（浅推荐），但是后续热图只能查看高变基因的变化
scale.genes <-  VariableFeatures(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)

• 什么是中心化处理？
  Scaling the data(数据比例伸缩)
  

  
  ScaleData函数:
  改变每个基因的表达，使细胞间的平均表达为0
  测量每个基因的表达，使细胞间的差异为1
  我们一般将平均值为0，方差值为1的数据认为是标准数据。改变每个基因的表达，使得跨细胞的平均表达为0
  缩放每个基因的表达，以便跨细胞的方差为1
  该步骤在下游分析中给予相同的权重，因此高表达的基因不占优势

• scRNA对象中原始表达矩阵经过标准化和中心化之后，已经产生了三套基因表达数据，可以通过以下命令获得

GetAssayData(scRNA,slot="counts",assay="RNA")
#原始表达矩阵
GetAssayData(scRNA,slot= "data",assay="RNA")
#标准化之后的表达矩阵
GetAssayData(scRNA,slot="scale.data",assay="RNA")
#中心化之后的表达矩阵

Normalize函数：会使用count数据（原始数据），将归一化后的结果放入data数据

• 什么是稀疏矩阵（dgcMatrix）
  0是用点进行表现的（节省内存），查看此类数据结构时应该用as.dataframe函数

3.检查细胞周期对降维聚类的影响

细胞周期回归: 上一步找到的高变基因，常常会包含一些细胞周期相关基因。它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。

• ?CaseMatch——Match the case of character vectors

#cc.genes包内自带数据；s期的基因有哪些（s.genes），g2m期的基因有哪些（g2m.genes）
cc.genes
#将cc.genes的数据组成一个新的character，与高变基因match（PCA降维只用了高变基因）
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA1))

#g2m基因与高变基因匹配
g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match =rownames(scRNA1))
#s基因与高变基因匹配
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = rownames(scRNA1))
#cellcyclescoring对每个细胞的细胞周期进行打分
scRNA1 <- CellCycleScoring(object=scRNA1,  g2m.features=g2m_genes,  s.features=s_genes)
#得到周期基因的PCA降维结果
scRNAa <- RunPCA(scRNA1, features = c(s_genes, g2m_genes))
p <- DimPlot(scRNAa, reduction = "pca", group.by = "Phase")
p

#scRNAb <- ScaleData(scRNA1, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA1))
#vars.to.regress 可消除周期对降维的全部影响

结果：

三个不同周期的细胞如果聚集在一起，说明细胞周期对降维的影响较小

4.PCA降维

PCA降维的信息放到reduction中；8684*50，50表示50个dimension

• PCA
  PCA (Principal Component Analysis，主成分分析)是一种降维方法，致力于解决三类问题:
  降维可以缓解维度灾难问题;
  降维可以在压缩数据的同时让信息损失最小化;
  理解几百个维度的数据结构很困难，两三个维度的数据通过可视化更容易理解。
• 随着维度的增加，数据的稀疏性会越来越高。在高维向量空间中探索同样的数据集比在同样稀疏的数据集中探索更加困难。主成分分析也称为卡尔胡宁-勒夫变换(Karhunen-Loeve Transform)，是一种用于探索高维数据结构的技术。PCA通常用于高维数据集的探索与可视化，还可以用于数据压缩，数据预处理等。PCA可以把可能具有相关性的高维变量合成线性无关的低维变量，称为主成分。新的低维数据集会尽可能的保留原始数据的变量。
• PCA把相关的基因合并到“metagene”或主成分(PC)中。PC1解释最大的数据差异，具有最大的标准差（例如对于一个实验，细胞之间30%的差异由定义了PC1的基因解释)，PC2解释了数据的第二大部分差异(例如，细胞之间20%的差异可归因于PC2中的基因，而8%则归因于PC3中的基因)，然后依此类推，简单来说PC的排名就是解释数据差异贡献的顺序，其中PC1是排名最高的PC，同时也说明PC排名越低，对解释数据差异的贡献就越小。

5.聚类

nature原文：Cluster-level quality control was performed after the standard Seurat clustering pipeline using the following functions in order ： FindVariableFeatures with 2,000 genes , ScaleData , RunPCA , FindNleighbors with the first 20 PCs and FindClusters with resolution 1/otherwise defaultsettings.

• 第一步：选择合适的PCs数
  为了克服scRNA序列数据单一特征中的广泛技术噪音，Seurat根据其PCA分数对细胞进行聚类，每个PC基本上表示一个“元特征”，该特征结合了相关特征集上的信息。因此，最主要的主成分代表数据集的强大压缩。但是, 我们应该选择多少个PC? 10个?20?还是100?
• 利用elbow point选择

ElbowPlot(scRNA1, ndims=20, reduction="pca") 

elbow point就是在它之前变化幅度很大，之后变化幅度很小，属于一个转折点。启发式评估方法生成一个Elbow plot图。在图中展示了每个主成分对数据方差的解释情况(百分比表示)，并进行排序。根据自己需要选择主成分。

纵坐标：维度对于整体数据结构特征的贡献度；横坐标：维度数目

• 每个PCs都能捕获一些生物差异，而且前面的PC比后面的PC包含的差异信息更多，更有价值。就像"二八准则"，仅有20%的变因操纵着80%的局面，PCA中也一样。于是当我们从头到尾观察每个PC的差异贡献时，会有一个明显的落差，落差之前的那些PC就是"二八准则"中20%的变因。
• 如果选的PC多，可以避免丢弃一些生物信号，但同时又增加了噪音的风险。很多有经验的人会将PC数设置在一个合理的区间，例如10-50之间
• 曲线下面积代表总的贡献度
• 注意: 鉴别数据的真实维度不是件容易的事情，除了上面两种方法，Seurat官方文档还建议将主成分(数据异质性的相关来源有关)与GSEA分析相结合。如Dendritic cell和NKaficionados可能识别的基因与主成分12和13相关，定义了罕见的免疫亚群(i.e.MZB1 is a marker for plasmacytoid DCs)。如果不是事先知道的情况下，很难发现这些问题。
  Seurat官方文档因此鼓励用户使用不同数量的PC(10、15甚至50)重复下游分析。其实观察到的结果通常没有显著差异。因此，在选择此参数时，可以尽量选大一点的维度，维度太小的话会对结果产生不好的影响。

第二步、聚类分析

• 常见的聚类如层次聚类、k-means等，以及为单细胞专门开发的SC3聚类。但是Seurat采用的是谱聚类(Spectral Clustering ) 。Seurat中的谱聚类基于共享最近邻图(SNN）和模块化优化的聚类算法来识别细胞簇。它首先计算k最近邻(k-nearestneighbors)并构造SNN图，然后优化模块化功能以确定具体类群。这个算法最早在2013年发表在物理学的期刊The European Physical Journal B上。

• 第一个难点是谱聚类是一种基于图论方法
  图论方法中的"图“和图像的"图""不一样，它属于离散数学。谱聚类是从图论中演化出来的算法，它的主要思想是把所有的细胞看为高维空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低; 距离较近的两个点之间的边权重值较高。然后通过对所有数据点组成的图进行"切图"，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。

• 第二个难点是如何把我们的单细胞的数据构建成为无向有权图
  谱聚类中的图采用的是无向有权图，也就是说，每个细胞（图中的点）之间的连接是没有方向性的。但是每条边上是有权重信息的。类比一下地图，两地之间的路是没有方向(双向的)，但是地点与地点之间路的远近距离是不同的。我们手中最初只有细胞的表达量矩阵，这个每个细胞中各个基因的表达情况，属于“节点信息"。

那么哪些细胞之间需要有边相连，边的权重又该如何计算呢?

• 这里引入了邻接矩阵(W)的概念，我们首先定义一种计算细胞之间距离度量的方法，来计算每条边之间的权重wij，又称相似性sij，然后获得整个的邻接矩阵W。

• 终于到了“切图"的环节，实际就是将图切成相互没有连接的k个子图。目标就是让子图内的点权重和高，子图间的点权重和低。实际就是一个最优化问题。可以理解为通过D，L这两个矩阵所给的信息，在图的部分边进行切分，让完整的图变成许多子图，每一个子图就构成了一个类群。
• resolution参数决定类群的“粒度”，官网建议将此参数设置在0.4-1.2之间，通常可为包含大约3000个cell的数据集会返回良好的聚类结果。对于较大的数据集，最佳resolution值通常会增加。
  （resolution越高，切的刀数越多）

nature原文：FindVariableFeatures with 2,000 genes, ScaleData, RunPCA, FindNeighbors with the first 20 PCs and FindClusters with resolution 1, otherwise default settings

pc.num=1:20#确定维度
#根据20个维度的信息构建无向有权图
scRNA1 <- FindNeighbors(scRNA1, dims = pc.num) 
#resolution=1.0进行切图，目的是寻找cluster，cluster的结果会存放到graph中
scRNA1 <- FindClusters(scRNA1, resolution = 1.0)

一共得到18个level，即18个cluster

系统发育分析(Phylogenetic Analysis of ldentity Classes)

【但是似乎很少有paper放这个】
nature原文：
Constructs a phylogenetic tree relating the ’average‘ cell from each identity class.
Tree is estimated based on a distance matrix constructed in either gene expression space or PCA spac

scRNA1<-BuildClusterTree(scRNA1)
PlotClusterTree(scRNA1)

系统发育分析

6.可视化

• 降维可以定义为将高维数据降低为低维同时尽可能保证数据的全局性，降维技术是机器学习中大数据预处理的一种重要技术。
• 降维的算法主要分为两类: (1)在数据中保持距离结构的算法——线性降维(PCA、MDS和Sammon映射)，(2)在全局距离上保持局部距离的算法（流形学习)——非线性降维（t-SNE、Jsomap、LargeVis、Laplacian特征映射和difusion映射)。
• 非线性降维——这个目的是为了可视化，而不是特征提取（PCA)，虽然它也可以用来做特征提取。tSNE和UMAP目前主要用于可视化。用于可视化的降维必然涉及信息丢失并改变细胞之间的距离。因此tSNE/UMAP图应仅只用于解释或传达基于更精确的、更多维度的定量分析结果。这样可以保证分析充分利用了压缩到二维空间时丢失的信息。假如二维图上呈现的细胞分布与使用更多数目的PC进行聚类获得的结果之间存在差异，应倾向于相信前者（聚类）的结果。

tSNE和UMAP的区别
t-分布式随机邻居嵌入(t-SNE: t-Distributed stochastic neighbor embedding)是一种常见的可视化方法。它使用机器学习的算法来降低维数，非常适合将高维数据放到二维或三维空间中可视化展示，并且不会丢失细胞之间的相对距离的信息。
例如，如果发现用七个PC可以很好地表示细胞的多样性，就得需要七个轴或维度来展示细胞的空间分布。t-SNE能维持细胞在七维空间的关系并在二维图上展示细胞，所以在七维图上相邻的细胞在二维图上仍然相邻。同时PCA分析是线性的。t-SNE是非线性降维方法。
统一流形逼近与投影(UMAP：Uniform Manifold Approximation and Projection)是一种新的降维技术，其理论基础是黎曼几何和代数拓扑。相对于T-SNE降维，UMAP的优点在于:(1)能够尽可能多的保留全局结构（2）耗时更短 (3）对嵌入维数没有限制，可以扩展到更大的维度的数据集。而T-SNE对数据的维度有所限制（Hinton的T-SNE实验中先用PCA降到50,再进一步的使用T-SNE，密集的数据才能达到较好的可视化效果)。T-SNE是建立在二维的基础上，而UMAP则是建立在三维的黎曼流行上，其相关定义、计算方法都会发生一定的变化。

• t-SNE图

#先选择多少个维度
scRNA1 = RunTSNE(scRNA1, dims = pc.num)
embed_tsne <- Embeddings(scRNA1, 'tsne')
write.csv(embed_tsne,'embed_tsne.csv')
plot1 = DimPlot(scRNA1, reduction = "tsne") 

plot1
#label = TRUE把注释展现在图中
DimPlot(scRNA1, reduction = "tsne",label = TRUE) 
#把cluster都标记在图中
ggsave("tSNE.pdf", plot = plot1, width = 8, height = 7)
#保存

结果：

t-SNE图

• UMAP图

scRNA1 <- RunUMAP(scRNA1, dims = pc.num)
embed_umap <- Embeddings(scRNA1, 'umap')
write.csv(embed_umap,'embed_umap.csv') 
plot2 = DimPlot(scRNA1, reduction = "umap") 

结果：

UMAP图

pipline之后：0/1/2...17cluster分别对应的是哪种细胞？——细胞注释