0
0

fastMNN|手把手教你理解和实现单细胞批次效应校正方法

再不更新就要忘记学习了，继续坚持分享，今天分享的是单细胞批次效应校正方法的fastMNN算法。

MNN原理介绍

先从MNN讲起，MNN是Haghverdi等人提出的一种批次校正算法，由R语言的batchelor包实现。算法原理及性能评测可自行阅读论文：Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors。大部分批次效应校正方法错误地假设细胞群的组成在各个批次之间是已知的或相同的。而MNN方法不依赖于各个批次之间预定义或相等的群体组成；相反，它只要求各个批次之间共享一个群体子集。

MNN使用假设:（i）至少有一个细胞群同时存在于两个批次中,（ii）批次效应几乎与生物子空间正交,（iii）批次效应变化远小于不同细胞类型之间的生物效应变化。

MNN pair：

首先根据两个批次样本的表达矩阵计算细胞之间的余弦距离（常用来表示单细胞转录特征相似性的度量值，可以在一定程度上屏蔽测序深度和捕获效率不同造成的差异）
然后为Batch1样本中的细胞（记作Cell-i）在Batch2样本寻找k个余弦距离最近的细胞集（记作Cell-set-i）
同样为Batch2样本中的细胞（记作Cell-j）在Batch1样本寻找k个距离最近的细胞集（记作Cell-set-j）
如果Cell-i在Cell-Set-j中，恰好Cell-j也在Cell-Set-i中，那么Cell-i与Cell-j就是相互最近邻细胞对（即MNN pair）
MNN pair 中细胞之间表达值的差异提供了批次效应的估计值，通过对许多这样的MNN pair对取平均值可以使其更加精确。从估计的批次效应中获得校正向量，并将其应用于表达值以执行批次校正
在这里插入图片描述
fastMNN原理介绍

fastMNN是MNN的升级版，主要改动是fastMNN采用PCA降维之后的低维空间计算细胞之间的距离，而MNN直接使用原始表达矩阵计算细胞之间的距离，因此分析速度会更快。

fastMNN简述：

对（余弦）标准化表达值执行多样本 PCA 以降低维数。

在参考批次和目标批次之间的低维空间中识别 MNN 对。

消除参考批次和目标批次中沿平均批次向量的变化。

使用局部加权校正向量，将目标批次中的细胞校正至参考值。

将更正后的目标批次与参考批次合并，然后重复下一个目标批次。

fastMNN实现：

实现方法源自：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1850-9#Sec13。

数据：dataset4，胰腺细胞数据

