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读取：
复制 # 下载
library (BiocFileCache)
bfc <- BiocFileCache( "raw_data" , ask = FALSE )
raw.path <- bfcrpath( bfc, file.path ( "http://cf.10xgenomics.com/samples" ,
"cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz" ) )
untar (raw.path, exdir = file.path ( tempdir (), "pbmc4k" ))
# 读取
suppressMessages ( library (DropletUtils))
fname <- file.path ( tempdir (), "pbmc4k/raw_gene_bc_matrices/GRCh38" )
sce.pbmc <- read10xCounts( fname, col.names = TRUE )
> dim (sce.pbmc)
[ 1 ] 33694 737280
# 基因注释之ID整合
library (scater)
rownames (sce.pbmc) <- uniquifyFeatureNames(
rowData( sce.pbmc ) $ ID, rowData( sce.pbmc ) $ Symbol )
# 基因注释之添加位置信息
library (EnsDb.Hsapiens.v86)
location <- mapIds( EnsDb.Hsapiens.v86, keys = rowData( sce.pbmc ) $ ID,
column = "SEQNAME" , keytype = "GENEID" ) 其中使用到了一个函数：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个基因在线粒体上
复制 # 接下来
# 空液滴检查
set.seed ( 100 )
e.out <- emptyDrops(counts( sce.pbmc ))
sce.pbmc <- sce.pbmc[, which (e.out $ FDR <= 0.001 )]
> dim (sce.pbmc)
[ 1 ] 33694 4233
# 可以看到，确实去除了很多空液滴，原来有737280个，留下百分之一不到
# 质控（尤其关注线粒体）
stats <- perCellQCMetrics( sce.pbmc, subsets =list (Mito = which (location == "MT" ) ) )
high.mito <- isOutlier( stats $ subsets_Mito_percent, type = "higher" )
sce.pbmc <- sce.pbmc[, ! high.mito]
> dim (sce.pbmc)
[ 1 ] 33694 3922
# 归一化
# 看到quickCluster和computeSumFactors就知道使用的是去卷积化方法
library (scran)
set.seed ( 1000 )
clusters <- quickCluster( sce.pbmc )
sce.pbmc <- computeSumFactors( sce.pbmc, cluster = clusters )
sce.pbmc <- logNormCounts( sce.pbmc )
sce.pbmc
## class: SingleCellExperiment
## dim: 33694 3922
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
## rowData names(2): ID Symbol
## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
## TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1
## colData names(3): Sample Barcode sizeFactor
## reducedDimNames(0):
## altExpNames(0): 416B数据
复制 # 加载数据
library (scRNAseq)
# sce.416b <- LunSpikeInData(which="416b")
sce.416b $ block <- factor (sce.416b $ block)
# 基因注释
library (AnnotationHub)
ens.mm.v97 <- AnnotationHub() [[ "AH73905" ]]
rowData( sce.416b ) $ ENSEMBL <- rownames (sce.416b)
rowData( sce.416b ) $ SYMBOL <- mapIds( ens.mm.v97, keys = rownames (sce.416b ) ,
keytype = "GENEID" , column = "SYMBOL" )
rowData( sce.416b ) $ SEQNAME <- mapIds( ens.mm.v97, keys = rownames (sce.416b ) ,
keytype = "GENEID" , column = "SEQNAME" )
library (scater)
rownames (sce.416b) <- uniquifyFeatureNames(rowData( sce.416b ) $ ENSEMBL,
rowData( sce.416b ) $ SYMBOL )
# 质控（添加线粒体、ERCC、批次信息）
mito <- which ( rowData( sce.416b ) $ SEQNAME == "MT" )
stats <- perCellQCMetrics( sce.416b, subsets =list (Mt = mito) )
qc <- quickPerCellQC( stats, percent_subsets = c ( "subsets_Mt_percent" ,
"altexps_ERCC_percent" ) , batch = sce.416b $ block)
sce.416b <- sce.416b[, ! qc $ discard]
# 归一化
# 看到computeSumFactors就知道使用的是去卷积化方法
library (scran)
sce.416b <- computeSumFactors( sce.416b )
sce.416b <- logNormCounts( sce.416b )
sce.416b
## class: SingleCellExperiment
## dim: 46604 185
## metadata(0):
## assays(2): counts logcounts
## rownames(46604): 4933401J01Rik Gm26206 ... CAAA01147332.1
## CBFB-MYH11-mcherry
## rowData names(4): Length ENSEMBL SYMBOL SEQNAME
## colnames(185): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
## colData names(10): Source Name cell line ... block sizeFactor
## reducedDimNames(0):
## altExpNames(2): ERCC SIRV 可以看到这里归一化没有使用quickCluster() ，这是因为416b数据自带了block，和cluster的效果一样
2 衡量每个基因的变化程度/方差(Variance)
2.1 方法一：使用log-counts衡量变化程度
这也是最简单的处理方法，基因log后的表达量在细胞间变化幅度越大，就说明细胞间距离越远（这个距离指欧氏距离Euclidean distances）
计算每个基因的变化程度很简单，但要基于这么多基因的结果去挑选最有价值的基因（feature），还需要构建一个模型。这个模型叫：mean-variance relationship，看上去与均值、方差有关
先看计算方法
复制 library (scran)
dec.pbmc <- modelGeneVar( sce.pbmc )
# 可视化一条线（下图的蓝线），这条线指所有的基因都会存在的一种偏差
fit.pbmc <- metadata( dec.pbmc )
plot (fit.pbmc $ mean, fit.pbmc $ var, xlab = "Mean of log-expression" ,
ylab = "Variance of log-expression" )
curve (fit.pbmc $ trend( x ) , col = "dodgerblue" , add = TRUE , lwd = 2 ) 因此，要衡量一个基因的生物因素偏差大小，就看对应的纵坐标减去对应的蓝线的值
最后可以得到一个表：
复制 dec.pbmc[ order (dec.pbmc $ bio, decreasing = TRUE ),] 可以看到生物因素（bio）导致的偏差从大到小排列，然而会有一些负值存在，难不成总体偏差比技术偏差还要小？这是因为拟合曲线时会有一些基因处于线以下（看上图也能知道），这部分其实可以忽略，反正这部分基因也用不到。我们实则更关心那些bio值更大的
2.2 方法二：使用变异系数(CV)衡量变化程度
这是另一种可选方案，和上面的log-counts是平行关系，它会使用squared coefficient of variation (CV^2)
关于概念： 可以看https://www.jianshu.com/p/f2df679843f3 和 https://www.jianshu.com/p/3525e624946a
Coefficient of variation （CV）中文翻译是变异系数 ，它是标准差与均值的比值 ，因此CV^2 就是方差与均值的平方比值
如果两组数据的观测值在一个数量级，那么可以用标准差来比较它们的离散度
如果两个数据差别太大（比如一群老鼠和一群大象的重量），然后就要用CV或CV^2，这样可以去除不同量纲均值差异的影响
我们这里可以利用modelGeneCV2()来计算
复制 dec.cv2.pbmc <- modelGeneCV2( sce.pbmc ) 同样的，我们这里的假设是： 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.
复制 fit.cv2.pbmc <- metadata( dec.cv2.pbmc )
plot (fit.cv2.pbmc $ mean, fit.cv2.pbmc $ cv2, log = "xy" )
curve (fit.cv2.pbmc $ trend( x ) , col = "dodgerblue" , add = TRUE , lwd = 2 ) 最后也能得到一个表：
其中偏离拟合曲线的数值用ratio表示，也是我们定义HVGs的途径
这个ratio = total (CV^2) / trend (CV^2) ，这样一相除，就抵消了均值的影响，最后实际上变成了方差的比值。但如果是相减，还是会带着均值【把公式从纸上一列就很清楚了】
复制 dec.cv2.pbmc[ order (dec.cv2.pbmc $ ratio, decreasing = TRUE ),] 总结：
log-counts与CV2都是衡量基因变化程度的有效途径
优先使用log-counts的结果，因为下游的分析也是基于log-counts值
2.3 考虑技术噪音
有两种方式：一种是有spike-in的，一种是没有的，都可以操作
之前画的技术偏差拟合线都是基于一个假设：大部分基因的表达量都是受外界技术误差的影响。但实际上，除了技术误差，还有一些生物因素会导致数据波动（例如transcriptional bursting）。因此之前估计的技术偏差可能”超出了实际的大小“，有些过度估计。
之前估计的偏差实际上是：真正的技术偏差 + 一些无关紧要的生物学因素偏差。因此如果有spike-in的话，它和内源基因不同，不会受到生物因素的影响。利用spike-in画的拟合线，更能代表技术偏差
如果有spike-in
复制 dec.spike.416b <- modelGeneVarWithSpikes( sce.416b, "ERCC" )
dec.spike.416b[ order (dec.spike.416b $ bio, decreasing = TRUE ),]
复制 # 作图
plot (dec.spike.416b $ mean, dec.spike.416b $ total, xlab = "Mean of log-expression" ,
ylab = "Variance of log-expression" )
fit.spike.416b <- metadata( dec.spike.416b )
points (fit.spike.416b $ mean, fit.spike.416b $ var, col = "red" , pch = 16 )
curve (fit.spike.416b $ trend( x ) , col = "dodgerblue" , add = TRUE , lwd = 2 ) 其中黑点是基因，红点是spike-in，蓝线是根据红点画的拟合线
如果没有spike-in
可以考虑利用数据分布来表示技术噪音，例如只考虑技术噪音的话，UMI counts通常会呈现近似泊松分布
复制 set.seed ( 0010101 )
dec.pois.pbmc <- modelGeneVarByPoisson( sce.pbmc )
dec.pois.pbmc <- dec.pois.pbmc[ order (dec.pois.pbmc $ bio, decreasing = TRUE ),]
head (dec.pois.pbmc)
## DataFrame with 6 rows and 6 columns
## mean total tech bio p.value FDR
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## LYZ 1.97770 5.11595 0.621547 4.49440 0 0
## S100A9 1.94951 4.58859 0.627306 3.96128 0 0
## S100A8 1.71828 4.45723 0.669428 3.78781 0 0
## HLA-DRA 2.09694 3.72690 0.596372 3.13053 0 0
## CD74 2.89840 3.30912 0.422624 2.88650 0 0
## CST3 1.49285 2.97369 0.695367 2.27833 0 0
# 作图
plot (dec.pois.pbmc $ mean, dec.pois.pbmc $ total, pch = 16 , xlab = "Mean of log-expression" ,
ylab = "Variance of log-expression" )
curve ( metadata( dec.pois.pbmc ) $ trend( x ) , col = "dodgerblue" , add = TRUE ) 有趣的是，如果真的单纯考虑技术偏差，我们会得到比原来更多的变化基因。其中就包含了之前没有的一些管家基因。不过这些基因在后面分析（比如研究细胞异质性）中又不太感兴趣，因此使用更准确的技术噪音拟合方法，可能不会得到更准确的HVGs排序
2.4 考虑批次效应
2.4.1 方法一：绘制拟合曲线时加入批次信息（Fitting block-specific trends）
有时我们看到的在不同细胞间高变化的基因，不一定是生物学因素驱动，也不一定是由于技术误差导致，而是由于不同细胞的批次不同，导致同一基因的表达量差异很大。
我们更想知道的是，在一个批次中的HVGs信息。一般可以一个批次一个批次地去处理，例如：
复制 # 这里用spike-in来拟合更准确的技术偏差，同时考虑上批次信息
dec.block.416b <- modelGeneVarWithSpikes( sce.416b, "ERCC" , block = sce.416b $ block )
head (dec.block.416b[ order (dec.block.416b $ bio, decreasing = TRUE ), 1 : 6 ])
## DataFrame with 6 rows and 6 columns
## mean total tech bio p.value FDR
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## Lyz2 6.61235 13.8619 1.58416 12.2777 0.00000e+00 0.00000e+00
## Ccl9 6.67841 13.2599 1.44553 11.8143 0.00000e+00 0.00000e+00
## Top2a 5.81275 14.0192 2.74571 11.2734 3.89855e-137 8.43398e-135
## Cd200r3 4.83305 15.5909 4.31892 11.2719 1.17783e-54 7.00721e-53
## Ccnb2 5.97999 13.0256 2.46647 10.5591 1.20380e-151 2.98405e-149
## Hbb-bt 4.91683 14.6539 4.12156 10.5323 2.52639e-49 1.34197e-47 看看不同批次的差异
复制 par (mfrow = c ( 1 , 2 ))
blocked.stats <- dec.block.416b $ per.block
for (i in colnames (blocked.stats)) {
current <- blocked.stats[[i]]
plot (current $ mean, current $ total, main = i, pch = 16 , cex = 0.5 ,
xlab = "Mean of log-expression" , ylab = "Variance of log-expression" )
curfit <- metadata( current )
points (curfit $ mean, curfit $ var, col = "red" , pch = 16 )
curve (curfit $ trend( x ) , col = 'dodgerblue' , add = TRUE , lwd = 2 )
} 其中黑点是基因，红点是spike-in，蓝线是根据红点画的拟合线。看到这里的两个批次之间差异很小，表示重复效果不错
2.4.2 利用实验设计矩阵（design matrix）
当只有一个批次信息时（例如只考虑sce.416b$block信息），使用上面的绘制拟合曲线方法是足够的，但如果再加上其他的批次信息，例如：
复制 > table ( colData( sce.416b ) $ phenotype)
induced CBFB - MYH11 oncogene expression
96
wild type phenotype
96 就需要利用实验设计矩阵（很像DESeq2、edgeR等差异分析做的），像是modelGeneVarWithSpikes()和modelGeneVar() 都支持这个参数
复制 design <- model.matrix ( ~ factor (block) + phenotype, colData( sce.416b ) )
dec.design.416b <- modelGeneVarWithSpikes( sce.416b, "ERCC" , design = design )
dec.design.416b[ order (dec.design.416b $ bio, decreasing = TRUE ),] 这种方法做起来很简单，但不太准确，因为没有考虑到各个批次中的平均表达水平。而这个表达水平又影响真正的技术偏差估计
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这一列提取
复制 hvg.pbmc.var <- getTopHVGs( dec.pbmc, n = 1000 )
str (hvg.pbmc.var)
## chr [1:1000] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" ... 对于modelGeneCV2()和modelGeneCV2WithSpikes() 基于CV2的，可以根据结果表格的ratio这一列提取
复制 hvg.pbmc.cv2 <- getTopHVGs( dec.cv2.pbmc, var.field = "ratio" , n = 1000 )
str (hvg.pbmc.cv2)
## chr [1:1000] "HIST1H2AC" "GNG11" "PRTFDC1" "TNNC2" "PF4" "HGD" "PPBP" ... 重点还是参数 n 的设定 ，实际上大部分基因对生物学异质性的贡献不大。如果我们知道一个假设：不超过5%的基因是差异表达的，那么这个问题就轻松许多。我们之间设置n小于等于基- 因总数的5%即可。但事实上，我们实现不知道差异基因的占比，这个5%也只是虚构。但给我们一些提示：
如果数据异质性比较低，可以降低n的数值，保留最多的生物信号同时，尽可能剔除不相干的噪音
如果数据异型性很高，就应该使用更大的n值，在寻求主要差异的同时，保留更多的次要差异因素
从上面的操作和解读来看，这个方法存在一定的弊端 ：设置一个固定的数值，就让基因们形成了竞争关系，能被n选择成为HVGs的，势必会把其他有信息价值的基因排除在外，如果这些被排除的基因还是一些细胞亚群的marker基因【这种问题在高度异质性群体中尤为严重】
另外这个n的选择很主观，从500-5000之间似乎都能说得通。上面选择1000，其实也没有什么特别的原因，只是一个直觉。如果选择这种方法，那么什么都不要管，根据直觉选择一个值 ，进行剩余的分析，通过结果来判断这个值选的是不是合理，而不是在这里辗转反侧思考用什么数值。
3.2 基于假设检验
原假设是：基因的方差与之前拟合的曲线相符（也就是对生物学角度的波动没什么贡献）。由此可以设置adjusted p-values的阈值，例如
复制 hvg.pbmc.var .2 <- getTopHVGs( dec.pbmc, fdr.threshold = 0.05 )
length (hvg.pbmc.var .2 )
## [1] 651 这个方法适用于HVGs的确大部分都是与研究点相关，使用FDR只会让结果更准确；但缺点是：比第一种指定数量的方法更难以预测。既然是假设检验，就可能存在假阳性，可能得到更多不相干的基因。另外，我们的关注重点实际上不在于这些基因本身，而是看它们对下游分析有没有作用，因此也不能认为5%的FDR阈值可以得到低噪音的基因集。
或许还会想，之前计算的表格中有p值这一列，那么可不可以用这一列来替代bio或者ratio那一列进行排序，从而挑选基因呢？
最好不要，因为还是那句话，表现显著的基因不一定对我们感兴趣的生物学问题有帮助。 很多基因确实p值很小，这是因为它们的技术偏差非常小，但同时总体偏差并不大，因此它看上去很显著，但实际我们不感兴趣。
因此这种方法还是慎用
3.3 保留之前拟合曲线上面全部的基因
这个方法的目的很单纯，拟合曲线以下的基因都是不感兴趣的，因为它们的生物学偏差很小，于是把它们全部剔除，留下剩下的。这也是走了一个极端：为了偏差最小化，选择了噪音最大化（保留了很多与研究问题无关的基因）。看下面计算的数量就知道
如果之前使用了modelGeneVar()：
复制 hvg.pbmc.var .3 <- getTopHVGs( dec.pbmc, var.threshold = 0 )
length (hvg.pbmc.var .3 )
## [1] 12791 如果之前使用了modelGeneCV2()：保留所有ratio大于1的
复制 hvg.pbmc.cv2 .3 <- getTopHVGs( dec.cv2.pbmc, var.field = "ratio" , var.threshold = 1 )
length (hvg.pbmc.cv2 .3 )
## [1] 9295 这种方法对于研究罕见的细胞亚群最有帮助，另外对于异质性非常严重的数据（包含许多不同类型的细胞类型，例如多个数据集合并后的数据）比较适合
4 特殊情况一：Keep or not
场景就是：我们筛选出来HVGs，剩下的基因是删除还是继续保留呢？ 主要有三种方法
首先选取前10%的HVGs
利用参数prop=0.1
复制 dec.pbmc <- modelGeneVar( sce.pbmc )
chosen <- getTopHVGs( dec.pbmc, prop = 0.1 )
str (chosen)
## chr [1:1279] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" ... 方法一：直接删除非HVGs
下游分析也全部基于HVGs，不过日后如果突然发现，一些非HVGs的基因也有点意思，想看看它们就变得很麻烦
复制 sce.pbmc.hvg <- sce.pbmc[chosen,]
dim (sce.pbmc.hvg)
## [1] 1279 3922 方法二：原数据不动，只是下游分析基于小部分HVGs
这样的好处是，如果以后调整了chosen这个HVGs基因集，那么还是可以继续进行下游分析的
复制 # 例如只根据HVGs进行PCA，使用参数subset_row
library (scater)
sce.pbmc <- runPCA( sce.pbmc, subset_row = chosen )
reducedDimNames( sce.pbmc )
## [1] "PCA" 方法三：原数据先备份，下游分析也是基于原数据，最后再复原
复制 # 首先保留一份原数据
altExp( sce.pbmc.hvg, "original" ) <- sce.pbmc
altExpNames( sce.pbmc.hvg )
## [1] "original"
# 然后按方法一得到处理的数据，并进行分析（这里就不需要指定参数subset_row）
sce.pbmc.hvg <- runPCA( sce.pbmc.hvg )
# 最后复原
sce.pbmc.original <- altExp( sce.pbmc.hvg, "original" , withColData = TRUE ) 5 特殊情况二：利用先验基因
如果想证明某个通路的基因不存在有意义的异质性，那么可以选取这部分基因作为先验基因，重复一遍分析，检查到底是否存在异质性。这个想法和荧光激活细胞分选（fluorescence activated cell sorting ，FACS）相似，只不多scRNA分析中可以更方便地更改先验基因
选取部分先验基因
例如PBMC数据集中，可以用MSigDB数据库的免疫相关基因集：C7 immunologic signatures
复制 # 首先选取基因集
library (msigdbr)
c7.sets <- msigdbr( species = "Homo sapiens" , category = "C7" )
head ( unique (c7.sets $ gs_name))
## [1] "GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_DN"
## [2] "GOLDRATH_EFF_VS_MEMORY_CD8_TCELL_UP"
## [3] "GOLDRATH_NAIVE_VS_EFF_CD8_TCELL_DN"
## [4] "GOLDRATH_NAIVE_VS_EFF_CD8_TCELL_UP"
## [5] "GOLDRATH_NAIVE_VS_MEMORY_CD8_TCELL_DN"
## [6] "GOLDRATH_NAIVE_VS_MEMORY_CD8_TCELL_UP"
# 利用Goldrath基因集来区分CD8细胞亚型
cd8.sets <- c7.sets[ grep ( "GOLDRATH" , c7.sets $ gs_name),]
cd8.genes <- rowData( sce.pbmc ) $ Symbol %in% cd8.sets $ human_gene_symbol
summary (cd8.genes)
## Mode FALSE TRUE
## logical 32869 825
# 利用GSE11924基因集来区分T helper细胞亚型
th.sets <- c7.sets[ grep ( "GSE11924" , c7.sets $ gs_name),]
th.genes <- rowData( sce.pbmc ) $ Symbol %in% th.sets $ human_gene_symbol
summary (th.genes)
## Mode FALSE TRUE
## logical 31786 1908
# 利用GSE11961基因集来区分B细胞亚型
b.sets <- c7.sets[ grep ( "GSE11961" , c7.sets $ gs_name),]
b.genes <- rowData( sce.pbmc ) $ Symbol %in% b.sets $ human_gene_symbol
summary (b.genes)
## Mode FALSE TRUE
## logical 28185 5509 以上就是我们自己筛选的先验基因集，然后利用这些作为chosen基因再进行下游分析，看看结果与我们推测的是否一致
剔除某些先验基因
和上面相反，我们还可以先剔除一些比较不感兴趣的基因，再进行下游分析
例如，可以剔除核糖体蛋白基因或线粒体基因
复制 # 匹配核糖体蛋白基因
ribo.discard <- grepl ( "^RP[SL]\\d+" , rownames (sce.pbmc))
sum (ribo.discard)
## [1] 99
# 另一种更准确的方法
c2.sets <- msigdbr( species = "Homo sapiens" , category = "C2" )
ribo.set <- c2.sets[c2.sets $ gs_name == "KEGG_RIBOSOME" ,] $ human_gene_symbol
ribo.discard <- rownames (sce.pbmc) %in% ribo.set
sum (ribo.discard)
## [1] 87 对于免疫细胞数据集，可以剔除免疫球蛋白基因和T细胞受体基因【当然这些都取决于生物背景】
复制 library (AnnotationHub)
edb <- AnnotationHub() [[ "AH73881" ]]
anno <- select( edb, keys = rowData( sce.pbmc ) $ ID, keytype = "GENEID" ,
columns = "TXBIOTYPE" )
# 去除免疫球蛋白基因
igv.set <- anno $ GENEID[anno $ TXBIOTYPE %in% c ( "IG_V_gene" , "IG_V_pseudogene" )]
igv.discard <- rowData( sce.pbmc ) $ ID %in% igv.set
sum (igv.discard)
## [1] 326
# 去除T细胞受体基因
tcr.set <- anno $ GENEID[anno $ TXBIOTYPE %in% c ( "TR_V_gene" , "TR_V_pseudogene" )]
tcr.discard <- rowData( sce.pbmc ) $ ID %in% tcr.set
sum (tcr.discard)
## [1] 138 不过以上的操作要慎重 ，如果下游分析没有发现问题，以及对先验基因存在十足的把握，还是不要随意使用先验基因
