单细胞交响乐
  • 前言:我与《单细胞交响乐》的缘分
  • 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 提供支持
在本页
  • 1 前言
  • 数据准备
  • ID转换
  • 2 质控
  • 数据备份
  • 然后进行过滤,并让函数注意批次信息
  • 作图
  • 【重点】发现了问题
  • 最后把过滤条件应用在原数据
  • 3 归一化
  • 4 找表达量高变化基因
  • 5 矫正批次效应
  • 6 降维 + 聚类
  • 降维
  • 聚类

这有帮助吗?

  1. 4 实战篇:活学活用

4.5 实战五 | CEL-seq2 | 人胰腺细胞

刘小泽写于2020.7.20

上一页4.4 实战四 | 10X | 过滤后的PBMC下一页4.6 实战六 | CEL-seq | 人胰腺细胞

最后更新于4年前

这有帮助吗?

1 前言

第一代CEL-Seq由以色列理工学院研究团队于

CEL-Seq的全称是:Cell expression by linear amplification and sequencing,采用线性扩增的测序方法,其主要优势在于错误率比较低,但是和PCR一样都存在序列偏向性。另外还具有以下优点:

  • 使用barcode允许多种类型的细胞同时分析;

  • 每个细胞在一个管中进行处理,避免了样本之间的污染

后来在2016年,CEL-seq2诞生,与早期版本相比具有低成本、高灵敏度,并缩短了实验操作实验时间

数据准备

这里使用的数据是: 中不同人类供体的胰腺细胞

library(scRNAseq)
sce.grun <- GrunPancreasData()

sce.grun
# class: SingleCellExperiment 
# dim: 20064 1728 
# metadata(0):
#   assays(1): counts
# rownames(20064): A1BG-AS1__chr19 A1BG__chr19 ...
# ZZEF1__chr17 ZZZ3__chr1
# rowData names(2): symbol chr
# colnames(1728): D2ex_1 D2ex_2 ... D17TGFB_95
# D17TGFB_96
# colData names(2): donor sample
# reducedDimNames(0):
#   altExpNames(1): ERCC

rowData(sce.grun)
# DataFrame with 20064 rows and 2 columns
# symbol         chr
# <character> <character>
#   A1BG-AS1__chr19    A1BG-AS1       chr19
# A1BG__chr19            A1BG       chr19
# A1CF__chr10            A1CF       chr10
# A2M-AS1__chr12      A2M-AS1       chr12
# A2ML1__chr12          A2ML1       chr12
# ...                     ...         ...
# ZYG11A__chr1         ZYG11A        chr1
# ZYG11B__chr1         ZYG11B        chr1
# ZYX__chr7               ZYX        chr7
# ZZEF1__chr17          ZZEF1       chr17
# ZZZ3__chr1             ZZZ3        chr1

ID转换

先将symbol转为Ensembl ID

library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=rowData(sce.grun)$symbol,
    keytype="SYMBOL", column="ENSEMBL")

接下来有两种处理方式:

第一种:将Ensembl ID加到原来数据中

library(scater)
rowData(sce.grun)$ensembl <- uniquifyFeatureNames(
  rowData(sce.grun)$symbol, gene.ids)
rowData(sce.grun)
# DataFrame with 20064 rows and 3 columns
# symbol         chr         ensembl
# <character> <character>     <character>
#   A1BG-AS1__chr19    A1BG-AS1       chr19 ENSG00000268895
# A1BG__chr19            A1BG       chr19 ENSG00000121410
# A1CF__chr10            A1CF       chr10 ENSG00000148584
# A2M-AS1__chr12      A2M-AS1       chr12 ENSG00000245105
# A2ML1__chr12          A2ML1       chr12 ENSG00000166535
# ...                     ...         ...             ...
# ZYG11A__chr1         ZYG11A        chr1 ENSG00000203995
# ZYG11B__chr1         ZYG11B        chr1 ENSG00000162378
# ZYX__chr7               ZYX        chr7 ENSG00000159840
# ZZEF1__chr17          ZZEF1       chr17 ENSG00000074755
# ZZZ3__chr1             ZZZ3        chr1 ENSG00000036549

