3.3 挑选表达量高变化基因
刘小泽写于2020.6.28

1 前言

我们经常使用scRNA数据探索细胞异质性。聚类分群和降维也是常用的操作,它们都是依赖于基因表达量(例如综合细胞中各个基因的表达量,然后把表达模式相近的细胞聚在一起)。因此挑选哪些基因进行聚类分群和降维分析,这是非常关键的,挑选的基因一定要有代表性,尽可能多的包含生物信息而剔除其他技术噪音。
单细胞分析的一个重点就是:用尽可能少的维度来展示数据真实的结构。降维的过程其实就是去繁存简,每个基因的变化都对整体有影响,对细胞来讲都是一个变化维度。但这么多维度我们看不了也处理不了,需要尽可能保证真实差异的前提下减少维度的数量,因此需要挑出那些更能代表整体差异的基因。
对每一个细胞来讲,2万多个基因都是它身上长向四面八方的维度,但其中有一些影响很大,一些影响很弱。我们试图把影响最大的几个维度找出来(比如Gene1、Gene2、Gene3),然后用这几个维度代表这个细胞的生物学特征。
上面说了这么多,重点还是:如何挑选有价值的基因(feature)
一般来说,如果一个基因在细胞群体中变化幅度很大,就是受关注对象,我们也会认为是生物因素导致了这么大的差异。这样的基因叫做:高度变化基因(highly variable genes ,HVGs)

数据准备

10X PBMC
其中的处理步骤值得学习,其中pbmc4k_raw_gene_bc_matrices.tar.gz数据已经帮大家下载好,链接:https://share.weiyun.com/s2wbVd7p 密码:3iwypw
下载的内容包含三件套,然后使用read10xCounts读取:
1
# 下载
2
library(BiocFileCache)
3
bfc <- BiocFileCache("raw_data", ask = FALSE)
4
raw.path <- bfcrpath(bfc, file.path("http://cf.10xgenomics.com/samples",
5
"cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"))
6
untar(raw.path, exdir=file.path(tempdir(), "pbmc4k"))
7
8
# 读取
9
suppressMessages(library(DropletUtils))
10
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
11
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
12
> dim(sce.pbmc)
13
[1] 33694 737280
14
15
# 基因注释之ID整合
16
library(scater)
17
rownames(sce.pbmc) <- uniquifyFeatureNames(
18
rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)
19
# 基因注释之添加位置信息
20
library(EnsDb.Hsapiens.v86)
21
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID,
22
column="SEQNAME", keytype="GENEID")
Copied!
其中使用到了一个函数:uniquifyFeatureNames() ,它的作用是
uniquifyFeatureNames() make gene name unique and valid,如果有symbol name的,它会将Ensembl ID替换为symbol name,没有name的仍然使用ID;如果name相同、ID不同(即两个Ensembl ID对应同一个Symbol name),它会将ID与name组合保证特异性(详见?uniquifyFeatureNames) 详见我之前写的:来自刘小泽理解的一份诚意满满的单细胞教程:https://www.jianshu.com/p/e90f5d4d0ab6
之后添加位置信息,得到了基因所在的染色体,可以看到存在13个基因在线粒体上
1
# 接下来
2
# 空液滴检查
3
set.seed(100)
4
e.out <- emptyDrops(counts(sce.pbmc))
5
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]
6
> dim(sce.pbmc)
7
[1] 33694 4233
8
# 可以看到,确实去除了很多空液滴,原来有737280个,留下百分之一不到
9
10
# 质控(尤其关注线粒体)
11
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=which(location=="MT")))
12
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
13
sce.pbmc <- sce.pbmc[,!high.mito]
14
> dim(sce.pbmc)
15
[1] 33694 3922
16
17
# 归一化
18
# 看到quickCluster和computeSumFactors就知道使用的是去卷积化方法
19
library(scran)
20
set.seed(1000)
21
clusters <- quickCluster(sce.pbmc)
22
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
23
sce.pbmc <- logNormCounts(sce.pbmc)
24
sce.pbmc
25
## class: SingleCellExperiment
26
## dim: 33694 3922
27
## metadata(1): Samples
28
## assays(2): counts logcounts
29
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
30
## rowData names(2): ID Symbol
31
## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
32
## TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1
33
## colData names(3): Sample Barcode sizeFactor
34
## reducedDimNames(0):
35
## altExpNames(0):
Copied!
416B数据
1
# 加载数据
2
library(scRNAseq)
3
# sce.416b <- LunSpikeInData(which="416b")
4
sce.416b$block <- factor(sce.416b$block)
5
6
# 基因注释
7
library(AnnotationHub)
8
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
9
rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)
10
rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
11
keytype="GENEID", column="SYMBOL")
12
rowData(sce.416b)$SEQNAME <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
13
keytype="GENEID", column="SEQNAME")
14
15
library(scater)
16
rownames(sce.416b) <- uniquifyFeatureNames(rowData(sce.416b)$ENSEMBL,
17
rowData(sce.416b)$SYMBOL)
18
19
# 质控(添加线粒体、ERCC、批次信息)
20
mito <- which(rowData(sce.416b)$SEQNAME=="MT")
21
stats <- perCellQCMetrics(sce.416b, subsets=list(Mt=mito))
22
qc <- quickPerCellQC(stats, percent_subsets=c("subsets_Mt_percent",
23
"altexps_ERCC_percent"), batch=sce.416b$block)
24
sce.416b <- sce.416b[,!qc$discard]
25
26
# 归一化
27
# 看到computeSumFactors就知道使用的是去卷积化方法
28
library(scran)
29
sce.416b <- computeSumFactors(sce.416b)
30
sce.416b <- logNormCounts(sce.416b)
31
sce.416b
32
## class: SingleCellExperiment
33
## dim: 46604 185
34
## metadata(0):
35
## assays(2): counts logcounts
36
## rownames(46604): 4933401J01Rik Gm26206 ... CAAA01147332.1
37
## CBFB-MYH11-mcherry
38
## rowData names(4): Length ENSEMBL SYMBOL SEQNAME
39
## colnames(185): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
40
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
41
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1
42
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
43
## colData names(10): Source Name cell line ... block sizeFactor
44
## reducedDimNames(0):
45
## altExpNames(2): ERCC SIRV
Copied!
可以看到这里归一化没有使用quickCluster() ,这是因为416b数据自带了block,和cluster的效果一样