代码复现：
1.使用batchelor包进行mnn，根据文献中的fastMNN进行分析

  1. library(Seurat)
  2. library(batchelor)
  3. expr_mat_mnn <- readRDS(file ="C:/Users/Wang/Documents/单细胞学习/单细胞蛋白组学/代码/MNNanalysis/data/pancreas_v3_files/pancreas_expression_matrix.rds")
  4. metadata <- readRDS(file ="C:/Users/Wang/Documents/单细胞学习/单细胞蛋白组学/代码/MNNanalysis/data/pancreas_v3_files/pancreas_metadata.rds")######数据前处理
  5. b_seurat <- CreateSeuratObject(expr_mat_mnn, meta.data = metadata, project ="MNNcorrect_benchmark")
  7. normData = T
  8. norm_method ="LogNormalize"
  9. scale_factor =10000
  10. batch_label ="tech"
  11. celltype_label ="celltype"#####对数据进行归一化if(normData){
  12.   b_seurat <- NormalizeData(object = b_seurat, normalization.method = norm_method, scale.factor = scale_factor)}else{
  13.   b_seurat@data = b_seurat@raw.data
  14. }
  15. b_seurat <- ScaleData(object = b_seurat)
  16. b_seurat <- FindVariableFeatures(object =  b_seurat)
  17. b_seurat <- ScaleData(object =  b_seurat)
  18. b_seurat <- RunPCA(object =  b_seurat)
  19. b_seurat <- FindNeighbors(object =  b_seurat, dims =1:30)
  20. b_seurat <- FindClusters(object =  b_seurat)
  21. b_seurat <- RunUMAP(object =  b_seurat, dims =1:30)
  22. DimPlot(object = b_seurat, reduction ="umap")
  23. meta_data <- b_seurat@meta.data
  24. data_filtered<-  b_seurat@assays$RNA$scale.data
  25. #split data into batches
  26. data_list <- list()
  27. meta_list <- list()
  28. list_batch <- unique(meta_data[,batch_label])for(i in1:length(list_batch)){
  29.   selected <- row.names(meta_data)[which (meta_data[, batch_label]== list_batch[i])]
  30.   data_list[[i]]<- data_filtered[,selected]
  31.   meta_list[[i]]<- meta_data[selected,]}  
  33.   b1_data <- data_list[[1]]
  34.   b2_data <- data_list[[2]]
  35.  
  36.   b1_meta <- meta_list[[1]]
  37.   b2_meta <- meta_list[[2]]
  38.   batch_data<- list(b1_data, b2_data, b1_meta, b2_meta)#####函数实现####
  39. call_mnn(data_list,meta_list, batch_label, celltype_label)#####call_mnn函数
  40. call_mnn <-function(batch_data, batch_label, celltype_label, npcs =20, 
  41.                          plotout_dir ="", saveout_dir ="", 
  42.                          outfilename_prefix ="", 
  43.                          visualize = T, save_obj = T){#plotting
  44.   k_seed =10
  45.   tsne_perplex =30
  46.   tsneplot_filename ="_fastmnn_tsne"
  47.   mnn_out_filename ="_fastmnn_out"
  48.   metadata_out_filename ="_fastmnn_metadata_out"
  49.   pca_filename ="_fastmnn_pca"
  51.   b1_data <- batch_data[[1]]
  52.   b2_data <- batch_data[[2]]
  53.   b1_meta <- batch_data[[3]]
  54.   b2_meta <- batch_data[[4]]
  55.   metadata <- rbind(b1_meta, b2_meta)##############################################run
  56.   t1 = Sys.time()
  57.   cosd1 <- cosineNorm(as.matrix(b1_data))
  58.   cosd2 <- cosineNorm(as.matrix(b2_data))
  60.   pcs_total <- multiBatchPCA(cosd1, cosd2)#defaut npca = 50
  61.   out_mnn_total <- fastMNN(pcs_total[[1]], pcs_total[[2]], pc.input =TRUE)
  62.   t2 = Sys.time()
  63.   print(t2-t1)######################################save
  65.   saveRDS(out_mnn_total, file=paste0(saveout_dir, outfilename_prefix, mnn_out_filename,".RDS"))
  66.   corrected_pca_mnn <- as.data.frame(out_mnn_total$corrected)
  67.   rownames(corrected_pca_mnn)<- colnames(cbind(b1_data,b2_data))
  68.   cells_use <- rownames(corrected_pca_mnn)
  69.   metadata <- metadata[cells_use,]
  71.   corrected_pca_mnn[rownames(metadata), batch_label]<- metadata[, batch_label]#corrected_pca_mnn[rownames(metadata), 'batch'] <- metadata[ , "batch"]
  72.   corrected_pca_mnn[rownames(metadata), celltype_label]<- metadata[, celltype_label]
  73.   write.table(corrected_pca_mnn, file=paste0(saveout_dir, outfilename_prefix, pca_filename,".txt"), quote=F, sep='\t', row.names = T, col.names =NA)
  75.   saveRDS(metadata, file = paste0(saveout_dir, outfilename_prefix, metadata_out_filename,".RDS"))#write.table(metadata, file=paste0(saveout_dir, outfilename_prefix, metadata_out_filename, ".txt"), quote=F, sep="\t", row.name=TRUE, col.name=NA)if(visualize){
  77.     set.seed(10)
  78.     out_tsne <- Rtsne(out_mnn_total$corrected, check_duplicates=FALSE, perplexity =30)
  79.     tsne_df<- as.data.frame(out_tsne$Y)
  81.     rownames(tsne_df)<- colnames(cbind(b1_data,b2_data))
  82.     dim(tsne_df)
  83.     colnames(tsne_df)<- c('tSNE_1','tSNE_2')
  84.     tsne_df$batchlb <- metadata[, batch_label]#tsne_df$batch <- metadata[ , "batch"]
  85.     tsne_df$CellType <- metadata[, celltype_label]
  87.     p01 <- ggplot(tsne_df, aes(x = tSNE_1, y = tSNE_2, colour = batchlb))+ geom_point(alpha =0.6)+ theme_bw()
  88.     p01 <- p01 + labs(x='tSNE_1',y='tSNE_2',title='fast MNN')
  89.     p01 <- p01 + theme(legend.title = element_text(size=17), 
  90.                        legend.key.size = unit(1.1,"cm"),
  91.                        legend.key.width = unit(0.5,"cm"), 
  92.                        legend.text = element_text(size=14), 
  93.                        plot.title = element_text(color="black", size=20, hjust =0.5))
  95.     p02 <- ggplot(tsne_df, aes(x = tSNE_1, y = tSNE_2, colour = CellType))+ geom_point(alpha =0.6)+ theme_bw()
  96.     p02 <- p02 + labs(x='tSNE_1',y='tSNE_2',title='fast MNN')
  97.     p02 <- p02 + theme(legend.title = element_text(size=17), 
  98.                        legend.key.size = unit(1.1,"cm"),
  99.                        legend.key.width = unit(0.5,"cm"), 
  100.                        legend.text = element_text(size=14), 
  101.                        plot.title = element_text(color="black", size=20, hjust =0.5))
  103.     png(paste0(plotout_dir,outfilename_prefix,tsneplot_filename,".png"),width =2*1000, height =800, res =2*72)
  104.     print(plot_grid(p01, p02))
  105.     dev.off()
  107.     pdf(paste0(plotout_dir,outfilename_prefix,tsneplot_filename,".pdf"),width=15,height=7,paper='special')
  108.     print(plot_grid(p01, p02))
  109.     dev.off()}}
  1. 使用seurat-wrappers包实现fastMNN
  1. rm(list=ls())
  2. remotes::install_github('satijalab/seurat-wrappers')
  3. library(Seurat)
  4. library(tidyverse)
  5. library(Seurat)
  6. library(tidyverse)
  7. library(patchwork)
  8. library(SeuratWrappers)
  9. dir <- dir("GSE139324/")#GSE139324是存放数据的目录，自行下载进行分析
  10. dir <- paste0("GSE139324/",dir)
  11. sample_name <- c('HNC01PBMC','HNC01TIL','HNC10PBMC','HNC10TIL','HNC20PBMC','HNC20TIL','PBMC1','PBMC2','Tonsil1','Tonsil2')
  12. scRNAlist <- list()for(i in1:length(dir)){
  13. counts <- Read10X(data.dir = dir[i])
  14. scRNAlist[[i]]<- CreateSeuratObject(counts, project=sample_name[i], min.cells=3, min.features =200)
  15. scRNAlist[[i]][["percent.mt"]]<- PercentageFeatureSet(scRNAlist[[i]], pattern ="^MT-")
  16. scRNAlist[[i]]<- subset(scRNAlist[[i]], subset = percent.mt <10)}   
  17. saveRDS(scRNAlist,"scRNAlist.rds")
  18. scRNAlist <- readRDS("scRNAlist.rds")##为了更好展现整合效果，只使用两个样本整合
  19. scRNAlist <- list(scRNAlist[[2]], scRNAlist[[10]])
  20. scRNAlist <- lapply(scRNAlist, FUN =function(x) NormalizeData(x))
  21. scRNAlist <- lapply(scRNAlist, FUN =function(x) FindVariableFeatures(x))
  22. scRNA <- RunFastMNN(object.list = scRNAlist)
  23. scRNA <- RunUMAP(scRNA, reduction ="mnn", dims =1:30)
  24. scRNA <- FindNeighbors(scRNA, reduction ="mnn", dims =1:30)
  25. scRNA <- FindClusters(scRNA)
  26. p1 <- DimPlot(scRNA, group.by ="orig.ident", pt.size=0.1)   
  27.       ggtitle("Integrated by fastMNN")
  28. p2 <- DimPlot(scRNA.orig, group.by="orig.ident", pt.size=0.1)   
  29.       ggtitle("No integrated")
  30. p = p1   p2   plot_layout(guides='collect')
  31. ggsave('fastMNN.png', p, width=8, height=4)

欢迎关注转发与交流，如有疑问请联系管理员：kriswcyYQ进行相关问题解决。

标签： 数据分析 数据挖掘 python
本文转载自: https://blog.csdn.net/wcy1995427/article/details/140159000
版权归原作者 王元启的生信记录 所有， 如有侵权，请联系我们删除。
登录后发布评论

“fastMNN|手把手教你理解和实现单细胞批次效应校正方法”的评论:

还没有评论
关于作者
...
overfit同步小助手
文章同步
相关阅读

网络安全法-网络运行安全

使用selenium/drissionpage时如何阻止chrome自动跳转http到https

docker desktop 里部署的Open WebUI 管理员密码忘记了的处理方法

在ubuntu20.04中搭建onsite比赛运行环境

利用开源的低代码表单设计器FcDesigner高效管理和渲染复杂表单结构

Kafka学习笔记

【前端】浏览器输入url到页面呈现发生了什么？

文章导航
没有找到文章内容分段。。。