批次效应指的是实验中由于批次之间的不同处理或操作（例如样本准备、测序技术、试剂批次等）引入的系统性误差，这些误差可能导致不同批次的样本在基因表达数据上产生偏差。当多个批次的数据被合并时，这种批次效应通常会影响到分析结果，可能掩盖真实的生物学差异或产生假阳性结果。批次效应也可能使某些基因的变异看起来很大，因此，这些基因可能会被误认为是具有生物学意义的高变异基因，而实际上它们的变异仅仅是由于批次效应所导致的。因此，我们通常不关心由批次效应驱动的 HVG，而是希望识别每个批次内的高变异基因，避免批次效应掩盖了真正的生物学差异。
为了有效地分析每个批次中的高变异基因，我们需要为每个批次分别执行趋势拟合和方差分解。这意味着我们会单独处理每个批次的数据，拟合它们的均值-方差趋势，识别批次内的高变异基因。在实践中，我们可以对每个批次进行单独的方差建模，这样在每个批次内我们都能够识别出高度变异的基因，而不会受到其他批次数据的影响。
dec.block.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC", block=sce.416b$block)head(dec.block.416b[order(dec.block.416b$bio, decreasing=TRUE),1:6])
modelGeneVarWithSpikes() 是用于基因变异分析的函数，通常用于估计基因的技术噪声和生物学变异，并将外源spike-in 基因（如ERCC基因）用于校正技术噪声。这个函数不仅会分析每个基因的总变异，还会使用 spike-in 基因来估计和校正实验中的技术噪声成分。 sce.416b$block 指定批次信息，使用block参数即对不同批次进行了初步处理。
## 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.58417 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.74575 11.2734 3.96145e-137 8.57006e-135## Cd200r3 4.83305 15.5909 4.31900 11.2719 1.18818e-54 7.06875e-53## Ccnb2 5.97999 13.0256 2.46649 10.5591 1.22233e-151 3.02999e-149## Hbb-bt 4.91683 14.6539 4.12163 10.5322 2.54730e-49 1.35308e-47
特定于批次的趋势拟合有助于更好地控制批次效应，尤其是在批次之间存在技术差异时。在每个批次中，趋势拟合可以分解出生物学成分和技术成分，这些成分有助于识别真实的基因表达变化，并消除技术噪声的影响。最终，通过将来自不同批次的信息进行整合，可以得到更准确的基因表达变异估计，并减少单个批次的技术误差对结果的影响。
可视化一下：
par(mfrow=c(1,2))blocked.stats <- dec.block.416b$per.blockfor (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) }
左边的图代表批次 20160113，右边是批次 20160325。你可以看到，在两个批次中，均值-方差关系的形状非常相似，这表明尽管可能存在一些技术差异，实验在两个批次之间保持了一定的重复性。如果在可视化中可见明显的批次效应，那么则需要用其他的方法。
注意：
同一数据集的批次效应主要来自实验中的技术差异，通常可以通过block 参数等方法控制。不同数据集的批次效应则可能由于实验平台、样本来源等因素的不同而更加复杂，需要使用更加复杂的批次效应修正方法（如 ComBat）来处理。
HVG如何选择？
如果选择较大的HVG子集，意味着保留了更多的基因。这有一个优点，就是它能减少“丢失”潜在有意义信号的风险，因为可能有更多与生物学过程相关的基因被保留。然而，更多的基因意味着数据中的噪声也会增加，尤其是那些与研究目标不相关的基因。这些噪声可能掩盖有用的生物学信号，从而影响下游分析的准确性。相对而言，选择较小的子集可以减少噪声，但也可能丢失一些重要的相关基因。因此，选择的基因子集大小需要在这些因素之间进行权衡。
在不同的研究中，噪声和信号的界限是模糊的。有些情况下，某些看似“噪声”的基因实际上可能包含有用的信息，特别是在揭示生物学异质性的研究中。例如，在T 细胞活化反应的异质性研究中，可能需要保留更多的基因来捕捉反应的多样性。但如果研究的目标仅仅是区分主要的免疫表型（例如不同类型的免疫细胞），那么这种异质性可能会被认为是“噪声”，应该被忽略。
选择策略是什么？
一种常见的选择方法是选择那些具有最大方差的基因。这些基因通常是在不同样本或条件下表现出显著变化的基因，往往是生物学上相关的基因。对于生物学变异的选择，通常选择那生物成分最大化的基因，即那些具有最大的生物学变异的基因。这些基因更有可能反映了细胞之间的真实生物学差异。通过选择前 n 个基因（通常是那些具有最大方差的基因），可以精确控制后续分析中保留的基因数量，从而使得计算复杂度可控。这对于处理大规模数据集时尤为重要。
hvg.pbmc.var <- getTopHVGs(dec.pbmc, n=1000)str(hvg.pbmc.var)
## chr [1:1000] "LYZ" "S100A9" "S100A8" "HLA-DRA" "CD74" "CST3" "TYROBP" ...
代码片段：可切换语言，无法单独设置文字格式
选择合适的 n 需要根据数据的异质性来决定。对于异质性较低的数据集，可以选择较小的 n，以减少噪声；对于异质性较高的数据集，可以选择较大的 n，以捕获更多生物学信号。选择 n 的一个缺点是它将基因选择过程变成了基因间的竞争，可能会导致重要的标记物被排除。特别是在高度异质性的数据集里，选择不当可能影响对不同亚群体的分辨能力。对于大多数分析，选择一个合理的 n 后进行分析，必要时再进行调整，而不是过于纠结于“最佳”值。
一种选择思路：
dec.pbmc <- modelGeneVar(sce.pbmc)chosen <- getTopHVGs(dec.pbmc, prop=0.1)str(chosen)
选择百分之十
我们可以进行子集化，以仅保留我们选择的 HVG。 这确保了下游方法将仅使用这些基因进行计算。 缺点是非 HVGs被丢弃，这使得查询非 HVG 的有趣基因的完整数据集稍微不方便。
sce.pbmc.hvg <- sce.pbmc[chosen,]dim(sce.pbmc.hvg)
## [1] 1274 3985
我们可以保留原始对象，并通过一个额外的参数（ 如果分析在不同步骤中使用多组 HVG，则这非常有用，因此在特定步骤中可以很容易地将一组 HVG 交换为另一组 HVG。
rowSubset(sce.pbmc, "HVGs.more") <- getTopHVGs(dec.pbmc, prop=0.2)rowSubset(sce.pbmc, "HVGs.less") <- getTopHVGs(dec.pbmc, prop=0.3)colnames(rowData(sce.pbmc))
## [1] "ID" "Symbol" "subset" "HVGs.more" "HVGs.less"
