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

刘小泽写于2020.7.20

1 前言

第一代CEL-Seq由以色列理工学院研究团队于2012发表在Cell Reports

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

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

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

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

数据准备

这里使用的数据是： Grun et al. (2016) 中不同人类供体的胰腺细胞

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
)