Seurat

Seurat整合流程与原理

1. 使用CCA分析将两个数据集降维到同一个低维空间，因为CCA降维之后的空间距离不是相似性而是相关性，所以相同类型与状态的细胞可以克服技术偏倚重叠在一起。

2. CCA降维之后细胞在低维空间有了可以度量的“距离”，MNN(mutual nearest neighbor)算法以此找到两个数据集之间互相“距离”最近的细胞，Seurat将这些相互最近邻细胞称为“锚点细胞”。我们用两个数据集A和B来说明锚点，假设：

   • A样本中的细胞A3与B样本中距离最近的细胞有3个（B1,B2,B3）
   • B样本中的细胞B1与A样本中距离最近的细胞有4个（A1,A2,A3,A4）
   • B样本中的细胞B2与A样本中距离最近的细胞有2个（A5,A6）
   • B样本中的细胞B3与A样本中距离最近的细胞有3个（A1,A2,A7）

MNN示意图

从图中可以看出，A3与B1是“相互“最近邻的细胞，而A3与B2、B3不是相互最近邻细胞，A3+B1就是A、B两个数据集中的锚点之一。实际数据中，两个数据集之间的锚点可能有几百上千个。

3. 理想情况下，相同类型和状态的细胞才能构成配对锚点细胞，但是异常的情况也会出现，如C图中query数据集中黑色的细胞团。它在reference数据集没有相同类型的细胞，但是它也找到了锚点配对细胞（红色连线）。Seurat会通过两步过滤这些不正确的锚点：

   • 在CCA低维空间找到的锚点，返回到基因表达数据构建的高维空间中验证，如果它们的转录特征相似性高则保留，否则过滤此锚点。
   • 检查锚点细胞所在数据集最邻近的30个细胞，查看它们重叠的锚点配对细胞的数量，重叠越多分值越高，代表锚点可靠性更高。

4. 经过层层过滤剩下的锚点细胞对，可以认为它们是相同类型和状态的细胞，它们之间的基因表达差异是技术偏倚引起的。Seurat计算它们的差异向量，然后用此向量校正这个锚点锚定的细胞子集的基因表达值。校正后的基因表达值即消除了技术偏倚，实现了两个单细胞数据集的整合。

Seurat工作流程图

Seurat工作流程

Seurat示例代码

pbmc.list <- SplitObject(pbmc, split.by="batch")
for(i in 1:length(pbmc.list)) {
  pbmc.list[[i]] <- NormalizeData(pbmc.list[[i]],verbose=FALSE)
  pbmc.list[[i]] <- FindVariableFeatures(pbmc.list[[i]],selection.method="vst",verbose=FALSE)
}
AllBatch.anchors <- FindIntegrationAnchors(object.list = pbmc.list, dims = 1:20,k.filter=80)
pbmc <- IntegrateData(anchorset = AllBatch.anchors, dims = 1:20)
DefaultAssay(pbmc) <- "integrated"
pbmc <- ScaleData(pbmc, verbose = T, do.center=F,split.by = 'batch')
pbmc <- RunPCA(pbmc, npcs = 50, verbose = T)
pbmc <- RunUMAP(pbmc, reduction = "pca", dims = 1:20,
                reduction.name = 'seurat_umap', reduction.key = 'seurat_umap')
DimPlot(pbmc, reduction='seurat_umap', group.by='batch', pt.size=0.3,label=F)

参考

Tutorial：https://satijalab.org/seurat/articles/integration_mapping.html
Paper：https://www.cell.com/cell/fulltext/S0092-8674(19)30559-8

LIGER

LIGER原理

• LIGER是从多个批次的数据中挑出一组共同的潜在因子，这些因子代表了每个数据集中相同的生物“信号”，但同时也保留了数据集的差异。为的是用这些因子联合识别细胞类型，以及识别和保留细胞类型特异性的差异。
• LIGER首先采用综合非负矩阵分解（iNMF）来学习低维空间，首先对每个批次的单细胞基因表达矩阵经过iNMF分解，可以得到三个矩阵：因子(行) x 细胞 (列)的H矩阵，数据集特有基因(行) x 因子 (列)的V矩阵，以及数据集共有基因(行) x 因子 (列)的W矩阵。

iNMF原理

