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)

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

> merged
class: SingleCellExperiment 
dim: 14699 19426 
metadata(2): merge.info pca.info
assays(3): reconstructed counts logcounts
rownames(14699): Xkr4 Rp1 ... Vmn2r122
  CAAA01147332.1
rowData names(3): rotation ENSEMBL SYMBOL
colnames(19426): cell_9769 cell_9770 ...
  cell_30701 cell_30702
colData names(13): batch cell ... sizeFactor
  label
reducedDimNames(3): corrected TSNE UMAP
altExpNames(0):

# 有哪些样本信息
> names(colData(merged))
 [1] "batch"           "cell"           
 [3] "barcode"         "sample"         
 [5] "stage"           "tomato"         
 [7] "pool"            "stage.mapped"   
 [9] "celltype.mapped" "closest.cell"   
[11] "doub.density"    "sizeFactor"     
[13] "label"  

# 看下三个批次
> table(merged$pool)

    3     4     5 
 3324  5644 10458 

# 看下两种处理
> table(merged$tomato)

FALSE  TRUE 
10331  9095 

# 看下分群结果
> table(colLabels(merged))

   1    2    3    4    5    6    7    8    9   10 
 947  112  868  680  606  507 2208  424 1259  252 
  11   12   13   14   15   16   17   18   19   20 
 104  727   58  423 1044  872  432 1264  454 1022 
  21   22   23   24   25   26 
 211  160  156 1640  860 2136 

# 当然也可以把三个批次和分群结果放一起看
table(colLabels(merged), merged$pool)
# 或者把两个处理和分群结果放一起看
table(colLabels(merged), merged$tomato)

# 再直观一点，看看数据整合的效果如何，也就是批次效应处理如何
gridExtra::grid.arrange(
    plotTSNE(merged, colour_by="tomato", text_by="label"),
    plotTSNE(merged, colour_by=data.frame(pool=factor(merged$pool))),
    ncol=2
)

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

by.label <- table(colLabels(merged), merged$celltype.mapped)
pheatmap::pheatmap(log2(by.label+1), cluster_cols=FALSE, cluster_rows=FALSE,
    color=viridis::viridis(101))

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

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

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

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

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

> colData(merged)[,c("celltype.mapped", "sample")]
DataFrame with 19426 rows and 2 columns
                        celltype.mapped    sample
                            <character> <integer>
cell_9769                    Mesenchyme         5
cell_9770                   Endothelium         5
cell_9771                     Allantois         5
cell_9772                    Erythroid3         5
cell_9773                    Erythroid1         5
...                                 ...       ...
cell_30698                   Erythroid2        10
cell_30699                   Erythroid3        10
cell_30700             Surface ectoderm        10
cell_30701 Forebrain/Midbrain/Hindbrain        10
cell_30702                   Erythroid3        10

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

summed <- aggregateAcrossCells(merged, 
    id=colData(merged)[,c("celltype.mapped", "sample")])

> dim(merged);dim(summed)
[1] 14699 19426
[1] 14699   186

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

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

2.2 进行差异分析

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

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

rm(list = ls())
options(stringsAsFactors = F)

load('airway.expreSet.Rdata')
library(edgeR)
##  第一步：创建edgeR对象
# 这里需要表达矩阵和分组信息
e <- DGEList(counts=expr,group=factor(grp))

## 第二步：预处理
# 进行一下过滤
keep <- rowSums(cpm(e)>1) >= 2
e <- e[keep, , keep.lib.sizes=FALSE]
# 进行一下校正
e$samples$lib.size <- colSums(e$counts)
e <- calcNormFactors(e)
e$samples

DEG=e
## 第三步：创建模型并分析
# 这里得到分组矩阵（实验设计矩阵）
design <- model.matrix(~0+factor(grp))
rownames(design)<-colnames(DEG)
colnames(design)<-levels(factor(grp))

DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)