2 衡量每个基因的变化程度/方差(Variance)

2.1 方法一:使用log-counts衡量变化程度

这也是最简单的处理方法,基因log后的表达量在细胞间变化幅度越大,就说明细胞间距离越远(这个距离指欧氏距离Euclidean distances)
计算每个基因的变化程度很简单,但要基于这么多基因的结果去挑选最有价值的基因(feature),还需要构建一个模型。这个模型叫:mean-variance relationship,看上去与均值、方差有关
先看计算方法
1
library(scran)
2
dec.pbmc <- modelGeneVar(sce.pbmc)
3
4
# 可视化一条线(下图的蓝线),这条线指所有的基因都会存在的一种偏差
5
fit.pbmc <- metadata(dec.pbmc)
6
plot(fit.pbmc$mean, fit.pbmc$var, xlab="Mean of log-expression",
7
ylab="Variance of log-expression")
8
curve(fit.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)
Copied!
  • 每个点表示一个基因
  • 图中蓝线指的是:技术因素导致的偏差
  • 纵坐标表示总偏差:它等于技术偏差+生物因素偏差
因此,要衡量一个基因的生物因素偏差大小,就看对应的纵坐标减去对应的蓝线的值
最后可以得到一个表:
1
dec.pbmc[order(dec.pbmc$bio, decreasing=TRUE),]
Copied!
可以看到生物因素(bio)导致的偏差从大到小排列,然而会有一些负值存在,难不成总体偏差比技术偏差还要小?这是因为拟合曲线时会有一些基因处于线以下(看上图也能知道),这部分其实可以忽略,反正这部分基因也用不到。我们实则更关心那些bio值更大的

