Ctrlk

3.9 多样本间差异分析

刘小泽写于2020.7.10+7.15

1 前言

scRNA-seq的一个强大之处在于：方便多个实验条件下的样本比较，比如可以检测药物处理后 (Richard et al. 2018) 或基因修饰后 (Scialdone et al. 2016) 的细胞类型变化，比起常规的单一实验因素能提供更多信息。

对于多样本的scRNA差异分析，主要分成两大类：差异表达（differential expression）和差异丰度（differential abundance），前者检测相同类型细胞在不同条件下表达的变化，后者检测不同条件下细胞类型(或状态等)组成的变化。下面就利用小鼠早期胚胎细胞数据(Pijuan-Sala et al. 2019) 来展示这两种分析。

数据介绍

数据来自小鼠E8.5时期的嵌合体胚胎，每个嵌合体胚胎都是将td-Tomato-positive的胚胎干细胞（ESCs）注射到野生型囊胚中得到的。这里注射的细胞和囊胚细胞除了td-Tomato报告基因，没有其他基因上的差异。利用野生型嵌合体进行研究的目的是为了确定注射过程本身是否引入了其他的差异。

在小鼠胚胎干细胞（ES 细胞）中进行DNA 同源重组，将ES 细胞重新注射到囊胚腔中，可以形成嵌合胚胎，并能在假孕小鼠体内发育，最终发育为嵌合体小鼠，是判定胚胎干细胞系是否具有种系分化能力的唯一方法

sce.chimera数据总共包含6个样本，包含3个重复批次，每个批次2个处理，一个处理样本属于 td-Tomato阳性细胞，另一个处理样本属于阴性。样本是从6-7个嵌合胚胎中利用荧光活化分选得到的。每个样本利用10X Genomics方法建库(Zheng et al. 2017)，得到2000-7000个细胞。

数据准备

scater+scran标准处理模式如下，这个需要下载很多数据，整个过程还是很慢的 数据已分享，可直接下载

其中merging这一步中指定出了批次信息，方便更好地进行数据整合

#--- loading ---#
library(MouseGastrulationData)
sce.chimera <- WTChimeraData(samples=5:10)
sce.chimera

#--- feature-annotation ---#
library(scater)
rownames(sce.chimera) <- uniquifyFeatureNames(
    rowData(sce.chimera)$ENSEMBL, rowData(sce.chimera)$SYMBOL)

#--- quality-control ---#
drop <- sce.chimera$celltype.mapped %in% c("stripped", "Doublet")
sce.chimera <- sce.chimera[,!drop]

#--- normalization ---#
sce.chimera <- logNormCounts(sce.chimera)

#--- variance-modelling ---#
library(scran)
dec.chimera <- modelGeneVar(sce.chimera, block=sce.chimera$sample)
chosen.hvgs <- dec.chimera$bio > 0

#--- merging ---#
library(batchelor)
set.seed(01001001)
merged <- correctExperiments(sce.chimera, 
    batch=sce.chimera$sample, 
    subset.row=chosen.hvgs,
    PARAM=FastMnnParam(
        merge.order=list(
            list(1,3,5), # WT (3 replicates)
            list(2,4,6)  # td-Tomato (3 replicates)
        )
    )
)

#--- clustering ---#
g <- buildSNNGraph(merged, use.dimred="corrected")
clusters <- igraph::cluster_louvain(g)
colLabels(merged) <- factor(clusters$membership)

#--- dimensionality-reduction ---#
merged <- runTSNE(merged, dimred="corrected", external_neighbors=TRUE)
merged <- runUMAP(merged, dimred="corrected", external_neighbors=TRUE)

拿到数据，先要有个初步的认识

其实现在有了分群结果，下一步应该进行marker基因检测以及细胞注释了，但这里的目的不是这些，所以这些步骤略去了，直接使用(Pijuan-Sala et al. 2019) 小鼠早期胚胎发育数据集做的的细胞类型。这个信息放在了merged$celltype.mapped中。不过为了看一下使用他人注释数据用在我们自己的数据质量如何，可以用热图来看二者重叠程度

