3.9 多样本间差异分析
刘小泽写于2020.7.10+7.15
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个细胞。其中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))
通过这个图可以看到，我们这里存在多对多的关系。也反映出：针对发育分化的细胞，细胞类型注释是个难题。
这也是最最常见的分析思路，一般拿到转录组表达量就会先想着把差异分析做一做 针对单细胞的差异分析方法有多种，一般的综述都会提及。这里使用的方法是“模拟bulk转录组表达量”的模式(Tung et al. 2017)
它的思路是：根据不同细胞类型和不同分群，把对应的细胞表达量加起来。比如根据上面得到的细胞分群信息和注释信息
> 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表达量
- 将多个细胞汇成一组，也不是随便挑的几个细胞，就像上面它是根据两个条件选的同组细胞。就像是生物体内的复制过程，虽然这些细胞之间也存在着一些差异，但相比于组外的细胞，它们之间更相似。而如果直接把每个细胞的表达量传给差异分析工具，它会认为每个细胞就是一个独立的生物样本，这也是和真实情况不符的
差异分析基于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.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
检查一下结果，例如看看每个细胞类型选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"
检测不同条件下细胞类型(或状态等)组成的变化 如果说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个处理进行比较）
首先为列的样本信息增加一些内容：
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 leastmin.countreads in a worthwhile number samples 依据：the filtering keeps genes that have count-per-million (CPM) above k in n samples 定义：k is determined bymin.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放在了前面前面说不再利用
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
最近更新 2yr ago