2.2 方法二:使用变异系数(CV)衡量变化程度

这是另一种可选方案,和上面的log-counts是平行关系,它会使用squared coefficient of variation (CV^2)
  • Coefficient of variation (CV)中文翻译是变异系数,它是标准差与均值的比值 ,因此CV^2 就是方差与均值的平方比值
  • 如果两组数据的观测值在一个数量级,那么可以用标准差来比较它们的离散度
  • 如果两个数据差别太大(比如一群老鼠和一群大象的重量),然后就要用CVCV^2,这样可以去除不同量纲均值差异的影响
我们这里可以利用modelGeneCV2()来计算
1
dec.cv2.pbmc <- modelGeneCV2(sce.pbmc)
Copied!
同样的,我们这里的假设是: most genes contain random noise and that the trend captures mostly technical variation,即大部分基因都包含外部噪音,并且拟合的曲线能捕获大部分技术偏差
如果CV^2 越大,偏离拟合曲线越远,就可能包含更多的生物偏差,更具有生物学意义
来自modelGeneCV2()的帮助文档: The ratio of the total CV2 to the trend is used as a metric to rank interesting genes, with larger ratios being indicative of strong biological heterogeneity.
1
fit.cv2.pbmc <- metadata(dec.cv2.pbmc)
2
plot(fit.cv2.pbmc$mean, fit.cv2.pbmc$cv2, log="xy")
3
curve(fit.cv2.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)
Copied!
最后也能得到一个表:
其中偏离拟合曲线的数值用ratio表示,也是我们定义HVGs的途径
这个ratio = total (CV^2) / trend (CV^2) ,这样一相除,就抵消了均值的影响,最后实际上变成了方差的比值。但如果是相减,还是会带着均值【把公式从纸上一列就很清楚了】
1
dec.cv2.pbmc[order(dec.cv2.pbmc$ratio, decreasing=TRUE),]
Copied!
总结:
  • log-counts与CV2都是衡量基因变化程度的有效途径
  • 优先使用log-counts的结果,因为下游的分析也是基于log-counts值

2.3 考虑技术噪音

有两种方式:一种是有spike-in的,一种是没有的,都可以操作
之前画的技术偏差拟合线都是基于一个假设:大部分基因的表达量都是受外界技术误差的影响。但实际上,除了技术误差,还有一些生物因素会导致数据波动(例如transcriptional bursting)。因此之前估计的技术偏差可能”超出了实际的大小“,有些过度估计。
之前估计的偏差实际上是:真正的技术偏差 + 一些无关紧要的生物学因素偏差。因此如果有spike-in的话,它和内源基因不同,不会受到生物因素的影响。利用spike-in画的拟合线,更能代表技术偏差
如果有spike-in
1
dec.spike.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC")
2
dec.spike.416b[order(dec.spike.416b$bio, decreasing=TRUE),]
Copied!
1
# 作图
2
plot(dec.spike.416b$mean, dec.spike.416b$total, xlab="Mean of log-expression",
3
ylab="Variance of log-expression")
4
fit.spike.416b <- metadata(dec.spike.416b)
5
points(fit.spike.416b$mean, fit.spike.416b$var, col="red", pch=16)
6
curve(fit.spike.416b$trend(x), col="dodgerblue", add=TRUE, lwd=2)
Copied!
其中黑点是基因,红点是spike-in,蓝线是根据红点画的拟合线
如果没有spike-in
可以考虑利用数据分布来表示技术噪音,例如只考虑技术噪音的话,UMI counts通常会呈现近似泊松分布
1
set.seed(0010101)
2
dec.pois.pbmc <- modelGeneVarByPoisson(sce.pbmc)
3
dec.pois.pbmc <- dec.pois.pbmc[order(dec.pois.pbmc$bio, decreasing=TRUE),]
4
head(dec.pois.pbmc)
5
## DataFrame with 6 rows and 6 columns
6
## mean total tech bio p.value FDR
7
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
8
## LYZ 1.97770 5.11595 0.621547 4.49440 0 0
9
## S100A9 1.94951 4.58859 0.627306 3.96128 0 0
10
## S100A8 1.71828 4.45723 0.669428 3.78781 0 0
11
## HLA-DRA 2.09694 3.72690 0.596372 3.13053 0 0
12
## CD74 2.89840 3.30912 0.422624 2.88650 0 0
13
## CST3 1.49285 2.97369 0.695367 2.27833 0 0
14
15
# 作图
16
plot(dec.pois.pbmc$mean, dec.pois.pbmc$total, pch=16, xlab="Mean of log-expression",
17
ylab="Variance of log-expression")
18
curve(metadata(dec.pois.pbmc)$trend(x), col="dodgerblue", add=TRUE)
Copied!
有趣的是,如果真的单纯考虑技术偏差,我们会得到比原来更多的变化基因。其中就包含了之前没有的一些管家基因。不过这些基因在后面分析(比如研究细胞异质性)中又不太感兴趣,因此使用更准确的技术噪音拟合方法,可能不会得到更准确的HVGs排序

