从所有基因中挑选出那些对细胞间异质性有代表性的基因，而排除掉那些只是由于技术噪音或无关的生物变异引起的变化的基因即寻找高变基因。在进行细胞间比较时（例如，进行聚类和降维），基因的选择会直接影响细胞间相似性度量的结果。这些度量通常会基于基因表达的差异来计算，而选择的基因会影响到计算出来的“相似性”或“不相似性”值。为了尽量保留有生物学意义的结构而去除无关的噪声，研究者通常会选择变异性较大的基因。这是因为，真正的生物学差异（例如，细胞类型、状态等差异）通常会导致某些基因在不同细胞之间表现出较大的变异，而这些基因更有可能提供关于生物学现象的信息。相反，其他基因的变异可能只是技术噪声或基线水平的波动，这些基因对生物学解释帮助不大。在选择变异性较大的基因时，通常会使用几种方法来量化每个基因的变异度，并从中挑选出那些表现出最大变异的基因（这些基因通常被称为“高度可变基因”）。这些基因将用于后续分析，以提高计算效率并增强分析结果的生物学相关性。
在这里，我们描述了对 10X Genomics 的外周血单核细胞 （PBMC） 数据集的简要分析。 这些数据可从 10X Genomics 网站公开获得。 我们从中下载原始基因/条形码计数矩阵，并读取：
library(DropletTestFiles)raw.path <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/raw.tar.gz")out.path <- file.path(tempdir(), "pbmc4k")untar(raw.path, exdir=out.path)library(DropletUtils)fname <- file.path(out.path, "raw_gene_bc_matrices/GRCh38")sce.pbmc <- read10xCounts(fname, col.names=TRUE)
getTestFile()从DropletTestFiles包中获取一个包含PBMC数据的压缩文件raw.tar.gz的路径。tempdir()返回一个临时目录路径，file.path()将此路径与指定的目录名（"pbmc4k"）组合，生成一个输出路径。untar()将压缩包解压到指定的out.path目录下。生成包含基因表达数据的文件路径fname，该路径指向解压后PBMC数据集的GRCh38文件夹。read10xCounts()函数将读取10x Genomics的基因表达矩阵，并将其加载到一个SingleCellExperiment对象（sce.pbmc）中。参数col.names=TRUE表示文件中包含列名称。
查看一下：
sce.pbmc
##: SingleCellExperiment ## dim: 33694 3985 ## metadata(1): Samples## assays(2): counts logcounts## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B## rowData names(2): ID Symbol## colnames(3985): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...## TTTGTCAGTTAAGACA-1 TTTGTCATCCCAAGAT-1## colData names(3): Sample Barcode sizeFactor## reducedDimNames(0):## mainExpName: NULL## altExpNames(0):
最简单的量化每个基因变异的方法是计算每个基因在所有细胞中的对数标准化表达值（对数计数）的方差。这意味着，对于每个基因，我们计算的是其在所有细胞中对数表达值的变化程度，而不是原始的计数数据。计算每个基因的方差可以衡量基因在不同细胞之间的变化程度。对于变异性较大的基因，表示该基因在不同细胞之间的表达差异较大，可能与细胞的状态或类型有关，而变异性较小的基因则可能反映出稳定的、无关的背景噪声。
尽管对数计数的方差计算方法简单，但在进行特征选择时，我们需要建模基因表达的均值-方差关系。
均值-方差关系：在基因表达数据中，通常基因的变异性（方差）与其丰度（均值）是相关的。通常，表达量高的基因会显示出更大的变异性，因为它们的信号更强，容易受噪声或技术波动的影响。这种现象导致了基因表达中的“均值-方差”关系。
对数变换的局限性：尽管取对数可以帮助处理基因表达的偏态分布，但对数变换本身并不是一种方差稳定变换。也就是说，基因的变异更多受到基因丰度的影响，而不是潜在的生物学异质性。这意味着，丰度较高的基因通常会有较高的方差，而这可能掩盖了真正的生物学差异。为了更好地理解和解释这种“均值-方差”关系，通常会对所有基因的方差进行建模，并拟合出一个趋势线。这一趋势线帮助我们理解不同丰度基因的变异性是如何变化的，并进一步指导我们在特征选择时识别出哪些基因可能是由于真正的生物学差异而具有较高的变异性，而不仅仅是因为其丰度较高。
library(scran)dec.pbmc <- modelGeneVar(sce.pbmc)
计算流程
总方差：对于每个基因，首先计算其在所有细胞中的方差，即总方差。这反映了基因表达的总体变异性。由于技术噪声（如采样误差、测序误差等）会影响基因的变异，modelGeneVar会拟合一个均值-方差趋势，用来估计基因表达量与方差之间的关系。这条拟合的趋势代表了基因的“技术成分”，即由技术因素引起的方差。生物学成分是通过从总方差中减去技术成分来获得的。这部分变异代表的是基因表达中的“有趣变异”，即由真正的生物学差异（如细胞类型、状态或环境变化等）引起的变异。
modelGeneVar(sce.pbmc)返回一个包含每个基因的以下信息的对象：
●var：基因的总方差。
●mean：基因的均值。
●bio：生物学成分，表示去除了技术噪声后的基因表达变异。
●tech：技术成分，表示由技术因素引起的变异。
可视化：
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)
fit.pbmc <- metadata(dec.pbmc)从dec.pbmc对象中提取元数据（metadata），该元数据包含了每个基因的表达均值、方差以及方差与均值的关系趋势。plot(fit.pbmc$mean, fit.pbmc$var, xlab="Mean of log-expression", ylab="Variance of log-expression")使用plot()函数绘制基因的均值（fit.pbmc$mean）与方差（fit.pbmc$var）之间的散点图。横坐标是基因表达的对数均值（log-expression），纵坐标是基因表达的方差（log-expression的方差）。curve(fit.pbmc$trend(x), col="dodgerblue", add=TRUE, lwd=2)使用curve()函数在现有的散点图上添加一条拟合曲线，表示基因的均值与方差之间的关系。fit.pbmc$trend(x)给出了拟合的方差趋势。add=TRUE表示将这条拟合曲线添加到已经存在的散点图中。
了解了基本的计算方法之后，我们还要考虑如何将基因表达拆分为生物因素和技术因素：
技术分量
在任何给定的基因丰度（即基因表达水平）下，大多数基因的表达变化是由采样噪声、技术变异等无关的过程驱动的。这些噪声通常不是生物学意义上的差异，而是与测序过程、样本处理或其他实验技术相关的。对于每个基因，根据其丰度和方差之间的关系，拟合出的趋势线代表了基因的“无趣变异”或“技术分量”。这一部分的方差反映的是由于技术噪声（而非生物学差异）引起的变化。例如，基因丰度较高时，其变异往往更多地由技术噪声而非真实的生物学差异所影响。
生物学分量
基因的总方差可以被分解为两个部分：一个是技术成分（技术噪声），另一个是生物学成分。生物学成分代表的是基因的“有趣”变异，即由真正的生物学差异（如细胞类型、状态等）引起的变异。生物学成分通过从基因的总方差中减去技术成分（由趋势拟合得出）来计算。具体来说，生物学成分 = 总方差 - 技术成分。这部分方差反映的是基因表达的真正生物学差异。
高度可变基因（HVG）选择
选择高度可变基因的关键在于识别那些表现出较大生物学成分（有趣变异）的基因，而不是仅仅依赖于技术噪声的变异。这些基因可能在不同细胞类型、状态或环境条件下表现出显著的变化，因此具有较高的生物学意义。
选择生物学变异：
dec.pbmc[order(dec.pbmc$bio, decreasing=TRUE),]
order(dec.pbmc$bio,decreasing=TRUE)中order()函数根据dec.pbmc$bio的值对基因进行排序。decreasing=TRUE表示按降序排列，即将生物学成分较大的基因排在前面。这样，我们就可以将具有最大生物学变异性的基因排在前面。dec.pbmc[order(dec.pbmc$bio,decreasing=TRUE), ]通过order()函数得到排序后的基因索引，使用该索引重新排列dec.pbmc对象中的基因。结果是一个按生物学变异性降序排列的基因数据框。
## DataFrame with 33694 rows and 6 columns## mean total tech bio p.value FDR## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>## LYZ 1.95605 5.05854 0.835339 4.22320 1.06165e-270 2.08816e-266## S100A9 1.93416 4.53551 0.835435 3.70007 2.62793e-208 7.38410e-205## S100A8 1.69961 4.41084 0.824339 3.58650 4.19449e-201 9.16683e-198## HLA-DRA 2.09785 3.75174 0.831235 2.92051 5.82146e-132 4.77093e-129## CD74 2.90176 3.36879 0.793182 2.57561 4.74972e-113 2.45848e-110## ... ... ... ... ... ... ...## TMSB4X 6.08142 0.441718 0.679233 -0.237515 0.992451 1## PTMA 3.82978 0.486454 0.731271 -0.244817 0.990003 1## HLA-B 4.50032 0.486130 0.739577 -0.253447 0.991377 1## EIF1 3.23488 0.482869 0.768938 -0.286070 0.995136 1## B2M 5.95196 0.314948 0.654261 -0.339313 0.999843 1
在进行均值-方差趋势拟合时，可能会出现一些基因的生物学成分为负值。这些基因的生物学成分为负值的情况在统计上是可能发生的，但通常没有明确的生物学解释。这些负值可能是由于数据中存在的噪声、拟合模型的局限性或一些极端的观测值。尽管存在负值，但它们通常不会对后续分析产生重大影响，特别是在大多数应用中，负的生物学成分可以被忽略或视为技术噪声的表现。在拟合每个基因的差异时，一些基因的表达差异低于拟合的趋势线，这可能导致负的生物学成分。这是因为趋势拟合的过程中，总体上基因的表达量可能受到技术噪声的影响，而在某些情况下，基因的变异性表现得较低。理论上，modelGeneVar函数假设大多数基因的表达谱变异主要由随机技术噪声主导。因此，拟合的趋势通常被解释为“技术成分”，即这种变异是由测序技术、样本处理等因素引起的。然而，实际上，大部分表达的基因都会展示一些生物学变异性，甚至在一些噪声较大的数据中也会出现这种变异。例如，转录爆发是基因表达的一个常见现象，其中即使是“无趣”的基因也会因为转录事件而表现出一定的变异。因此，假设大多数基因的变异主要由技术噪声加上“无趣的生物变异”（即不影响重要生物学异质性的变异）来解释可能更加合适。也就是说，尽管技术噪声占主导，但基因表达的变异性也受到一些无关的生物学因素的影响，比如转录爆发、偶发的转录活动等。将拟合趋势解释为技术噪声加上无关的生物变异比将其完全视为纯技术噪声更加符合实际数据的特点。这种方法可以帮助我们更好地理解基因的变异性，尤其是在高度可变基因选择和生物学分析中。
什么是转录爆发？
转录爆发 是指在细胞中，基因的转录活动并不是平稳持续的，而是以短时间内的高频次活动（即“爆发”）和长时间的低频或静默期交替进行。这个现象意味着，在某些时间点，基因的表达水平会突然增加，然后迅速回落，进入低表达或静默状态。