第二种:将没有匹配的NA去掉,并且去掉重复的行(因为可能存在一个symbol对应多个Ensembl ID的情况)【这里就试一下这种方法】

keep <- !is.na(gene.ids) & !duplicated(gene.ids)
> sum(is.na(gene.ids));sum(duplicated(gene.ids))
[1] 2487
[1] 2515

sce.grun <- sce.grun[keep,]
rownames(sce.grun) <- gene.ids[keep]
> dim(sce.grun) # 过滤掉2000多基因
[1] 17548  1728

2 质控

前几次实战的质控都很顺利,但并不意味着实际操作中这一步就可以跑跑代码解决 这一次,就遇到了一个常规代码解决不了的问题,需要思考+调整

数据备份

把unfiltered数据主要用在质控的探索上

unfiltered <- sce.grun

看到这里的数据中没有线粒体基因

table(rowData(sce.grun)$chr)
# 
# chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 
# 1766   672  1071   947   328   551   532   740  1050   256 
# chr19  chr2 chr20 chr21 chr22  chr3  chr4  chr5  chr6  chr7 
# 1245  1142   484   188   400  1014   693   791   907   825 
# chr8  chr9  chrX  chrY 
# 604   664   658    20

倒是有几个donor的批次信息

table(colData(sce.grun)$donor)
# 
# D10 D17  D2  D3  D7 
# 288 480  96 480 384

我们之前进行质控的时候,一般采用的先使用perCellQCMetrics 指定线粒体计算,然后用 quickPerCellQC 指定ERCC+线粒体进行过滤的策略。现在既然没有线粒体的信息,那么在第一步的perCellQCMetrics就不需要加subset参数了

stats <- perCellQCMetrics(sce.grun)
# 参数subsets的含义: used to identify interesting feature subsets such as ERCC spike-in transcripts or mitochondrial genes.
# > colnames(stats)
# [1] "sum"                   "detected"             
# [3] "percent_top_50"        "percent_top_100"      
# [5] "percent_top_200"       "percent_top_500"      
# [7] "altexps_ERCC_sum"      "altexps_ERCC_detected"
# [9] "altexps_ERCC_percent"  "total"

然后进行过滤,并让函数注意批次信息

qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
                     batch=sce.grun$donor)
colSums(as.matrix(qc), na.rm=TRUE)
# low_lib_size            low_n_features high_altexps_ERCC_percent 
# 85                        93                        88 
# discard 
# 129

作图

colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc$discard

gridExtra::grid.arrange(
  plotColData(unfiltered, x="donor", y="sum", colour_by="discard") +
    scale_y_log10() + ggtitle("Total count"),
  plotColData(unfiltered, x="donor", y="detected", colour_by="discard") +
    scale_y_log10() + ggtitle("Detected features"),
  plotColData(unfiltered, x="donor", y="altexps_ERCC_percent",
              colour_by="discard") + ggtitle("ERCC percent"),
  ncol=3
)

【重点】发现了问题

质控并非是一帆风顺的,遇到问题应该知道如何探索解决

看到ERCC那个图,很明显D10和D3没有被正确过滤,它们中高ERCC的细胞还是没有被去掉

质控过滤有一个前提条件:大部分的细胞都是高质量的,但显然这里的D10、D3不符合这个要求,因此我们需要再次只根据D17、D2、D7这三个“好一点的样本”进行过滤

# 重新计算过滤条件
qc2 <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
    batch=sce.grun$donor,
    subset=sce.grun$donor %in% c("D17", "D7", "D2"))

colSums(as.matrix(qc2), na.rm=TRUE)
# low_lib_size            low_n_features high_altexps_ERCC_percent                   discard 
# 450                       511                       606                       665

这次再作图看看

colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- qc2$discard

gridExtra::grid.arrange(
  plotColData(unfiltered, x="donor", y="sum", colour_by="discard") +
    scale_y_log10() + ggtitle("Total count"),
  plotColData(unfiltered, x="donor", y="detected", colour_by="discard") +
    scale_y_log10() + ggtitle("Detected features"),
  plotColData(unfiltered, x="donor", y="altexps_ERCC_percent",
              colour_by="discard") + ggtitle("ERCC percent"),
  ncol=3
)

