单细胞交响乐
  • 前言:我与《单细胞交响乐》的缘分
  • 1 准备篇:背景知识
    • 1.1 数据结构
    • 1.2 总览 | 从实验到分析
  • 2 积累篇:文献阅读
    • 2.1.1 综述 | 2019-单细胞转录组分析最佳思路
    • 2.1.2 综述 | 2018-单细胞捕获平台
    • 2.1.3 综述 | 2017-scRNA中的细胞聚类分群
    • 2.1.4 综述 | scRNA已经开发出超过1000款工具了,你用过几种?
    • 2.1.5 综述 | 2021-单细胞测序的微流控技术应用
    • 2.2.1 研究 | 2018-单细胞转录组探索癌症免疫治疗获得性抗性机理
    • 2.2.2 研究 | 2018-人类结直肠癌单细胞多组学分析
    • 2.2.3 研究 | 2020-单细胞分析揭示葡萄膜黑色素瘤新的进化复杂性
    • 2.2.4 研究 | 2020-COVID-19病人支气管免疫细胞单细胞测序分析
    • 2.2.5 研究 | 2020-原汁原味读--单细胞肿瘤免疫图谱
    • 2.2.6 研究 | 2021-多发性骨髓瘤发展过程中肿瘤和免疫细胞的共同进化
    • 2.2.7 研究 | 2021-多个组织的成纤维细胞图谱
    • 2.2.8 研究 | 2021-多组学分析肺结核队列的记忆T细胞状态
    • 2.2.9 研究 | 2021-CancerSCEM: 人类癌症单细胞表达图谱数据库
    • 2.2.10 研究| 2021-单细胞转录组分析COVID-19重症患者肺泡巨噬细胞亚型
    • 2.2.11 研究 |2021-单细胞转录组揭示肺腺癌特有的肿瘤微环境
    • 2.2.12 研究 | 2021-单细胞转录组揭示乳头状甲状腺癌起始与发展
    • 2.2.13 研究 | 2021-解析食管鳞癌化疗病人的单细胞转录组
    • 2.2.14 研究 | 2021-单细胞水平看骨髓瘤的细胞状态和基因调控
    • 2.3.1 算法|2020-BatchBench比较scRNA批次矫正方法
    • 2.3.2 算法 | 2021-scPhere——用地球仪来展示降维结果
    • 2.3.3 算法 | 2021-单细胞差异分析方法评测
    • 2.3.4 算法 | 2021-细胞分群新方法——CNA(co-varying neighborhood analysis)
    • 2.3.5 工具 | 2018-iSEE:单细胞数据可视化辅助网页工具
    • 2.3.6 工具 | 2021-MACA: 一款自动注释细胞类型的工具
    • 2.3.7 工具 | 2021-一个很有想法的工具——Ikarus,想要在单细胞水平直接鉴定肿瘤细胞
  • 3 流程篇:分析框架
    • 3.1 质控
    • 3.2 归一化
    • 3.3 挑选表达量高变化基因
    • 3.4 降维
    • 3.5 聚类
    • 3.6 Marker/标记基因检测
    • 3.7 细胞类型注释
    • 3.8 批次效应处理
    • 3.9 多样本间差异分析
    • 3.10 检测Doublet
    • 3.11 细胞周期推断
    • 3.12 细胞轨迹推断
    • 3.13 与蛋白丰度信息结合
    • 3.14 处理大型数据
    • 3.15 不同R包数据的相互转换
  • 4 实战篇:活学活用
    • 4.1 实战一 | Smart-seq2 | 小鼠骨髓
    • 4.2 实战二 | STRT-Seq | 小鼠大脑
    • 4.3 实战三 | 10X | 未过滤的PBMC
    • 4.4 实战四 | 10X | 过滤后的PBMC
    • 4.5 实战五 | CEL-seq2 | 人胰腺细胞
    • 4.6 实战六 | CEL-seq | 人胰腺细胞
    • 4.7 实战七 | SMARTer | 人胰腺细胞
    • 4.8 实战八 | Smart-seq2 | 人胰腺细胞
    • 4.9 实战九 | 不同技术数据整合 | 人胰腺细胞
    • 4.10 实战十 | CEL-seq | 小鼠造血干细胞
    • 4.11 实战十一 | Smart-seq2 | 小鼠造血干细胞
    • 4.12 实战十二 | 10X | 小鼠嵌合体胚胎
    • 4.13 实战十三 | 10X | 小鼠乳腺上皮细胞
    • 4.14 | 实战十四 | 10X | HCA计划的38万骨髓细胞
  • 5 补充篇:开拓思路
    • 5.1 10X Genomics概述
      • 5.1.1 10X Genomics 问题集锦
    • 5.2 CellRanger篇
      • 5.2.1 CellRanger实战(一)数据下载
      • 5.2.2 CellRanger实战(二) 使用前注意事项
      • 5.2.3 CellRanger实战(三) 使用初探
      • 5.2.4 CellRanger实战(四)流程概览
      • 5.2.5 CellRanger实战(五) 理解count输出的结果
    • 5.3 Seurat的使用
      • 5.3.1 Seurat V3 | 实战之2700 PBMCs分析
      • 5.3.2 Seurat V3 | 如何改造Seurat包的DoHeatmap函数?
      • 5.3.3 scRNA的3大R包对比
      • 5.3.4 Seurat两种数据比较:integrated vs RNA assay
      • 5.3.5 seurat 的几种findmaker比较
    • 5.4 Monocle的使用
      • 5.4.1 Monocle V3实战
    • 5.5 多个数据集的整合
      • 5.5.1 使用Seurat的merge功能进行整合
      • 5.5.2 如何使用sctransform去除批次效应