fit <- glmFit(DEG, design)
# edgeR User guide (Page. 30 => "GLM Approach")
lrt <- glmLRT(fit,  contrast=c(-1,1)) # accoding to design to modify
# or we can use another way (glmQLFTest):
# CasevsCtrl <- makeContrasts(Case-Ctrl=trt-untrt, levels=design)
# lrt <- glmQLFTest(fit,contrast=CasevsCtrl)

nrDEG=topTags(lrt, n=nrow(DEG))
nrDEG=as.data.frame(nrDEG)

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

label <- "Mesenchyme"
current <- summed[,label==summed$celltype.mapped]
# 于是从原来的186组样本又再次过滤为了6组样本
> dim(current)
[1] 14699     6

第一步：创建edgeR对象

library(edgeR)
y <- DGEList(counts(current), samples=colData(current))
y

第二步：预处理

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

discarded <- current$ncells < 20
y <- y[,!discarded]
summary(discarded)
##    Mode   FALSE 
## logical       6

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

keep <- filterByExpr(y, group=current$tomato)
y <- y[keep,]
summary(keep)
##    Mode   FALSE    TRUE 
## logical    9011    5688

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

y <- calcNormFactors(y)

第三步：模型并分析

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

design <- model.matrix(~factor(pool) + factor(tomato), y$samples)
design
##         (Intercept) factor(pool)4 factor(pool)5 factor(tomato)TRUE
## Sample1           1             0             0                  1
## Sample2           1             0             0                  0
## Sample3           1             1             0                  1
## Sample4           1             1             0                  0
## Sample5           1             0             1                  1
## Sample6           1             0             1                  0
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$`factor(pool)`
## [1] "contr.treatment"
## 
## attr(,"contrasts")$`factor(tomato)`
## [1] "contr.treatment"

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

y <- estimateDisp(y, design)
summary(y$trended.dispersion)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0103  0.0167  0.0213  0.0202  0.0235  0.0266

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

fit <- glmQLFit(y, design, robust=TRUE)

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

res <- glmQLFTest(fit, coef=ncol(design))
summary(decideTests(res))
##        factor(tomato)TRUE
## Down                    8
## NotSig               5672
## Up                      8

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

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

topTags(res)
## Coefficient:  factor(tomato)TRUE 
##            logFC logCPM       F    PValue       FDR
## Phlda2   -4.3874  9.934 1638.59 1.812e-16 1.031e-12
## Erdr1     2.0691  8.833  356.37 1.061e-11 3.017e-08
## Mid1      1.5191  6.931  120.15 1.844e-08 3.497e-05
## H13      -1.0596  7.540   80.80 2.373e-07 2.527e-04
## Kcnq1ot1  1.3763  7.242   83.31 2.392e-07 2.527e-04
## Akr1e1   -1.7206  5.128   79.31 2.665e-07 2.527e-04
## Zdbf2     1.8008  6.797   83.66 6.809e-07 5.533e-04
## Asb4     -0.9235  7.341   53.45 2.918e-06 2.075e-03
## Impact    0.8516  7.353   50.31 4.145e-06 2.620e-03
## Lum      -0.6031  9.275   41.67 1.205e-05 6.851e-03

2.3 循环操作

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

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

首先是过滤细胞

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

summed.filt <- summed[,summed$ncells >= 20]
> dim(summed.filt)
[1] 14699   126

然后构建分组矩阵

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

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

targets <- colData(merged)[!duplicated(merged$sample),]
> dim(targets)
[1]  6 13

# 再根据这个分组信息构建矩阵
design <-  model.matrix(~factor(pool) + factor(tomato), data=targets)
rownames(design) <- targets$sample

最后就可以使用pseudoBulkDGE()

library(scran)
de.results <- pseudoBulkDGE(summed.filt, 
    sample=summed.filt$sample,
    label=summed.filt$celltype.mapped,
    design=design,
    coef=ncol(design),
    condition=targets$tomato 
)