2.4 考虑批次效应

2.4.1 方法一:绘制拟合曲线时加入批次信息(Fitting block-specific trends)
有时我们看到的在不同细胞间高变化的基因,不一定是生物学因素驱动,也不一定是由于技术误差导致,而是由于不同细胞的批次不同,导致同一基因的表达量差异很大。
我们更想知道的是,在一个批次中的HVGs信息。一般可以一个批次一个批次地去处理,例如:
1
# 这里用spike-in来拟合更准确的技术偏差,同时考虑上批次信息
2
dec.block.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC", block=sce.416b$block)
3
head(dec.block.416b[order(dec.block.416b$bio, decreasing=TRUE),1:6])
4
## DataFrame with 6 rows and 6 columns
5
## mean total tech bio p.value FDR
6
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
7
## Lyz2 6.61235 13.8619 1.58416 12.2777 0.00000e+00 0.00000e+00
8
## Ccl9 6.67841 13.2599 1.44553 11.8143 0.00000e+00 0.00000e+00
9
## Top2a 5.81275 14.0192 2.74571 11.2734 3.89855e-137 8.43398e-135
10
## Cd200r3 4.83305 15.5909 4.31892 11.2719 1.17783e-54 7.00721e-53
11
## Ccnb2 5.97999 13.0256 2.46647 10.5591 1.20380e-151 2.98405e-149
12
## Hbb-bt 4.91683 14.6539 4.12156 10.5323 2.52639e-49 1.34197e-47
Copied!
看看不同批次的差异
1
par(mfrow=c(1,2))
2
blocked.stats <- dec.block.416b$per.block
3
for (i in colnames(blocked.stats)) {
4
current <- blocked.stats[[i]]
5
plot(current$mean, current$total, main=i, pch=16, cex=0.5,
6
xlab="Mean of log-expression", ylab="Variance of log-expression")
7
curfit <- metadata(current)
8
points(curfit$mean, curfit$var, col="red", pch=16)
9
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
10
}
Copied!
其中黑点是基因,红点是spike-in,蓝线是根据红点画的拟合线。看到这里的两个批次之间差异很小,表示重复效果不错
2.4.2 利用实验设计矩阵(design matrix)
当只有一个批次信息时(例如只考虑sce.416b$block信息),使用上面的绘制拟合曲线方法是足够的,但如果再加上其他的批次信息,例如:
1
> table(colData(sce.416b)$phenotype)
2
3
induced CBFB-MYH11 oncogene expression
4
96
5
wild type phenotype
6
96
Copied!
就需要利用实验设计矩阵(很像DESeq2、edgeR等差异分析做的),像是modelGeneVarWithSpikes()modelGeneVar() 都支持这个参数
1
design <- model.matrix(~factor(block) + phenotype, colData(sce.416b))
2
dec.design.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC", design=design)
3
dec.design.416b[order(dec.design.416b$bio, decreasing=TRUE),]
Copied!
这种方法做起来很简单,但不太准确,因为没有考虑到各个批次中的平均表达水平。而这个表达水平又影响真正的技术偏差估计
the true technical component is the average of the fitted values at the per-block means, which may be quite different for strong batch effects and non-linear mean-variance relationships
相比之下,第一种方法的block= 更稳妥一些