• 我们知道NMF的优化目标是在迭代中优化以下误差

NMF优化目标

• 而当我们有多个数据集并想将其整合起来的时候，就需要使用iNMF。设有k个待整合的数据集，iNMF对NMF的形式进行如下修正：

iNMF的优化目标

• 我们可以假定待整合的不同数据集之间必然存在着共有的一些特征，也必然存在着一些独特的特征。因此我们可以把原始W矩阵写成W矩阵与Vk矩阵的叠加，新的W矩阵则相当于不同数据集之间共有的特征，而Vk矩阵则是第k个数据集独特的特征。另外引入λ正则化项，λ的大小反映了数据集之间的相似性。当λ等于0时，VH的二范数对目标函数没有影响，因此在优化的过程中我们得到的W矩阵（也就是共有的特征）会趋向于0，而数据集独特的特征则会使得Error项尽可能的小，因此数据集之间的差异（也就是数据集独特的特征）就会尽可能的大。同理，当λ趋于正无穷大时，减小VH的二范数对优化目标函数的权重极大，在这种情况下，W矩阵，也就是数据集共有的特征，就会尽可能的大。

• 因此在进行非负矩阵分解的过程中，我们需要选择合适的参数K和λ，注意，这里的K指的是W矩阵的维度。在最理想的状态下，数据集中的每一个真实存在的cluster将占据W矩阵的一个纬度，此时数据之间分散得越开，H矩阵的KL散度就越大。因此参数K的选择可以基于H矩阵的KL散度进行优化，当KL散度刚刚饱和时的K应当是最优的K。同理，一般而言，我们希望在进行数据整合之后，不同数据集之间能够尽可能均匀的map在一起，这种均匀的程度可以用Alignment score进行衡量，随着λ的增大，我们假定的数据集之间的相似性越高，通常Alignment score就会越高。因此参数λ的选择可以基于聚类的Alignment score进行优化，当Alignment score刚刚饱和时的λ应当是最优的λ。但必须要指出的是，这种参数的优化仅仅基于通常无先验知识的情况，我们更建议以此为参考，根据我们对数据的理解和先验知识选择合适的参数。

• iNMF分解矩阵后的“因子”，相当于PCA降维之后的"PC"轴；但是iNMF的因子比 PCA的PC轴可解释性更强，在单细胞数据分析中，每个因子常常对应一种细胞类型。V矩阵和W矩阵中loading值靠前的基因，往往具有显著的生物学意义。
• Liger的聚类算法也比较有意思，它并不基于单个细胞的坐标（label）进行聚类，而是用单个细胞周围k个邻居细胞的坐标（label）为该细胞赋予一个新的坐标，再根据曼哈顿距离或其他的范数距离进行聚类。Liger作者的观点是单个细胞被错误分配坐标（label）的概率要远大于其周围k个邻居细胞坐标都被系统性地错误分配的概率要高得多，因此使用周围k个邻居细胞的坐标来替代该细胞本身的坐标会使聚类的结果更加鲁棒。

LIGER工作流程图

LIGER工作流程

LIGER示例代码

library(liger)
library(SeuratWrappers)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc, split.by = "batch", do.center = FALSE)
pbmc <- RunOptimizeALS(pbmc, k = 20, lambda = 5, split.by = "batch")
pbmc <- RunQuantileNorm(pbmc, split.by = "batch")
pbmc <- RunUMAP(pbmc, dims = 1:ncol(pbmc[["iNMF"]]), reduction = "iNMF",
                reduction.name = 'liger_umap', reduction.key = 'liger_umap')
DimPlot(pbmc, group.by = "batch", label=F, reduction = 'liger_umap')

参考

Tutorial：https://htmlpreview.github.io/?https://github.com/satijalab/seurat.wrappers/blob/master/docs/liger.html
Paper：https://www.cell.com/cell/fulltext/S0092-8674%2819%2930504-5
https://www.plob.org/article/24538.html

Harmony

Harmony原理

