在某些情况下，特定基因的丰度可能会受到某些生物过程的影响，比如某些细胞类型特异性基因的表达可能会显著增加。这种强烈的上调可能导致这些基因在高丰度区间中的表现被夸大，从而影响到整个数据集的分析结果，尤其是对这些基因的变异估计。高变异基因是通过对基因表达数据进行筛选，选择那些变异程度较大的基因。在细胞类型特异性基因的强烈上调情况下，某些基因可能在丰度较高的区域集中，这样会使得这些基因的变异表现出过度增强的趋势，从而影响到拟合的准确性和其他基因的变异检测。为了解决这个问题，可以通过拟合 spike-ins 的均值依赖趋势来避免生物过程对丰度区间拟合的影响。这里的 spike-ins 是指在实验过程中添加的外源性基因，它们的表达量已知且不受生物变异的影响。因为这些外源性基因不受生物学因素的干扰，所以它们的表达趋势能够更准确地反映技术成分（如测序偏差等），而不是生物学变化。

数据集选用之前进行标准化的数据集：

sce.416b

##: 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(9): cell line cell type ... block sizeFactor## reducedDimNames(0):## mainExpName: endogenous## altExpNames(2): ERCC SIRV

在之前的介绍中我们知道，该数据集在实验设计中添加了外源性基因，因此我们直接用于计算：

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

modelGeneVarWithSpikes的作用是根据输入的单细胞基因表达数据，结合已知的 spike-ins，来估计每个基因的生物学变异和技术变异。该函数返回一个包含每个基因的估计变异信息的对象，通常包括：bio生物学变异部分，即基因在不同样本之间的变异，反映基因的生物学差异tech：技术变异部分，即测序技术的噪音或偏差。

第二段代码对 dec.spike.416b 对象进行排序，按照 bio列的值进行降序排列。order(dec.spike.416b$bio, decreasing=TRUE)将 dec.spike.416b 中的基因按照其生物变异从高到低排序，即变异最大的基因会排在最前面。[,]：这种语法用于将排序后的数据框返回。

## DataFrame with 46604 rows and 6 columns##               mean     total      tech       bio      p.value          FDR##          <numeric> <numeric> <numeric> <numeric>    <numeric>    <numeric>## Lyz2       6.61097   13.8497   1.57131   12.2784 1.51773e-186 1.57032e-183## Ccl9       6.67846   13.1869   1.50035   11.6866 2.25908e-185 2.23998e-182## Top2a      5.81024   14.1787   2.54778   11.6309  3.83585e-65  1.14102e-62## Cd200r3    4.83180   15.5613   4.22990   11.3314  9.50102e-24  6.11070e-22## Ccnb2      5.97776   13.1393   2.30179   10.8375  3.72198e-69  1.21331e-66## ...            ...       ...       ...       ...          ...          ...## Rpl5-ps2   3.60625  0.612623   6.32854  -5.71591     0.999616     0.999726## Gm11942    3.38768  0.798570   6.51471  -5.71614     0.999459     0.999726## Gm12816    2.91276  0.838670   6.57355  -5.73488     0.999422     0.999726## Gm13623    2.72844  0.708071   6.45440  -5.74633     0.999544     0.999726## Rps12l1    3.15420  0.746615   6.59326  -5.84664     0.999522     0.999726

查看一下：

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)

416B 数据集中作为平均值函数的方差。每个点代表一个基因（黑色）或加标转录本（红色），蓝线代表拟合到所有加标的趋势。

如果没有spike-in数据怎么办？

