单细胞分析十八般武艺3：fastMNN

生信技能树jimmy

发布于 2021-04-16 10:38:53

6.8K0

发布于 2021-04-16 10:38:53

文章被收录于专栏：单细胞天地

单细胞测序技术的发展日新月异，新的分析工具也层出不穷。每个工具都有它的优势与不足，在没有权威工具和流程的单细胞生信江湖里，多掌握几种分析方法和工具，探索数据时常常会有意想不到的惊喜。

fastMNN简介

MNN是Haghverdi等人提出的一种批次校正算法，由R语言的batchelor包实现。算法原理及性能评测的细节可参阅2018发表在Nature Biotechnology上的论文：Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. 此算法是Seurat3锚点整合算法的核心部分，也被Monocle3采纳作为批次校正的算法。fastMNN是MNN的升级版，主要改动是fastMNN采用PCA降维之后的低维空间计算细胞之间的距离，而MNN直接使用原始表达矩阵计算细胞之间的距离，因此分析速度会更快。

MNN的工作原理

先解释一下MNN pairs是什么意思：

首先根据两个批次样本的表达矩阵计算细胞之间的余弦距离（常用来表示单细胞转录特征相似性的度量值，可以在一定程度上屏蔽测序深度和捕获效率不同造成的差异）；
然后为Batch1样本中的细胞（记作Cell-i）在Batch2样本寻找k个余弦距离最近的细胞集（记作Set-i）；
同样为Batch2样本中的细胞（记作Cell-j）在Batch1样本寻找k个距离最近的细胞集（记作Set-j）；
如果Cell-i在Set-j中，恰好Cell-j也在Set-i中，那么Cell-i与Cell-j就是相互最近邻细胞对（即MNN pair）。

（a） Batch1样本和Batch2样本的细胞转录谱具有与生物学意义无关的差异，即批次效应产生的差异。（b）MNN算法通过寻找相互最近邻的细胞对来识别两个批次样本中相同的细胞类型。（c）在MNN pairs之间计算批次修正向量。（d） Batch1作为参考，Batch2通过减去校正向量整合到Batch1中。（e）整合后的数据作为参考，任何新批次的整合都会重复该程序。

R包的安装

fastMNN算法由batchelor包提供，安装seurat包时已经自动安装过此包。要单独使用fastMNN算法整合数据，需要使用SeuratWrappers包，其安装命令如下：

remotes::install_github('satijalab/seurat-wrappers')

10X单细胞数据实测

数据依然使用《单细胞转录组高级分析一：多样本合并与批次校正》一文中使用的那10个样本的数据，没有数据的朋友可以添加我的微信后索取，微信二维码可以点击文末“阅读原文”找到。

数据准备

library(Seurat)
library(tidyverse)
rm(list=ls())
setwd("~/project/harmony")
dir <- dir("GSE139324/")           #GSE139324是存放数据的目录
dir <- paste0("GSE139324/",dir)
sample_name <- c('HNC01PBMC', 'HNC01TIL', 'HNC10PBMC', 'HNC10TIL', 'HNC20PBMC',
                 'HNC20TIL',  'PBMC1', 'PBMC2', 'Tonsil1', 'Tonsil2')
scRNAlist <- list()
for(i in 1:length(dir)){
counts <- Read10X(data.dir = dir[i])
scRNAlist[[i]] <- CreateSeuratObject(counts, project=sample_name[i], min.cells=3, min.features = 200)
scRNAlist[[i]][["percent.mt"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^MT-")
scRNAlist[[i]] <- subset(scRNAlist[[i]], subset = percent.mt < 10) 
}   
saveRDS(scRNAlist, "scRNAlist.rds")

数据整合

library(Seurat)
library(tidyverse)
library(patchwork)
library(SeuratWrappers)
scRNAlist <- readRDS("scRNAlist.rds")
##为了更好展现整合效果，只使用两个样本整合
scRNAlist <- list(scRNAlist[[2]], scRNAlist[[10]])
scRNAlist <- lapply(scRNAlist, FUN = function(x) NormalizeData(x))
scRNAlist <- lapply(scRNAlist, FUN = function(x) FindVariableFeatures(x))
scRNA <- RunFastMNN(object.list = scRNAlist)
scRNA <- RunUMAP(scRNA, reduction = "mnn", dims = 1:30)
scRNA <- FindNeighbors(scRNA, reduction = "mnn", dims = 1:30)
scRNA <- FindClusters(scRNA)
p1 <- DimPlot(scRNA, group.by = "orig.ident", pt.size=0.1) + 
      ggtitle("Integrated by fastMNN")
p2 <- DimPlot(scRNA.orig, group.by="orig.ident", pt.size=0.1) + 
      ggtitle("No integrated")
p = p1 + p2 + plot_layout(guides='collect')
ggsave('fastMNN.png', p, width=8, height=4)