# 3.3 挑选表达量高变化基因

### 1 前言

我们经常使用scRNA数据探索细胞异质性。聚类分群和降维也是常用的操作，它们都是依赖于基因表达量（例如综合细胞中各个基因的表达量，然后把表达模式相近的细胞聚在一起）。因此**挑选哪些基因进行聚类分群和降维分析，这是非常关键的**，挑选的基因一定要有代表性，尽可能多的包含生物信息而剔除其他技术噪音。

单细胞分析的一个重点就是：用尽可能少的维度来展示数据真实的结构。降维的过程其实就是去繁存简，每个基因的变化都对整体有影响，对细胞来讲都是一个变化维度。但这么多维度我们看不了也处理不了，需要尽可能保证真实差异的前提下减少维度的数量，因此需要挑出那些更能代表整体差异的基因。

对每一个细胞来讲，2万多个基因都是它身上长向四面八方的维度，但其中有一些影响很大，一些影响很弱。我们试图把影响最大的几个维度找出来（比如Gene1、Gene2、Gene3），然后用这几个维度代表这个细胞的生物学特征。

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-23-012726.png)

**上面说了这么多，重点还是：如何挑选有价值的基因（feature）**

一般来说，如果一个基因在细胞群体中变化幅度很大，就是受关注对象，我们也会认为是生物因素导致了这么大的差异。这样的基因叫做：**高度变化基因（highly variable genes ，HVGs）**

### **数据准备**

**10X PBMC**

其中的处理步骤值得学习，其中`pbmc4k_raw_gene_bc_matrices.tar.gz`数据已经帮大家下载好，链接：<https://share.weiyun.com/s2wbVd7p> 密码：3iwypw

下载的内容包含三件套，然后使用`read10xCounts`读取：

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-005621.png)

```r
# 下载
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个基因在线粒体上

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-010242.png)

```r
# 接下来
# 空液滴检查
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数据**

```r
# 加载数据
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的效果一样

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-013155.png)

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

### **2.1 方法一：使用log-counts衡量变化程度**

> 这也是最简单的处理方法，基因log后的表达量在细胞间变化幅度越大，就说明细胞间距离越远（这个距离指欧氏距离Euclidean distances）

计算每个基因的变化程度很简单，但要基于这么多基因的结果去挑选最有价值的基因（feature），还需要构建一个模型。这个模型叫：mean-variance relationship，看上去与均值、方差有关

**先看计算方法**

```r
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)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-020520.png)

* 每个点表示一个基因
* 图中蓝线指的是：技术因素导致的偏差
* 纵坐标表示总偏差：它等于技术偏差+生物因素偏差

因此，要衡量一个基因的生物因素偏差大小，就看对应的纵坐标减去对应的蓝线的值

最后可以得到一个表：

```r
dec.pbmc[order(dec.pbmc$bio, decreasing=TRUE),]
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-020940.png)

可以看到生物因素（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()`来计算

```r
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.**

```r
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)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-021748.png)

**最后也能得到一个表：**

其中偏离拟合曲线的数值用`ratio`表示，也是我们定义HVGs的途径

这个`ratio = total (CV^2) / trend (CV^2)` ，这样一相除，就抵消了均值的影响，最后实际上变成了方差的比值。但如果是相减，还是会带着均值【把公式从纸上一列就很清楚了】

```r
dec.cv2.pbmc[order(dec.cv2.pbmc$ratio, decreasing=TRUE),]
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-022009.png)

**总结：**

* log-counts与CV2都是衡量基因变化程度的有效途径
* 优先使用log-counts的结果，因为下游的分析也是基于log-counts值

### **2.3 考虑技术噪音**

> 有两种方式：一种是有spike-in的，一种是没有的，都可以操作

之前画的技术偏差拟合线都是基于一个假设：大部分基因的表达量都是受外界技术误差的影响。但实际上，除了技术误差，还有一些生物因素会导致数据波动（例如transcriptional bursting）。因此之前估计的技术偏差可能”超出了实际的大小“，有些过度估计。

之前估计的偏差实际上是：真正的技术偏差 + 一些无关紧要的生物学因素偏差。因此如果有spike-in的话，它和内源基因不同，不会受到生物因素的影响。利用spike-in画的拟合线，更能代表技术偏差

**如果有spike-in**

```r
dec.spike.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC")
dec.spike.416b[order(dec.spike.416b$bio, decreasing=TRUE),]
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-032639.png)