通过这个图可以看到，我们这里存在多对多的关系。也反映出：针对发育分化的细胞，细胞类型注释是个难题。

2 不同处理间差异基因表达分析 DE analysis

这也是最最常见的分析思路，一般拿到转录组表达量就会先想着把差异分析做一做 针对单细胞的差异分析方法有多种，一般的综述都会提及。这里使用的方法是“模拟bulk转录组表达量”的模式(Tung et al. 2017)

2.1 首先构造一个“拟bulk转录组”的样本

它的思路是：根据不同细胞类型和不同分群，把对应的细胞表达量加起来。比如根据上面得到的细胞分群信息和注释信息

然后根据这两个信息作为分组，将各自的组内细胞表达量分别加起来。

可见大大降低了样本数量，原来是一个细胞当做一个样本对待；现在是整合后的一组细胞当做一个样本。这样做的原因有以下几个：

  • bulk转录组软件的标准差异分析流程需要比较大的表达量，而单细胞中存在很多0表达量

  • 将多个细胞汇成一组，也不是随便挑的几个细胞，就像上面它是根据两个条件选的同组细胞。就像是生物体内的复制过程，虽然这些细胞之间也存在着一些差异，但相比于组外的细胞，它们之间更相似。而如果直接把每个细胞的表达量传给差异分析工具，它会认为每个细胞就是一个独立的生物样本，这也是和真实情况不符的

2.2 进行差异分析

差异分析基于edgeR的 quasi-likelihood (QL)方法，使用负二项广义线性模型(NB GLM) 来处理过度分散的表达量数据。

先来回忆一下bulk转录组的edegR分析怎么做吧：

由于我们前面使用了细胞类型+批次信息两个分组条件，因此这里需要先锁定一个再分析另一个。意思就是我们可以先指定一个细胞类型，然后去进行这3个批次（各2个重复）共6个样本的差异分析，比如可以指定细胞类型是：间质细胞

第一步：创建edgeR对象

第二步：预处理

首先是去除低质量样本（担心影响后面的归一化）：不过这里的文库大小不是像原来一样计算每个样本的表达量加和。因为这里已经将细胞分组，看到的6不是6个细胞，而是6组细胞。因此这里的计算标准是：每组内含有不少于多少的细胞。 例如(Crowell et al. 2019)就定义了如果一组包含少于20个细胞，就是非常低的文库

然后去除低质量基因（目的是增加后面模型构建的准确度，减少多重矫正的压力）：先得到一个logCPM阈值，然后过滤。可以看到下面过滤掉了不少基因

最后，归一化矫正：使用TMM（trimmed mean of M-values）方法，也算得上是calcNormFactors()函数的标志了

第三步：模型并分析

这部分的重点就是这个分组矩阵了，下面锁定样本中批次（pool）的差异，而保留注射导致的差异（tomato）

接下来就是新增部分了，使用estimateDisp()来估计负二项分布negative binomial (NB)

接着用glmQLFit() 估计 quasi-likelihood分布【具体的统计知识可以日后再理解，这里不是重点】

用glmQLFTest()看一下表达量差异，差异基因是logFC不为0（logFC=1也就是treat/control=1）并且FDR小于5%的

看到很少有基因是差异表达的，说明注射这个条件对间质细胞的表达量影响不大

最后，得到一个整合的结果

2.3 循环操作

上面2.2部分是针对间质细胞这一种细胞类型，发现基本没有和注射这个因素相关的差异基因。但是还存在其他很多种细胞类型，因此想到了：如何进行批处理？

最简单的办法是利用scran包的pseudoBulkDGE() ，它会对每个细胞类型进行操作，但依然是需要准备好数据

首先是过滤细胞

直接基于“拟bulk转录组”的数据

然后构建分组矩阵

在2.2中，因为指定了某个细胞类型，所以最后就是得到了这个批次下的6个处理

这里也是，即使使用了全部的细胞类型，但还是要限制为6个处理，因此可以只获得每个处理出现第一次的细胞

最后就可以使用pseudoBulkDGE()

这个结果是这样的：

至于这5列是什么内容，可以看看，里面既有logFC判断上下调，又有FDR判断显著性：