3 挑选出高变化基因HVGs

既然前面找到了每个基因的变化程度,并且排了个序。但挑多少基因呢,这是个问题。我们首先想是不是多多益善?
的确,一个大的基因集保留了更多的基因信息,可以防止丢失一些有趣的生物信号,但代价是增加一些和研究不相干的基因信号,造成干扰。而哪些是不相干,这个不好权衡,因为一种情况下的不相干,可能在另一种情况下又是有用的。
例如,T细胞活化反应中存在的异质性是有趣的,但如果我们关注点是区分主要的免疫表型,那么前面这个异质性就是无关的噪音。简言之,虽然都是有趣的生物学信号,但我们关注点只有一个,研究一个势必另一个会成为干扰因素。
关于HVGs的选择,有以下几种方案:

3.1 基于方差最大化

这个方法是最简单理解的,就是提取在生物学因素中方差最大的(也就是对波动贡献最大)基因。优点就是可以自己控制基因的数量,方便对后面的分析复杂度有一个大体的估计。
比如可以挑选前1000个变化最大的基因:
对于modelGeneVar()modelGeneVarWithSpikes() 基于方差的,直接根据结果表格bio这一列提取
1
hvg.pbmc.var <- getTopHVGs(dec.pbmc, n=1000)
2
str(hvg.pbmc.var)
3
## chr [1:1000] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" ...
Copied!
对于modelGeneCV2()modelGeneCV2WithSpikes() 基于CV2的,可以根据结果表格的ratio这一列提取
1
hvg.pbmc.cv2 <- getTopHVGs(dec.cv2.pbmc, var.field="ratio", n=1000)
2
str(hvg.pbmc.cv2)
3
## chr [1:1000] "HIST1H2AC" "GNG11" "PRTFDC1" "TNNC2" "PF4" "HGD" "PPBP" ...
Copied!
重点还是参数n的设定,实际上大部分基因对生物学异质性的贡献不大。如果我们知道一个假设:不超过5%的基因是差异表达的,那么这个问题就轻松许多。我们之间设置n小于等于基- 因总数的5%即可。但事实上,我们实现不知道差异基因的占比,这个5%也只是虚构。但给我们一些提示:
  • 如果数据异质性比较低,可以降低n的数值,保留最多的生物信号同时,尽可能剔除不相干的噪音
  • 如果数据异型性很高,就应该使用更大的n值,在寻求主要差异的同时,保留更多的次要差异因素
从上面的操作和解读来看,这个方法存在一定的弊端:设置一个固定的数值,就让基因们形成了竞争关系,能被n选择成为HVGs的,势必会把其他有信息价值的基因排除在外,如果这些被排除的基因还是一些细胞亚群的marker基因【这种问题在高度异质性群体中尤为严重】
另外这个n的选择很主观,从500-5000之间似乎都能说得通。上面选择1000,其实也没有什么特别的原因,只是一个直觉。如果选择这种方法,那么什么都不要管,根据直觉选择一个值,进行剩余的分析,通过结果来判断这个值选的是不是合理,而不是在这里辗转反侧思考用什么数值。

3.2 基于假设检验

