seurat3的merge功能和cellranger的aggr整合多个10X单细胞转录组对比

昨天我在单细胞天地的教程：使用seurat3的merge功能整合8个10X单细胞转录组样本 完完整整的展示了如何使用seurat3的merge功能整合8个10X单细胞转录组样本，因为这个数据集的文章作者使用的是cellranger流程，而且我们在单细胞天地多次分享过流程笔记，大家可以自行前往学习，如下：

• 单细胞实战(一)数据下载

• 单细胞实战(二) cell ranger使用前注意事项

• 单细胞实战(三) Cell Ranger使用初探

• 单细胞实战(四) Cell Ranger流程概览

• 单细胞实战(五) 理解cellranger count的结果

我们得比较一下，作者的ellranger的aggr整合多个10X单细胞转录组得到的表达矩阵，跟我们使用seurat3的merge功能整合8个10X单细胞转录组样本后的表达矩阵是否有差异。

幸运的是，作者在GEO数据库也放上去了自己ellranger的aggr整合多个10X单细胞转录组得到的表达矩阵，所以我们只需要下载走seurat3的降维分群即可跟昨天的教程：使用seurat3的merge功能整合8个10X单细胞转录组样本 聚类分群结果进行对比啦！

首先走seurat3的降维分群流程

library(Seurat)
sce=CreateSeuratObject(counts = Read10X('input/'), 
                     project = 'merge' )

# 步骤 ScaleData 的耗时取决于电脑系统配置（保守估计大于一分钟）
sce <- ScaleData(object = sce, 
                     vars.to.regress = c('nCount_RNA'), 
                     model.use = 'linear', 
                     use.umi = FALSE)
sce <- FindVariableFeatures(object = sce, 
                                mean.function = ExpMean, 
                                dispersion.function = LogVMR, 
                                x.low.cutoff = 0.0125, 
                                x.high.cutoff = 4, 
                                y.cutoff = 0.5)
length(VariableFeatures(sce)) 
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))
# 下面只是展现不同降维算法而已，并不要求都使用一次。
sce <- RunICA(sce )
sce <- RunTSNE(sce )
#sce <- RunUMAP(sce,dims = 1:10)
#VizPCA( sce, pcs.use = 1:2)
DimPlot(object = sce, reduction = "pca") 
DimPlot(object = sce, reduction = "ica")
DimPlot(object = sce, reduction = "tsne")
#DimPlot(object = sce, reduction = "umap")

# 针对PCA降维后的表达矩阵进行聚类 FindNeighbors+FindClusters 两个步骤。
sce <- FindNeighbors(object = sce, dims = 1:20, verbose = FALSE) 
sce <- FindClusters(object = sce, resolution = 0.3,verbose = FALSE)
# 继续tSNE可视化
DimPlot(object = sce, reduction = "tsne",group.by = 'RNA_snn_res.0.3')
DimPlot(object = sce, reduction = "tsne",
        group.by = 'orig.ident')
head(rownames(sce@meta.data))
tail(rownames(sce@meta.data))
library(stringr)
table(str_split(rownames(sce@meta.data),'-',simplify = T)[,2],
      sce@meta.data$RNA_snn_res.0.3)

可以看到4个发育时期的8个样本，仍然是同一个发育时期的2个样本非常类似，不同发育时期的样本就差异很大。所以也可以理解为10X单细胞转录组的样本批次效应被去除了。

我们下一讲会详细说明，其实大家也可以自己看看试试看更多数据整合，去除批次效应算法，比如

• Batch correction: integrative non-negative matrix factorization (NMF) using LIGER  

• Batch correction: canonical correlation analysis (CCA) using Seurat

• 三个10X单细胞转录组样本CCA整合

这本书基本上都介绍了：https://broadinstitute.github.io/2019_scWorkshop/correcting-batch-effects.html

仔细检查两次分类

前面我们虽然得到的结论都是同一个发育时期的2个样本非常类似，不同发育时期的样本就差异很大。实际上，我们有着8个样本的详细信息，是可以深入比较每个细胞的：

写在后面

我的教程承接生信技能树的一贯优良传统，都是可以完完整整的复现的，而且是免费提供给所有人看，在国庆节假期，花了4个小时，虽然经过统计，公众号粉丝只有不到3%的人会打赏，不过我本来就不靠文字吃饭，走的是全名偶像路线。

虽然教程是免费，但并不意味着读者有疑问我就会为他/她写教程， 这样不合理，这个教程是写给购买了单细胞课程的学员，教程配套的数据代码，都是会上传到售后微云给到学员！

如果你也想享受这样的VIP待遇，请自助购买，然后找小助手拉你进入答疑群：

• 手机微信观看（推荐）：http://jm.grazy.cn/wechat/course/detail?cid=55  

• 电脑浏览器的链接 http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55  

点击阅读原文，就可以自助购买啦，请不要咨询，一旦购买，概不退款，直接下单就好，不要磨磨唧唧讨人厌！