这个结果是这样的：

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

> de.results$Allantois[1:3,1:5]
DataFrame with 3 rows and 5 columns
          logFC    logCPM         F    PValue       FDR
      <numeric> <numeric> <numeric> <numeric> <numeric>
Xkr4         NA        NA        NA        NA        NA
Rp1          NA        NA        NA        NA        NA
Sox17 -0.139844     4.909  0.335911  0.562224  0.891031

2.4 看一下循环操作的结果

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

is.de <- decideTestsPerLabel(de.results, threshold=0.05)
> dim(is.de)
[1] 14699    26

# 记录了每个基因在每个类型中是否差异
> is.de[1:3,1:3]
      Allantois Blood progenitors 2 Cardiomyocytes
Xkr4         NA                  NA             NA
Rp1          NA                  NA             NA
Sox17         0                  NA             NA

summarizeTestsPerLabel(is.de)
##                                -1    0  1    NA
## Allantois                      69 5048 66  9516
## Blood progenitors 2             1 2472  2 12224
## Cardiomyocytes                  6 4361  5 10327
## Caudal epiblast                 0    0  0 14699
## Caudal Mesoderm                 0    0  0 14699
## Def. endoderm                   0    0  0 14699
## Endothelium                     3 3222  6 11468
## Erythroid1                     12 3035 25 11627
## Erythroid2                      5 3389  8 11297
## Erythroid3                     13 5048 16  9622
## ExE ectoderm                    0    0  0 14699
## ExE mesoderm                    2 5097 10  9590
## Forebrain/Midbrain/Hindbrain    8 6226 11  8454
## Gut                             5 4482  6 10206
## Haematoendothelial progenitors  7 4347 17 10328
## Intermediate mesoderm           6 3256  8 11429
## Mesenchyme                      8 5672  8  9011
## Neural crest                    6 3311  8 11374
## NMP                             6 4107 10 10576
## Paraxial mesoderm               4 4756  5  9934
## Parietal endoderm               0    0  0 14699
## Pharyngeal mesoderm             2 5082  9  9606
## Rostral neurectoderm            0    0  0 14699
## Somitic mesoderm                7 2948 13 11731
## Spinal cord                     7 4591  7 10094
## Surface ectoderm                9 5556  8  9126

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

# 上调
up.de <- is.de > 0 & !is.na(is.de)
head(sort(rowMeans(up.de), decreasing=TRUE), 10)
##     Mid1    Erdr1   Impact    Mcts2     Nnat Kcnq1ot1  Slc38a4    Zdbf2 
##   0.7692   0.6538   0.5385   0.5000   0.5000   0.5000   0.3846   0.3462 
##     Hopx     Peg3 
##   0.3462   0.2308

# 下调
down.de <- is.de < 0 & !is.na(is.de)
head(sort(rowMeans(down.de), decreasing=TRUE), 10)
##        Akr1e1          Xist        Cdkn1c        Phlda2           H13 
##       0.61538       0.57692       0.57692       0.57692       0.46154 
##         Wfdc2         Hbb-y         Grb10 B930036N10Rik         Pink1 
##       0.19231       0.11538       0.11538       0.07692       0.07692

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

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

# 首先复制一份差异分析的结果
remotely.de <- decideTestsPerLabel(de.results, threshold=0.5)
# 从中把非差异基因选出来
not.de <- remotely.de==0 | is.na(remotely.de)
# 如果我们只关注Allantois这个细胞类型，那就把剩余细胞类型设为“其他”
other.labels <- setdiff(colnames(not.de), "Allantois")
# 从差异基因中找Allantois的，同时需要排除那些也存在于其他细胞类型的基因：即选择非差异基因中其他细胞类型均为True的，也就是下面rowMeans(not.de[,other.labels])==1
unique.degs <- is.de[,"Allantois"]!=0 & rowMeans(not.de[,other.labels])==1
# 提取基因名
unique.degs <- names(which(unique.degs))