```r
# 作图
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)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-032709.png)

其中黑点是基因，红点是spike-in，蓝线是根据红点画的拟合线

**如果没有spike-in**

可以考虑利用数据分布来表示技术噪音，例如只考虑技术噪音的话，UMI counts通常会呈现近似泊松分布

```r
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)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-032808.png)

有趣的是，如果真的单纯考虑技术偏差，我们会得到比原来更多的变化基因。其中就包含了之前没有的一些管家基因。不过这些基因在后面分析（比如研究细胞异质性）中又不太感兴趣，因此**使用更准确的技术噪音拟合方法，可能不会得到更准确的HVGs排序**

### **2.4 考虑批次效应**

**2.4.1 方法一：绘制拟合曲线时加入批次信息（Fitting block-specific trends）**

有时我们看到的在不同细胞间高变化的基因，不一定是生物学因素驱动，也不一定是由于技术误差导致，而是由于不同细胞的批次不同，导致同一基因的表达量差异很大。

我们更想知道的是，在一个批次中的HVGs信息。一般可以一个批次一个批次地去处理，例如：

```r
# 这里用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
```

看看不同批次的差异

```r
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) 
}
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-040016.png)

其中黑点是基因，红点是spike-in，蓝线是根据红点画的拟合线。看到这里的两个批次之间差异很小，表示重复效果不错

**2.4.2 利用实验设计矩阵（design matrix）**

当只有一个批次信息时（例如只考虑`sce.416b$block`信息），使用上面的绘制拟合曲线方法是足够的，但如果再加上其他的批次信息，例如：

```r
> table(colData(sce.416b)$phenotype)

induced CBFB-MYH11 oncogene expression 
                                    96 
                   wild type phenotype 
                                    96
```

就需要利用实验设计矩阵（很像DESeq2、edgeR等差异分析做的），像是`modelGeneVarWithSpikes()`和`modelGeneVar()` 都支持这个参数

```r
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),]
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-28-113935.png)

这种方法做起来很简单，但不太准确，因为没有考虑到各个批次中的平均表达水平。而这个表达水平又影响真正的技术偏差估计

> 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`这一列提取

```r
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`这一列提取

```r
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的阈值，例如

```r
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`那一列进行排序，从而挑选基因呢？

最好不要，因为还是那句话，**表现显著的基因不一定对我们感兴趣的生物学问题有帮助。**&#x5F88;多基因确实p值很小，这是因为它们的技术偏差非常小，但同时总体偏差并不大，因此它看上去很显著，但实际我们不感兴趣。

因此这种方法还是慎用

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

这个方法的目的很单纯，拟合曲线以下的基因都是不感兴趣的，因为它们的生物学偏差很小，于是把它们全部剔除，留下剩下的。这也是走了一个极端：为了偏差最小化，选择了噪音最大化（保留了很多与研究问题无关的基因）。看下面计算的数量就知道

如果之前使用了`modelGeneVar()`：

```r
hvg.pbmc.var.3 <- getTopHVGs(dec.pbmc, var.threshold=0)
length(hvg.pbmc.var.3)
## [1] 12791
```

如果之前使用了`modelGeneCV2()`：保留所有ratio大于1的

```r
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`

```r
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的基因也有点意思，想看看它们就变得很麻烦

```r
sce.pbmc.hvg <- sce.pbmc[chosen,]
dim(sce.pbmc.hvg)
## [1] 1279 3922
```

### **方法二：原数据不动，只是下游分析基于小部分HVGs**

这样的好处是，如果以后调整了chosen这个HVGs基因集，那么还是可以继续进行下游分析的

```r
# 例如只根据HVGs进行PCA，使用参数subset_row
library(scater)
sce.pbmc <- runPCA(sce.pbmc, subset_row=chosen)
reducedDimNames(sce.pbmc)
## [1] "PCA"
```

### **方法三：原数据先备份，下游分析也是基于原数据，最后再复原**

```r
# 首先保留一份原数据
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

```r
# 首先选取基因集
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`基因再进行下游分析，看看结果与我们推测的是否一致

### **剔除某些先验基因**

和上面相反，我们还可以先剔除一些比较不感兴趣的基因，再进行下游分析

例如，可以剔除核糖体蛋白基因或线粒体基因

```r
# 匹配核糖体蛋白基因
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细胞受体基因【当然这些都取决于生物背景】

```r
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
```

**不过以上的操作要慎重**，如果下游分析没有发现问题，以及对先验基因存在十足的把握，还是不要随意使用先验基因


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://jieandze1314.osca.top/03/03-3.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.