看到这次对D10、D3的过滤力度大了很多,基本满意

最后把过滤条件应用在原数据

sce.grun <- sce.grun[,!qc$discard]

3 归一化

也是使用预分群+去卷积计算size factor的方法

library(scran)
set.seed(1000)
clusters <- quickCluster(sce.grun)
sce.grun <- computeSumFactors(sce.grun, cluster=clusters)
sce.grun <- logNormCounts(sce.grun)

summary(sizeFactors(sce.grun))
# Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 0.09581  0.50459  0.79505  1.00000  1.22212 12.16491

看看两种归一化方法的差异

plot(librarySizeFactors(sce.grun), sizeFactors(sce.grun), pch=16,
     col=as.integer(factor(sce.grun$donor)),
     xlab="Library size factors", ylab="Deconvolution factors", log="xy")
abline(a=0, b=1, col="red")

也是再一次证明了去卷积方法比常规的方法更加精细,对批次层面也会纳入考量

4 找表达量高变化基因

既然有批次的信息,那就加到构建模型中去,而且有ERCC,就选用modelGeneVarWithSpikes这个方法

block <- paste0(sce.grun$sample, "_", sce.grun$donor)
dec.grun <- modelGeneVarWithSpikes(sce.grun, spikes="ERCC", block=block)
top.grun <- getTopHVGs(dec.grun, prop=0.1)

5 矫正批次效应

使用FastMNN方法:Correct for batch effects in single-cell expression data using a fast version of the mutual nearest neighbors (MNN) method.

library(batchelor)
set.seed(1001010)
merged.grun <- fastMNN(sce.grun, subset.row=top.grun, batch=sce.grun$donor)
merged.grun
# class: SingleCellExperiment 
# dim: 1467 1063 
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(1467): ENSG00000172023 ENSG00000172016 ...
# ENSG00000144867 ENSG00000179820
# rowData names(1): rotation
# colnames(1063): D2ex_1 D2ex_2 ... D17TGFB_94
# D17TGFB_95
# colData names(1): batch
# reducedDimNames(2): corrected
# altExpNames(0):

检查一下结果,使用lost.var ,值越大表示丢失的真实生物异质性越多

metadata(merged.grun)$merge.info$lost.var
# D10         D17          D2         D3         D7
# [1,] 0.030843209 0.030956515 0.000000000 0.00000000 0.00000000
# [2,] 0.007129994 0.011275730 0.036836491 0.00000000 0.00000000
# [3,] 0.003528589 0.005027722 0.006846244 0.05020422 0.00000000
# [4,] 0.012070814 0.014673302 0.013972521 0.01267586 0.05281354
  • It contains a matrix of the variance lost in each batch (column) at each merge step (row).

  • Large proportions of lost variance (>10%) suggest that correction is removing genuine biological heterogeneity.

6 降维 + 聚类

降维

set.seed(100111)
merged.grun <- runTSNE(merged.grun, dimred="corrected")

聚类

snn.gr <- buildSNNGraph(merged.grun, use.dimred="corrected")
colLabels(merged.grun) <- factor(igraph::cluster_walktrap(snn.gr)$membership)

# 聚类的分群与批次之间的对比
table(Cluster=colLabels(merged.grun), Donor=merged.grun$batch)
##        Donor
## Cluster D10 D17  D2  D3  D7
##      1   32  71  33  80  29
##      2    5   7   3   4   4
##      3    3  44   0   0  13
##      4   11 119   0   0  55
##      5   12  73  30   3  78
##      6   14  29   3   2  64
##      7    1   9   0   0   7
##      8    2   5   2   3   4
##      9    5  13   0   0  10
##      10  15  33  11  10  37
##      11   5  18   0   1  33
##      12   4  13   0   0   1

最后作图看看批次效应矫正如何

gridExtra::grid.arrange(
    plotTSNE(merged.grun, colour_by="label"),
    plotTSNE(merged.grun, colour_by="batch"),
    ncol=2
)
2012发表在Cell Reports
Grun et al. (2016)