# 看一下差别：如果只看Allantois中存在的是14699个，如果看它特有的是44个
> length(is.de[,"Allantois"]!=0)
[1] 14699
> length(unique.degs)
[1] 44

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

de.allantois <- de.results$Allantois
de.allantois <- de.allantois[order(de.allantois$PValue),]
de.allantois <- de.allantois[rownames(de.allantois) %in% unique.degs,]

> head(de.allantois)
DataFrame with 6 rows and 5 columns
             logFC    logCPM         F      PValue         FDR
         <numeric> <numeric> <numeric>   <numeric>   <numeric>
Slc22a18 -4.629906   3.87065  118.8863 2.19531e-27 1.42228e-24
Eif2s3y   1.757650   5.84561   74.0124 1.01564e-17 4.04928e-15
Rbp4      1.978988   4.38295   32.7443 1.11016e-08 1.74363e-06
Cfc1     -0.884383   5.66713   23.7376 1.13690e-06 1.33922e-04
Pdgfa    -0.423246   8.70475   23.4935 1.28991e-06 1.48569e-04
H3f3b     0.269312  12.07841   22.8530 1.79717e-06 2.02494e-04

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

metadata(de.results)$failed
## [1] "Caudal epiblast"      "Caudal Mesoderm"      "Def. endoderm"       
## [4] "ExE ectoderm"         "Parietal endoderm"    "Rostral neurectoderm"

3 差异细胞丰度分析 DA analysis

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

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

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

abundances <- table(merged$celltype.mapped, merged$sample) 
> class(abundances)
[1] "table"

abundances <- unclass(abundances) #它的作用显而易见
> class(abundances)
[1] "matrix" "array" 

> dim(abundances)
[1] 34  6
head(abundances)
##                      
##                        5  6   7   8   9  10
##   Allantois           97 15 139 127 318 259
##   Blood progenitors 1  6  3  16   6   8  17
##   Blood progenitors 2 31  8  28  21  43 114
##   Cardiomyocytes      85 21  79  31 174 211
##   Caudal epiblast      2  2   0   0  22  45
##   Caudal Mesoderm     10 10   9   3  10  29

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

3.1 基于edgeR做DA analysis

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

extra.info <- colData(merged)[match(colnames(abundances), merged$sample),]
y.ab <- DGEList(abundances, samples=extra.info)
y.ab

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

keep <- filterByExpr(y.ab, group=y.ab$samples$tomato)
y.ab <- y.ab[keep,]
summary(keep)
##    Mode   FALSE    TRUE 
## logical      10      24

关于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()这么复杂的归一化方法，我们只需要简单矫正一个文库大小即可，具体原因下面再讨论

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

# 分组矩阵
design <- model.matrix(~factor(pool) + factor(tomato), y.ab$samples)
# 对每个细胞类型估计负二项分布（由于细胞数量远不如基因数量，点不够就不用设置趋势线）
y.ab <- estimateDisp(y.ab, design, trend="none")
# 再进行一个QL（ quasi-likelihood） 分布
fit.ab <- glmQLFit(y.ab, design, robust=TRUE, abundance.trend=FALSE)
# 检验
res <- glmQLFTest(fit.ab, coef=ncol(design))
summary(decideTests(res))
##        factor(tomato)TRUE
## Down                    1
## NotSig                 22
## Up                      1

# 看结果
topTags(res)
## Coefficient:  factor(tomato)TRUE 
##                                  logFC logCPM      F    PValue       FDR
## ExE ectoderm                   -6.5663  13.02 66.267 1.352e-10 3.245e-09
## Mesenchyme                      1.1652  16.29 11.291 1.535e-03 1.841e-02
## Allantois                       0.8345  15.51  5.312 2.555e-02 1.621e-01
## Cardiomyocytes                  0.8484  14.86  5.204 2.701e-02 1.621e-01
## Neural crest                   -0.7706  14.76  4.106 4.830e-02 2.149e-01
## Endothelium                     0.7519  14.29  3.912 5.371e-02 2.149e-01
## Erythroid3                     -0.6431  17.28  3.604 6.367e-02 2.183e-01
## Haematoendothelial progenitors  0.6581  14.72  3.124 8.351e-02 2.505e-01
## ExE mesoderm                    0.3805  15.68  1.181 2.827e-01 6.258e-01
## Pharyngeal mesoderm             0.3793  15.72  1.169 2.850e-01 6.258e-01