spike-ins通常用于校正技术噪声和生物噪声。在缺乏这些数据的情况下，我们可以通过对噪声进行一些合理的分布假设来构建技术噪声的趋势。例如，UMI计数通常表现出近泊松分布的变异，这意味着在没有其他因素的干扰下，技术噪声的变异呈现出一种基于计数的随机分布。我们可以基于这种假设，使用对数计数值来建立均值-方差关系。具体来说，使用对数计数可以帮助平滑数据，减少大丰度基因对拟合结果的影响，并且使数据在不同丰度区间内更加稳定。在这种均值-方差趋势构建过程中，通常会观察到在高丰度基因（即在样本中表达量很高的基因）的残差增大。残差是指实际数据点与趋势线之间的差异，反映了基因的变异度。高丰度基因的残差增大可以解释为，在使用均值-方差趋势进行拟合时，这些基因的“无趣”生物学变异。这是因为这些基因的表达量已经很高，其变异可能更多是由技术噪声或其他实验误差引起的，而不是由生物学因素引起的。因此，这些基因的残差被认为是“噪声”而非生物学上有意义的变异。

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.95605   5.05854  0.631216   4.42732         0         0## S100A9    1.93416   4.53551  0.635128   3.90038         0         0## S100A8    1.69961   4.41084  0.671512   3.73933         0         0## HLA-DRA   2.09785   3.75174  0.604479   3.14727         0         0## CD74      2.90176   3.36879  0.444953   2.92384         0         0## CST3      1.47546   2.95646  0.691408   2.26505         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)

PBMC 数据集中每个基因的标准化对数表达值的方差，与平均对数表达作对比。蓝线表示对应于泊松噪声的均值-方差关系。

基因表达数据中的噪声通常来源于测序、文库制备等技术步骤。通过拟合这种噪声趋势，可以帮助我们识别哪些基因的变异主要是由技术问题造成的，而非生物学因素。然而，纯粹基于技术噪声的趋势往往会对高表达基因（特别是一些“看家基因”）产生较大的生物学成分。这是因为这些基因在不同细胞中通常表现出较高且相对稳定的表达水平，因此，它们的表达变异可能会被误认为是技术噪声以外的“生物学”差异。这些基因通常编码细胞基本功能所需的组成成分，如核糖体蛋白等。由于它们在细胞中普遍存在并保持较高的表达水平，因此它们的表达在不同的细胞类型或状态中往往并不会发生剧烈变化，被认为是“无趣”的基因。传统上，这些看家基因被认为对表征细胞异质性并不重要，因为它们的表达较为一致，差异较小。然而实际上，这些基因的生物学变异也可能反映出细胞之间的某些差异，尤其是在不同的实验条件下。虽然更准确的噪声模型能够更好地识别技术噪声，但它可能并不一定会产生更好的高变异基因排名。这是因为某些高表达的看家基因，尽管其在不同细胞之间的表达变化较小，但其变化仍可能与细胞的某些生物学特征相关。“看家基因”也可能是差异表达基因，尽管看家基因通常被认为在细胞中没有显著变化，但它们在某些条件下可能仍然表现出差异表达。即使是这些传统上“无趣”的基因，其生物学成分的变化也可能揭示细胞之间的强烈差异。因此，不能简单地认为这些基因的变化毫无意义。

例如：

RPS16 是一个编码小亚基核糖体蛋白的基因，是典型的看家基因，通常在所有细胞中都有较为稳定的表达。在某些疾病状态下（如癌症），RPS16 的表达可能发生变化。例如，肿瘤细胞可能需要更高的蛋白合成能力，因此RPS16的表达可能相对于正常细胞发生上调。在这种情况下，尽管 RPS16 是看家基因，但它的变化反映了细胞的代谢需求变化，这可能是肿瘤细胞与正常细胞的一个差异标志。ACTB 是一个常见的细胞骨架蛋白基因，广泛用于作为内参基因。它通常在细胞中稳定表达，维持细胞形态和结构。在不同的生理或病理条件下（如肌肉细胞的生长或损伤），ACTB 的表达水平可能发生变化。例如，在肌肉损伤后，ACTB 的表达水平可能增加，以支持细胞的修复和再生过程。这种变化虽然在传统意义上被认为是“无趣”的，但它实际上可以反映出细胞的适应性反应，特别是在损伤或恢复过程中。因此，尽管这些基因常被用作内参基因，但它们的差异表达也能为我们提供关于细胞功能变化的重要信息，不应简单地忽视它们的生物学意义。