由 GitBook 提供支持
在本页
  • 前言
  • 开始操作
  • 第一步:准备原始测序数据
  • 第二步:整理数据
  • 第三步:批量读取成10X对象
  • 第四步:组合
  • 补充

这有帮助吗?

  1. 5 补充篇:开拓思路
  2. 5.5 多个数据集的整合

5.5.1 使用Seurat的merge功能进行整合

刘小泽写于19.10.8

上一页5.5 多个数据集的整合下一页5.5.2 如何使用sctransform去除批次效应

最后更新于4年前

这有帮助吗?

前言

单细胞数据未来会朝着多样本发展,因此数据整合是一项必备技能。cellranger中自带了aggr的整合功能,而这篇文章(Differentiation dynamics of mammary epithelial cells revealed by single-cell RNA-sequencing)的作者也正是这么做得到的组合后的表达矩阵,然后用Read10X读入

关于文章

这是发表在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这些对象

开始操作

第一步:准备原始测序数据

感觉手机下载速度就是比电脑快,这个文件在手机上下载3分钟,电脑预计时间1小时

下载后解压,整个过程直接在Rstudio中的Terminal直接完成

第二步:整理数据

思路:根据中间的分组信息(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脚本

# 列出当前目录下所有开头是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:

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】

  • 说找不到基因文件,那么就要看看测序数据是不是解压后的

Error in Read10X() : Gene name or features file missing

第三步:批量读取成10X对象

Read10X() + CreateSeuratObject()

# 因为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)

第四步:组合

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两个seurat对象(原始数据)

需要注意的是,组合数据时需要注明每个数据的名称,使用add.cell.ids参数指定

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 设成一个向量,就可以指定其他的数据

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指定

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对比:

#################
# 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 . . .

上面的组合多个数据就结束了,接下来是检查组合后的分群结果

首先检查原样本分群结果

# 归一化+标准化(移除了不想要的差异来源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 )

然后鉴定亚群,看看它们的分群结果

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')
table(sce.big$orig.ident,sce.big@meta.data$RNA_snn_res.0.23)

看到,NP有三个主群、G有三个主群、L有三个主群、PI有两个主群

一点小差别就是:第一张图的L中0组在第二张图中拆分成了0、1组

对比原文的数据

它得到了3个NP、3个G、2个L、3个PI,其余的分给了Basal 4群

我们下载第一个:GSE106273_RAW.tar(183.5 Mb)

官网的merge教程在:

在之前中,得到了:

为了保持参数一致,设置tsne的resolution=0.3再进行测试 但下面第三张图中L组得到的也是一大群,而且PI得到了3个主群

https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE106273&format=file
https://satijalab.org/seurat/v3.1/merge_vignette.html
单细胞天地的推送
image-20191008224728646