原假设是:基因的方差与之前拟合的曲线相符(也就是对生物学角度的波动没什么贡献)。由此可以设置adjusted p-values的阈值,例如
1
hvg.pbmc.var.2 <- getTopHVGs(dec.pbmc, fdr.threshold=0.05)
2
length(hvg.pbmc.var.2)
3
## [1] 651
Copied!
这个方法适用于HVGs的确大部分都是与研究点相关,使用FDR只会让结果更准确;但缺点是:比第一种指定数量的方法更难以预测。既然是假设检验,就可能存在假阳性,可能得到更多不相干的基因。另外,我们的关注重点实际上不在于这些基因本身,而是看它们对下游分析有没有作用,因此也不能认为5%的FDR阈值可以得到低噪音的基因集。
或许还会想,之前计算的表格中有p值这一列,那么可不可以用这一列来替代bio或者ratio那一列进行排序,从而挑选基因呢?
最好不要,因为还是那句话,表现显著的基因不一定对我们感兴趣的生物学问题有帮助。很多基因确实p值很小,这是因为它们的技术偏差非常小,但同时总体偏差并不大,因此它看上去很显著,但实际我们不感兴趣。
因此这种方法还是慎用

3.3 保留之前拟合曲线上面全部的基因

这个方法的目的很单纯,拟合曲线以下的基因都是不感兴趣的,因为它们的生物学偏差很小,于是把它们全部剔除,留下剩下的。这也是走了一个极端:为了偏差最小化,选择了噪音最大化(保留了很多与研究问题无关的基因)。看下面计算的数量就知道
如果之前使用了modelGeneVar()
1
hvg.pbmc.var.3 <- getTopHVGs(dec.pbmc, var.threshold=0)
2
length(hvg.pbmc.var.3)
3
## [1] 12791
Copied!
如果之前使用了modelGeneCV2():保留所有ratio大于1的
1
hvg.pbmc.cv2.3 <- getTopHVGs(dec.cv2.pbmc, var.field="ratio", var.threshold=1)
2
length(hvg.pbmc.cv2.3)
3
## [1] 9295
Copied!
这种方法对于研究罕见的细胞亚群最有帮助,另外对于异质性非常严重的数据(包含许多不同类型的细胞类型,例如多个数据集合并后的数据)比较适合

4 特殊情况一:Keep or not

场景就是:我们筛选出来HVGs,剩下的基因是删除还是继续保留呢? 主要有三种方法

首先选取前10%的HVGs

利用参数prop=0.1
1
dec.pbmc <- modelGeneVar(sce.pbmc)
2
chosen <- getTopHVGs(dec.pbmc, prop=0.1)
3
str(chosen)
4
## chr [1:1279] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" ...
Copied!

方法一:直接删除非HVGs

下游分析也全部基于HVGs,不过日后如果突然发现,一些非HVGs的基因也有点意思,想看看它们就变得很麻烦
1
sce.pbmc.hvg <- sce.pbmc[chosen,]
2
dim(sce.pbmc.hvg)
3
## [1] 1279 3922
Copied!

方法二:原数据不动,只是下游分析基于小部分HVGs

这样的好处是,如果以后调整了chosen这个HVGs基因集,那么还是可以继续进行下游分析的
1
# 例如只根据HVGs进行PCA,使用参数subset_row
2
library(scater)
3
sce.pbmc <- runPCA(sce.pbmc, subset_row=chosen)
4
reducedDimNames(sce.pbmc)
5
## [1] "PCA"
Copied!

方法三:原数据先备份,下游分析也是基于原数据,最后再复原

1
# 首先保留一份原数据
2
altExp(sce.pbmc.hvg, "original") <- sce.pbmc
3
altExpNames(sce.pbmc.hvg)
4
## [1] "original"
5
6
# 然后按方法一得到处理的数据,并进行分析(这里就不需要指定参数subset_row)
7
sce.pbmc.hvg <- runPCA(sce.pbmc.hvg)
8
9
# 最后复原
10
sce.pbmc.original <- altExp(sce.pbmc.hvg, "original", withColData=TRUE)
Copied!

5 特殊情况二:利用先验基因

如果想证明某个通路的基因不存在有意义的异质性,那么可以选取这部分基因作为先验基因,重复一遍分析,检查到底是否存在异质性。这个想法和荧光激活细胞分选(fluorescence activated cell sorting ,FACS)相似,只不多scRNA分析中可以更方便地更改先验基因

选取部分先验基因

