查看原文
其他

手把手教你做单细胞测序(四)——多样本整合

Biomamba Biomamba 生信基地 2024-04-03

上次的视频中已花费大量时间讲解过单样本分析的基本流程,所以这节课的学习需要有上节课的基础,希望大家按顺序观看。此次的内容较简单、篇幅也较小,代码与视频请看下文,测试数据代码存于文末,可以投微信豆或回复”多样本整合“获取。由于测试数据比较特殊,并没有展示出去批次的精妙之处,留一个悬念给大家吧,可以用自己的数据集测试一下。


视频:


代码:


###########单纯的merge################# library(Seurat) library(multtest) library(dplyr) library(ggplot2) library(patchwork)  ##########准备用于拆分的数据集###########pbmc <- subset(pbmc, downsample = 50)ifnb <- readRDS('pbmcrenamed.rds')ifnb.list <- SplitObject(ifnb, split.by = "group")C57 <- ifnb.list$C57AS1 <- ifnb.list$AS1######简单merge######## #不具有去批次效应功能pbmc <- merge(C57, y = c(AS1), add.cell.ids = c("C57", "AS1"), project = "ALL")pbmchead(colnames(pbmc))unique(sapply(X = strsplit(colnames(pbmc), split = "_"), FUN = "[", 1))table(pbmc$orig.ident)##############anchor###############library(Seurat)library(tidyverse)### testA ----myfunction1 <- function(testA.seu){ testA.seu <- NormalizeData(testA.seu, normalization.method = "LogNormalize", scale.factor = 10000) testA.seu <- FindVariableFeatures(testA.seu, selection.method = "vst", nfeatures = 2000) return(testA.seu)}C57 <- myfunction1(C57)AS1 <- myfunction1(AS1)
### Integration ----testAB.anchors <- FindIntegrationAnchors(object.list = list(C57,AS1), dims = 1:20)testAB.integrated <- IntegrateData(anchorset = testAB.anchors, dims = 1:20)
#需要注意的是:上面的整合步骤相对于harmony整合方法,对于较大的数据集(几万个细胞)#非常消耗内存和时间,大约9G的数据32G的内存就已经无法运行;#当存在某一个Seurat对象细胞数很少(印象中200以下这样子),#会报错,这时建议用第二种整合方法
DefaultAssay(testAB.integrated) <- "integrated"
# # Run the standard workflow for visualization and clusteringtestAB.integrated <- ScaleData(testAB.integrated, features = rownames(testAB.integrated))testAB.integrated <- RunPCA(testAB.integrated, npcs = 50, verbose = FALSE)testAB.integrated <- FindNeighbors(testAB.integrated, dims = 1:30)testAB.integrated <- FindClusters(testAB.integrated, resolution = 0.5)testAB.integrated <- RunUMAP(testAB.integrated, dims = 1:30)testAB.integrated <- RunTSNE(testAB.integrated, dims = 1:30)p1<- DimPlot(testAB.integrated,label = T,split.by = 'group')#integrated
DefaultAssay(testAB.integrated) <- "RNA"testAB.integrated <- ScaleData(testAB.integrated, features = rownames(testAB.integrated))testAB.integrated <- RunPCA(testAB.integrated, npcs = 50, verbose = FALSE)testAB.integrated <- FindNeighbors(testAB.integrated, dims = 1:30)testAB.integrated <- FindClusters(testAB.integrated, resolution = 0.5)testAB.integrated <- RunUMAP(testAB.integrated, dims = 1:30)testAB.integrated <- RunTSNE(testAB.integrated, dims = 1:30)
p2 <- DimPlot(testAB.integrated,label = T,split.by = 'group')p1|p2
###########harmony 速度快、内存少################if(!require(harmony))devtools::install_github("immunogenomics/harmony")test.seu <- pbmctest.seu <- test.seu%>% Seurat::NormalizeData() %>% FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>% ScaleData()test.seu <- RunPCA(test.seu, npcs = 50, verbose = FALSE)

#####run 到PCA再进行harmony,相当于降维########test.seu=test.seu %>% RunHarmony("group", plot_convergence = TRUE)
test.seu <- test.seu %>% RunUMAP(reduction = "harmony", dims = 1:30) %>% FindNeighbors(reduction = "harmony", dims = 1:30) %>% FindClusters(resolution = 0.5) %>% identity()
test.seu <- test.seu %>% RunTSNE(reduction = "harmony", dims = 1:30) p3 <- DimPlot(test.seu, reduction = "tsne", group.by = "group", pt.size=0.5)+theme( axis.line = element_blank(), axis.ticks = element_blank(),axis.text = element_blank())p4 <- DimPlot(test.seu, reduction = "tsne", group.by = "ident", pt.size=0.5, label = TRUE,repel = TRUE)+theme( axis.line = element_blank(), axis.ticks = element_blank(),axis.text = element_blank())p3|p4





往期回顾


单细胞数据分析系列教程:

B站视频,先看一遍视频再去看推送操作,建议至少看三遍:https://www.bilibili.com/video/BV1S44y1b76Z/
单细胞测序基础数据分析保姆级教程,代码部分整理在往期推送之中:手把手教你做单细胞测序数据分析(一)——绪论
手把手教你做单细胞测序数据分析(二)——各类输入文件读取
手把手教你做单细胞测序数据分析(三)——单样本分析



如何联系我们


公众号后台中消息更新不及时,超过48h后便不允许回复读者消息,这里给大家留一下领取资料、咨询服务器的微信号,方便大家随时交流、提建议(由助理接待)。此外呼声一直很高的交流群也建好了,欢迎大家入群讨论:你们要的微信、交流群,来了!
大家可以阅读完这几篇之后添加笑一笑也就算了有关提高"Biomamba 生信基地"运行效率事宜


代码与测试文件链接:

链接:

继续滑动看下一个
向上滑动看下一个

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存