2.4 看一下循环操作的结果

检查一下结果，例如看看每个细胞类型选FDR小于5%（即显著）的差异基因数量。注意到这里有一列NA，它指的是：要么这个细胞类型中基因由于表达量太低被过滤掉，要么就是在这个细胞类型中没办法比较。另外1表示显著上调，-1表示显著下调，0表示差异不显著

还可以看看哪些上、下调差异基因在各种类型细胞中比较多

还可以看：哪些差异基因是在某个细胞类型中特有的

也就是这些基因仅仅在某一个类型的细胞中是差异表达的，在其他类型细胞中不是差异表达。

还可以将这些基因排个序（例如下面按Pvalue排序）

如果说某些细胞类型没有重复或没有对照，它们就不好做差异分析，这些“失败”的类型可以提取出来：

3 差异细胞丰度分析 DA analysis

检测不同条件下细胞类型(或状态等)组成的变化 如果说DE analysis重点在基因的表达量变化，那么DA analysis就着眼于细胞数量的变化

差异丰度分析在流式细胞术中应用较多，常被用于检查不同条件对复杂细胞群体组成的影响。我们这里也类比FACS技术，对scRNA的细胞也做一个“标记”，然后看看处理前后标记的丰度变化程度。

还是先将细胞按两个维度（细胞类型+处理）归类

这个abundances结果，依然可以类比为原来的差异表达分析，即行为基因列为样本，我们要对样本间进行差异比较（这里也就是对6个处理进行比较）

3.1 基于edgeR做DA analysis

首先为列的样本信息增加一些内容：

然后也还是按行过滤，只不过这里不是根据基因的表达量了，而是按照细胞数量

关于filterByExpr函数是如何过滤的【参考帮助文档】： 宗旨：keeps genes that have at least min.count reads in a worthwhile number samples 依据：the filtering keeps genes that have count-per-million (CPM) above k in n samples 定义：k is determined by min.countand by the sample library sizes； n is determined by the design matrix.

这里不再利用calcNormFactors()这么复杂的归一化方法，我们只需要简单矫正一个文库大小即可，具体原因下面再讨论

下面的步骤其实和差异基因分析很像：

注意：

这几个步骤里面的重点其实一直都是model.matrix的设置，其他的统计计算我们并不关心，根据edgeR的官方文档，设置的前后顺序很重要：

也就是说，这里的6个处理样本同时具有两个属性，即pool和tomato，因此要对这6个样本进行差异分析，那就必须把这两个属性都考虑进去，但哪个更重要我们需要掂量一下

pool 即每个处理的来源是不关心的，我们关心的是是否有tomato处理。因此pool放在了前面

3.2 处理细胞组成的影响

前面说不再利用calcNormFactors()这么复杂的归一化方法，因为这种方法的假设是：基本每一行（feature）在各个列都是存在的，也就是说大部分基因在所有细胞中都是存在的。但是对于我们这里的DA分析来讲，行（细胞类型）不一定在所有的列（不同处理）中都有。很多实验中只包含少数几个细胞类型，如果还是按照这种计算方法，会对真实的细胞数量产生较大的干扰。

因此，默认的操作就是基于不同处理样本的细胞总数进行归一化处理。这个在技术和统计上说得过去，但落实到生物意义上又有争议。如果一个处理中某种类型的细胞丰度很大，在归一化后，就会导致其他几个类型的细胞占比下降，从而可能得到错误结论：认为这个处理下的其他几种类型的细胞数量会减少。

下面有几种处理方式：

方法一：首先假设大部分的细胞类型在各个处理中都是存在的

这个假设对于野生型/对照组来说是合理的，因为它不需要考虑注射处理后产生的不同结果。 这样就可以利用calcNormFactors()去计算

方法二：移除那些细胞数量很多的细胞类型

采用的常规的归一化处理方法，只是考虑到那些丰度大的细胞类型对整体的影响，于是可以去除它们。

比如我们认为ExE ectoderm这个细胞类型中细胞数量太多

上一页3.8 批次效应处理下一页3.10 检测Doublet

最后更新于