例如PBMC数据集中,可以用MSigDB数据库的免疫相关基因集:C7 immunologic signatures
1
# 首先选取基因集
2
library(msigdbr)
3
c7.sets <- msigdbr(species = "Homo sapiens", category = "C7")
4
head(unique(c7.sets$gs_name))
5
## [1] "GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_DN"
6
## [2] "GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_UP"
7
## [3] "GOLDRATH_NAIVE_VS_EFF_CD8_TCELL_DN"
8
## [4] "GOLDRATH_NAIVE_VS_EFF_CD8_TCELL_UP"
9
## [5] "GOLDRATH_NAIVE_VS_MEMORY_CD8_TCELL_DN"
10
## [6] "GOLDRATH_NAIVE_VS_MEMORY_CD8_TCELL_UP"
11
12
# 利用Goldrath基因集来区分CD8细胞亚型
13
cd8.sets <- c7.sets[grep("GOLDRATH", c7.sets$gs_name),]
14
cd8.genes <- rowData(sce.pbmc)$Symbol %in% cd8.sets$human_gene_symbol
15
summary(cd8.genes)
16
## Mode FALSE TRUE
17
## logical 32869 825
18
19
# 利用GSE11924基因集来区分T helper细胞亚型
20
th.sets <- c7.sets[grep("GSE11924", c7.sets$gs_name),]
21
th.genes <- rowData(sce.pbmc)$Symbol %in% th.sets$human_gene_symbol
22
summary(th.genes)
23
## Mode FALSE TRUE
24
## logical 31786 1908
25
26
# 利用GSE11961基因集来区分B细胞亚型
27
b.sets <- c7.sets[grep("GSE11961", c7.sets$gs_name),]
28
b.genes <- rowData(sce.pbmc)$Symbol %in% b.sets$human_gene_symbol
29
summary(b.genes)
30
## Mode FALSE TRUE
31
## logical 28185 5509
Copied!
以上就是我们自己筛选的先验基因集,然后利用这些作为chosen基因再进行下游分析,看看结果与我们推测的是否一致

剔除某些先验基因

和上面相反,我们还可以先剔除一些比较不感兴趣的基因,再进行下游分析
例如,可以剔除核糖体蛋白基因或线粒体基因
1
# 匹配核糖体蛋白基因
2
ribo.discard <- grepl("^RP[SL]\\d+", rownames(sce.pbmc))
3
sum(ribo.discard)
4
## [1] 99
5
6
# 另一种更准确的方法
7
c2.sets <- msigdbr(species = "Homo sapiens", category = "C2")
8
ribo.set <- c2.sets[c2.sets$gs_name=="KEGG_RIBOSOME",]$human_gene_symbol
9
ribo.discard <- rownames(sce.pbmc) %in% ribo.set
10
sum(ribo.discard)
11
## [1] 87
Copied!
对于免疫细胞数据集,可以剔除免疫球蛋白基因和T细胞受体基因【当然这些都取决于生物背景】
1
library(AnnotationHub)
2
edb <- AnnotationHub()[["AH73881"]]
3
anno <- select(edb, keys=rowData(sce.pbmc)$ID, keytype="GENEID",
4
columns="TXBIOTYPE")
5
6
# 去除免疫球蛋白基因
7
igv.set <- anno$GENEID[anno$TXBIOTYPE %in% c("IG_V_gene", "IG_V_pseudogene")]
8
igv.discard <- rowData(sce.pbmc)$ID %in% igv.set
9
sum(igv.discard)
10
## [1] 326
11
12
# 去除T细胞受体基因
13
tcr.set <- anno$GENEID[anno$TXBIOTYPE %in% c("TR_V_gene", "TR_V_pseudogene")]
14
tcr.discard <- rowData(sce.pbmc)$ID %in% tcr.set
15
sum(tcr.discard)
16
## [1] 138
Copied!
不过以上的操作要慎重,如果下游分析没有发现问题,以及对先验基因存在十足的把握,还是不要随意使用先验基因
最近更新 1yr ago