• Harmony需要输入低维空间的坐标值（embedding），一般使用PCA的降维结果。Harmony导入PCA的降维数据后，会采用soft k-means clustering算法将细胞聚类。常用的聚类算法仅考虑细胞在低维空间的距离，但是soft clustering算法会考虑我们提供的校正因素。
• 例如，细胞c2距离cluster1有点远，本来不能算作cluster1的一份子；但是c2和cluster1的细胞来自不同的数据集，因为我们期望不同的数据集融合，所以破例让它加入cluster1了。
• 聚类之后，先计算每个cluster内各个数据集的细胞的中心点，然后根据这些中心点计算各个cluster的中心点。最后通过算法让cluster内的细胞向中心聚集，实在收敛不了的离群细胞就过滤掉。调整之后的数据重复：聚类—计算cluster中心点—收敛细胞—聚类的过程，不断迭代直至聚类效果趋于稳定。

Harmony工作流程图

harmony工作流程

Harmony示例代码

library(harmony)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc, split.by = "batch", do.center = FALSE)
pbmc <- RunPCA(pbmc, npcs = 50, verbose = FALSE)
pbmc <- RunHarmony(pbmc, group.by.vars = "batch")
pbmc <- RunUMAP(pbmc, reduction = "harmony", dims = 1:20,
                reduction.name = 'harmony_umap', reduction.key = 'harmony_umap')
DimPlot(pbmc, reduction='harmony_umap', label=F, group.by = 'batch')

参考

Tutorial：http://htmlpreview.github.io/?https://github.com/immunogenomics/harmony/blob/master/docs/SeuratV3.html
Paper：https://www.nature.com/articles/s41592-019-0619-0

BEER

BEER的主要步骤包含两步

1. 首先对每个批次的单细胞矩阵进行降维，对相似的10个细胞进行分组，产生一个代表性的表达谱，然后找出两个批次中所有相似的细胞组（Mutual Nearest, MN）。
2. 将两个批次整合，然后进行PCA。对于一对MN，如果没有批次效应，那么他们在某个PC上应该具有相似的值；BEER删除那些在一对MN具有低相关性的PC，以达到去除批次效应的目的。

BEER工作流程图

BEER工作流程

BEER示例代码

source("../BEER-master/BEER.R")
DATA <- pbmc@assays$RNA@counts
DATA.BATCH <- pbmc$batch
RM1=read.table('../data/Apoptosis.txt',sep='\t',header=FALSE)[,1]
RM2=read.table('../data/CellCycle.txt',sep='\t',header=FALSE)[,1]
RM3=read.table('../data/Ribosome.txt',sep='\t',header=FALSE)[,1]
RMG=c(as.character(RM1),as.character(RM2),as.character(RM3))
mybeer <- BEER(DATA, DATA.BATCH, GNUM=30, PCNUM=50, ROUND=1, GN=2000, SEED=1, COMBAT=TRUE, RMG=RMG)
mybeer <- ReBEER(mybeer, GNUM=30, PCNUM=50, ROUND=1, SEED=1, RMG=RMG) 
PCUSE <- mybeer$select
PCUSE <- PCUSE[which(PCUSE <= 20)]
# PCUSE=.selectUSE(mybeer, CUTR=0.8, CUTL=0.8, RR=0.5, RL=0.5)
mybeer$seurat@meta.data <- pbmc@meta.data
pbmc0 <- mybeer$seurat
pbmc0 <- RunUMAP(object = pbmc0, reduction='pca',dims = PCUSE, check_duplicates=FALSE,
                 reduction.name = "beer_umap",reduction.key = "beer_umap",
                 min.dist = 0.3, n.neighbors = 200L)
DimPlot(pbmc0, reduction='beer_umap', group.by='batch', pt.size=0.1) + NoLegend()

参考

Tutorial：https://github.com/jumphone/BEER
Paper：https://www.nature.com/articles/s41421-019-0114-x

参考

• https://mp.weixin.qq.com/s/XUm-1Xi4StA6N7sz9CA1IA
• https://mp.weixin.qq.com/s/txnz-wfBk5vLpBiuGBnNTQ
• https://mp.weixin.qq.com/s?__biz=MzIyMzMwNDQ2MA==&mid=2247484164&idx=1&sn=9da8c20591b95cc6663cea2c5cb195bd&chksm=e8210d67df568471612b6426a599e618d0158e734a0f8c3fafa4b1eaa2b7adaec373a30a2ee6&scene=178&cur_album_id=1593332494622359552#rd