注意：

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

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

# pool是哪两个样本属于重复
> y.ab$samples$pool
[1] 3 3 4 4 5 5
# tomato指处理还是对照
> y.ab$samples$tomato
[1]  TRUE FALSE  TRUE FALSE  TRUE FALSE

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

3.2 处理细胞组成的影响

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

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

下面有几种处理方式：

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

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

y.ab2 <- calcNormFactors(y.ab)
y.ab2$samples$norm.factors
y.ab2 <- estimateDisp(y.ab2, design, trend="none")
fit.ab2 <- glmQLFit(y.ab2, design, robust=TRUE, abundance.trend=FALSE)
res2 <- glmQLFTest(fit.ab2, coef=ncol(design))
topTags(res2, n=10)
## Coefficient:  factor(tomato)TRUE 
##                                  logFC logCPM      F    PValue       FDR
## ExE ectoderm                   -6.9215  13.17 70.364 5.738e-11 1.377e-09
## Mesenchyme                      0.9513  16.27  6.787 1.219e-02 1.143e-01
## Neural crest                   -1.0032  14.78  6.464 1.429e-02 1.143e-01
## Erythroid3                     -0.8504  17.35  5.517 2.299e-02 1.380e-01
## Cardiomyocytes                  0.6400  14.84  2.735 1.047e-01 4.809e-01
## Allantois                       0.6054  15.51  2.503 1.202e-01 4.809e-01
## Forebrain/Midbrain/Hindbrain   -0.4943  16.55  1.928 1.713e-01 5.178e-01
## Endothelium                     0.5482  14.27  1.917 1.726e-01 5.178e-01
## Erythroid2                     -0.4818  16.00  1.677 2.015e-01 5.373e-01
## Haematoendothelial progenitors  0.4262  14.73  1.185 2.818e-01 6.240e-01

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

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

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

offenders <- "ExE ectoderm"
# keep.lib.sizes=FALSE参数意为重置文库大小：reset the total number of cells for all samples
y.ab3 <- y.ab[setdiff(rownames(y.ab), offenders),, keep.lib.sizes=FALSE]
y.ab3 <- estimateDisp(y.ab3, design, trend="none")
fit.ab3 <- glmQLFit(y.ab3, design, robust=TRUE, abundance.trend=FALSE)
res3 <- glmQLFTest(fit.ab3, coef=ncol(design))
topTags(res3, n=10)
## Coefficient:  factor(tomato)TRUE 
##                                  logFC logCPM      F   PValue     FDR
## Mesenchyme                      1.1274  16.32 11.501 0.001438 0.03308
## Allantois                       0.7950  15.54  5.231 0.026836 0.18284
## Cardiomyocytes                  0.8104  14.90  5.152 0.027956 0.18284
## Neural crest                   -0.8085  14.80  4.903 0.031798 0.18284
## Erythroid3                     -0.6808  17.32  4.387 0.041743 0.19202
## Endothelium                     0.7151  14.32  3.830 0.056443 0.21636
## Haematoendothelial progenitors  0.6189  14.76  2.993 0.090338 0.29683
## Def. endoderm                   0.4911  12.43  1.084 0.303347 0.67818
## ExE mesoderm                    0.3419  15.71  1.036 0.314058 0.67818
## Pharyngeal mesoderm             0.3407  15.76  1.025 0.316623 0.67818

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

最后更新于4年前