# 5.5.1 使用Seurat的merge功能进行整合 ### 前言 单细胞数据未来会朝着多样本发展，因此数据整合是一项必备技能。cellranger中自带了aggr的整合功能，而这篇文章（Differentiation dynamics of mammary epithelial cells revealed by single-cell RNA-sequencing）的作者也正是这么做得到的组合后的表达矩阵，然后用`Read10X`读入 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-055124.png) **关于文章** 这是发表在2017年10月的NC文章。 作者的论点是：乳腺上皮细胞对研究乳腺癌的发展很重要，但目前只有很少的marker可以追踪这群细胞，因此有必要探索乳腺发育不同阶段的乳腺上皮细胞变化。 实验涉及了四个时期：8 weeks virgin =》nulliparous (NP) 未怀孕时期、14.5d gestation (G) 妊娠期第14.5天、6d lactation (L) 哺乳期第6天、11d post involution (PI) 完全退化第11天。其中每个时期都采集两只老鼠的组织细胞，所以一共8个样本，然后使用10X建库，那么最后的测序文件就是`8*3 = 24`个 如果要对这24个文件分别去整合，使用seurat的`merge`函数即可，不过**问题的关键是：如何用代码将这些样本区分开，然后分别构建对象，最后merge这些对象** ### 开始操作 ### **第一步：准备原始测序数据** 我们下载第一个：GSE106273\_RAW\.tar（183.5 Mb） > 感觉手机下载速度就是比电脑快，这个文件在手机上下载3分钟，电脑预计时间1小时 下载后解压，整个过程直接在Rstudio中的Terminal直接完成 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-062817.png) ### **第二步：整理数据** 思路：根据中间的分组信息（NP、G）将包含的文件放到相应的文件夹中 **方法一：shell脚本** ``` # 将同一组数据放在同一目录下 ls GSM* | awk -F '_' '{print $2"_"$3}'| uniq | while read i;do mkdir $i;mv *$i*gz $i;done # 各自重命名 find -name "*barcodes.tsv.gz" | while read i;do mv $i $(dirname $i)/barcodes.tsv.gz;done find -name "*genes.tsv.gz" | while read i;do mv $i $(dirname $i)/genes.tsv.gz;done find -name "*matrix.mtx.gz" | while read i;do mv $i $(dirname $i)/matrix.mtx.gz;done ``` **方法二：R脚本** ```r # 列出当前目录下所有开头是GSM的文件 fs=list.files('./','^GSM') # 然后获取四个样本信息 library(stringr) samples=str_split(fs,'_',simplify = T)[,1] # 设置一个循环，对每个样本信息做同样的事： #（1）找到包含这个样本的文件(用grepl) # (2)设置对应的目录名（str_split+paste）然后创建目录（用dir.create） # (3)将文件放到对应目录(采用的是file.rename)并重命名文件 lapply(unique(samples),function(x){ y=fs[grepl(x,fs)] folder=paste(str_split(y[1],'_',simplify = T)[,2:3], collapse = '') dir.create(folder,recursive = T) file.rename(y[1],file.path(folder,"barcodes.tsv.gz")) file.rename(y[2],file.path(folder,"genes.tsv.gz")) file.rename(y[3],file.path(folder,"matrix.mtx.gz")) }) ``` 需要注意的是，`Read10X`函数需要读取解压后的文件，于是还要对所有的数据文件进行解压 ``` find ./ -name "*gz" |xargs gunzip ``` > **常见错误：** > > * 说找不到Barcode文件，但明明存在Barcode： > > ```r > Error in Read10X() : Barcode file missing > ``` > > 那很有可能是因为三个10X数据的命名出了问题，一定要命名成"barcodes.tsv" "genes.tsv""matrix.mtx" > > 【补充： > > **cellranger的V2版本**得到的结果分别是：barcodes.tsv、genes.tsv、matrix.mtx； > > **V3版本**得到的结果分别是：matrix.mtx.gz、features.tsv.gz、barcodes.tsv.gz】 > > * 说找不到基因文件，那么就要看看测序数据是不是解压后的 > > ```r > Error in Read10X() : Gene name or features file missing > ``` ### **第三步：批量读取成10X对象** **Read10X() + CreateSeuratObject()** ```r # 因为Read10X函数需要对目录进行操作，所以先把目录名提取出来 folders=list.files('./',pattern='[12]$') > folders [1] "G_1" "G_2" "L_1" "L_2" [5] "NP_1" "NP_2" "PI_1" "PI_2" # 然后使用lapply进行循环(看下lapply的帮助文档就知道，它是对列表或向量进行循环，而apply是对数据框或矩阵操作) library(Seurat) sceList = lapply(folders,function(folder){ CreateSeuratObject(counts = Read10X(folder), project = folder ) }) # 此时的sceList仅仅是一个堆砌了8个10X对象的集合，下一步就要真正合并起来 > sceList [[1]] An object of class Seurat 27998 features across 2915 samples within 1 assay Active assay: RNA (27998 features) [[2]] An object of class Seurat 27998 features across 3106 samples within 1 assay Active assay: RNA (27998 features) ``` ### **第四步：组合** ```r sce.big <- merge(sceList[[1]], y = c(sceList[[2]],sceList[[3]],sceList[[4]], sceList[[5]],sceList[[6]], sceList[[7]],sceList[[8]]), add.cell.ids = folders, project = "mouse8") > table(sce.big$orig.ident) G_1 G_2 L_1 L_2 NP_1 NP_2 PI_1 PI_2 2915 3106 5906 3697 2249 2127 1500 4306 save(sce.big,file = 'sce.big.merge.mouse8.Rdata') # 保存的数据是1.4G ``` ### 补充 > 官网的merge教程在： 描述了三种情况 **第一种： merge两个seurat对象（原始数据）** 需要注意的是，组合数据时需要注明每个数据的名称，使用`add.cell.ids`参数指定 ```r pbmc.combined <- merge(pbmc4k, y = pbmc8k, add.cell.ids = c("4K", "8K"), project = "PBMC12K") pbmc.combined ## An object of class Seurat ## 33694 features across 12721 samples within 1 assay ## Active assay: RNA (33694 features) # 之后的组合数据就会出现列名的标识 head(colnames(pbmc.combined)) ## [1] "4K_AAACCTGAGAAGGCCT" "4K_AAACCTGAGACAGACC" "4K_AAACCTGAGATAGTCA" ## [4] "4K_AAACCTGAGCGCCTCA" "4K_AAACCTGAGGCATGGT" "4K_AAACCTGCAAGGTTCT" table(pbmc.combined$orig.ident) ## ## PBMC4K PBMC8K ## 4340 8381 ``` **第二种：merge两个以上（原始数据）** 将参数`y` 设成一个向量，就可以指定其他的数据 ```r pbmc.big <- merge(pbmc3k, y = c(pbmc4k, pbmc8k), add.cell.ids = c("3K", "4K", "8K"), project = "PBMC15K") unique(sapply(X = strsplit(colnames(pbmc.big), split = "_"), FUN = "[", 1)) ## [1] "3K" "4K" "8K" table(pbmc.big$orig.ident) ## ## pbmc3k PBMC4K PBMC8K ## 2638 4340 8381 ``` **第三种：merge归一化、标准化数据** 默认情况，只会组合原始数据，但如果有的数据时标准化之后的呢？ 其实可以通过一个参数`merge.data = TRUE`指定 ```r pbmc4k <- NormalizeData(pbmc4k) pbmc8k <- NormalizeData(pbmc8k) pbmc.normalized <- merge(pbmc4k, y = pbmc8k, add.cell.ids = c("4K", "8K"), project = "PBMC12K", merge.data = TRUE) ``` 看看第一种组合raw data和第三种组合normalized data对比： ```r ################# # raw data ################# GetAssayData(pbmc.combined)[1:10, 1:15] ## 10 x 15 sparse Matrix of class "dgCMatrix" ## ## RP11-34P13.3 . . . . . . . . . . . . . . . ## FAM138A . . . . . . . . . . . . . . . ## OR4F5 . . . . . . . . . . . . . . . ## RP11-34P13.7 . . . . . . . . . . . . . . . ## RP11-34P13.8 . . . . . . . . . . . . . . . ## RP11-34P13.14 . . . . . . . . . . . . . . . ## RP11-34P13.9 . . . . . . . . . . . . . . . ## FO538757.3 . . . . . . . . . . . . . . . ## FO538757.2 . . . . . . . . . 1 . . . . . ## AP006222.2 . . . . . . . . . . . 1 . . . ################# # normalized data ################# GetAssayData(pbmc.normalized)[1:10, 1:15] ## 10 x 15 sparse Matrix of class "dgCMatrix" ## ## RP11-34P13.3 . . . . . . . . . . . . . . . ## FAM138A . . . . . . . . . . . . . . . ## OR4F5 . . . . . . . . . . . . . . . ## RP11-34P13.7 . . . . . . . . . . . . . . . ## RP11-34P13.8 . . . . . . . . . . . . . . . ## RP11-34P13.14 . . . . . . . . . . . . . . . ## RP11-34P13.9 . . . . . . . . . . . . . . . ## FO538757.3 . . . . . . . . . . . . . . . ## FO538757.2 . . . . . . . . . 0.7721503 . . . . . ## AP006222.2 . . . . . . . . . . . 1.087928 . . . ``` > 上面的组合多个数据就结束了，接下来是检查组合后的分群结果 **首先检查原样本分群结果** ```r # 归一化+标准化（移除了不想要的差异来源nCount_RNA） sce.big <- NormalizeData(sce.big) sce.big <- ScaleData(sce.big, vars.to.regress = c('nCount_RNA'), model.use = 'linear', use.umi = FALSE) # 默认选2000个HVGs sce.big <- FindVariableFeatures(object = sce.big, mean.function = ExpMean, dispersion.function = LogVMR, mean.cutoff = c(0.0125,4), dispersion.cutoff = c(0.5,Inf)) # 降维（PCA+tSNE） sce.big <- RunPCA(object = sce.big, pc.genes = VariableFeatures(sce.big)) sce.big <- RunTSNE(object = sce.big, dims.use = 1:10) DimPlot(object = sce.big, reduction = "tsne") # 当然也有ICA的选择 # sce.big <- RunICA(sce.big ) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-130236.png) **然后鉴定亚群，看看它们的分群结果** ```r ElbowPlot(sce.big) # 官方建议，下游分析时可以多用几个PCs试试 sce.big <- FindNeighbors(sce.big, dims = 1:20) # 保持和原文一样的15个亚群 sce.big <- FindClusters(sce.big, resolution = 0.23) head(Idents(sce.big), 5) # 新的亚群结果 DimPlot(object = sce.big, reduction = "tsne", group.by = 'RNA_snn_res.0.23') # 原样本分群结果 DimPlot(object = sce.big, reduction = "tsne", group.by = 'orig.ident') ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-131058.png) ```r table(sce.big$orig.ident,sce.big@meta.data$RNA_snn_res.0.23) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-131654.png) 看到，NP有三个主群、G有三个主群、L有三个主群、PI有两个主群 **在之前**[**单细胞天地的推送**](https://mp.weixin.qq.com/s?__biz=MzI1Njk4ODE0MQ==\&mid=2247485216\&idx=1\&sn=ad88d057acfc5e0fefd5ab69ffb46ee8\&scene=21#wechat_redirect)**中，得到了：** ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-143621.png) 一点小差别就是：第一张图的L中0组在第二张图中拆分成了0、1组 > **为了保持参数一致，设置tsne的resolution=0.3再进行测试** 但下面第三张图中L组得到的也是一大群，而且PI得到了3个主群 ![image-20191008224728646](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-144729.png) **对比原文的数据** 它得到了3个NP、3个G、2个L、3个PI，其余的分给了Basal 4群 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2019-10-08-132826.png) --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://jieandze1314.osca.top/05/5.5-duo-ge-shu-ju-ji-de-zheng-he/5.5.1-shi-yong-seurat-de-merge-gong-neng-jin-hang-zheng-he.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.