基于概率白化变换的超高维遗传数据快速处理方法

《PLOS Genetics》：Fast probabilistic whitening transformation for ultra-high dimensional genetic data

　　

： 统计方法通常对数据集的样本或特征之间的独立性作出假设，然而相关性结构在真实数据中普遍存在。在统计遗传学中，由邻近变异之间的连锁不平衡（linkage disequilibrium, LD）导致的回归汇总统计量之间的相关性使下游分析复杂化。许多方法使用LD参考面板来考虑这种相关性。白化变换（whitening transformation）是一类去除相关性结构的方法，被许多常见模型隐式使用。本文提出了一种用于数据白化的概率模型，并提出了参数估计和高效算法以提升统计与计算性能。在将精神分裂症全基因组关联研究（genome-wide association study, GWAS）中未观测遗传变异的z统计量进行插补的应用中，概率白化变换在均方误差最低的同时，比其他方法快一个数量级。利用该方法，研究人员还识别出能够解释疾病相关基因遗传调控信号的串联重复序列（tandem repeats）。分析已实现于研究人员新开发的开源R包decorrelate和imputez中。

论文解读文章

一、研究背景与问题

在真实世界数据集中，相关性结构普遍存在，然而许多统计模型往往假设数据集的样本或特征之间相互独立，这一假设在实践中常常无法满足。白化变换是一类广泛使用的预处理方法，旨在去除数据中的相关性结构后再进行其他分析。对于包含n个样本和p个特征的数据矩阵，白化变换通过计算特征间的样本协方差矩阵，并对观测数据施加一个非唯一的白化矩阵，使变换后数据的协方差矩阵近似为单位矩阵。白化变换既可被显式地用作预处理步骤，也更多地被隐式地嵌入到多元回归、典型相关分析、线性判别分析、高斯混合模型、马氏距离、多元离群点检测以及线性混合模型等广泛使用的方法中。

统计遗传学领域大量使用隐式白化变换来处理邻近遗传变异之间的相关性，即连锁不平衡（LD）。LD是GWAS回归汇总统计量（如系数估计或z统计量）下游分析中的主要复杂因素。实践中通常从外部参考面板估计相关矩阵，这些样本不包含在目标研究中，随后这些样本外估计被用作预处理步骤或统计模型内部组件。该方法已成为表型预测、统计精细定位、z统计量插补、异常z统计量检测、多重检验校正、风险变异基因组注释富集分析以及模拟具有相似相关结构的敲除数据集等应用的标准做法。

然而，标准白化变换存在两方面的显著局限。计算上，评估协方差矩阵并构造白化矩阵需要p2量级的内存和p3量级的时间，对于大规模p值而言变得不可行。统计上，在高维数据（n < p）情形下，样本协方差矩阵不再满秩或不可逆，变换必须加以调整。直接估计协方差矩阵和精度矩阵是一个庞大的研究领域，许多方法通过正则化确保估计矩阵具有特定性质，如收缩样本特征值以保证正定性、通过收缩偏相关系数对应的回归系数保证稀疏性、正则化条件数等。但白化变换作为预处理步骤或统计模型内部组件时，并不需要直接估计协方差矩阵，协方差矩阵仅仅是冗余参数（nuisance parameters）。这一重要区别促使研究人员提出隐式协方差方法，在不直接估计协方差或白化矩阵的情况下执行白化变换。

二、研究内容与主要结论

研究人员提出了高斯逆威沙特经验贝叶斯（Gaussian Inverse Wishart empirical Bayes, GIW-EB）概率白化变换模型。该模型将观测数据建模为矩阵正态分布，并对列间协方差矩阵施加逆威沙特（inverse Wishart, IW）先验分布。后验期望协方差是样本协方差与先验目标矩阵的凸组合，其中目标矩阵设为缩放的单位矩阵，缩放因子为所有特征的平均方差。研究人员通过经验贝叶斯方法边际化协方差矩阵并最大化观测数据的似然函数来估计收缩参数，数据仅通过样本奇异值进入对数似然。

该概率白化变换具有重要的统计性质。后验期望协方差与样本协方差共享相同的特征向量，但特征值向平均方差方向收缩，因此始终正定且可逆，即使样本协方差矩阵不可逆时亦然。白化变换保留原始奇异向量，仅修改奇异值，这一性质带来了显著的计算优势。研究人员推导出关键恒等式，使得白化变换仅需非零奇异值及其对应奇异向量即可精确计算，从而将计算复杂度从特征数量的三次方降低到线性。在超高维情形下，通过秩k部分奇异值分解（singular value decomposition, SVD），计算复杂度进一步在样本数和特征数上均达到线性。

在模拟数据研究中，研究人员比较了11种白化变换的样本外性能。GIW-EB方法在四种相关性场景中均一致地给出最低的样本外均方根误差（root mean squared error, rMSE），尤其在协方差矩阵为低秩、特征谱由前几个特征值主导时改进最为显著。在计算性能方面，n = 5000样本且特征数不断增加时，GIW-EB方法运行时间小于2分钟，低秩GIW-EB（k = 50）小于1分钟，而标准方法在时间和内存上均变得不可行。

在1000 Genomes Project遗传数据上，研究人员分别对欧洲、亚洲和非洲血统个体进行了独立LD区块的白化变换评估。尽管遗传数据并非连续或多元高斯分布，收缩方法在欧洲个体的大多数窗口中表现良好。GIW-EB与从数据估计收缩参数的其他方法相当或更具竞争力，低秩GIW-EB（k = 50）在少于1500个变异的窗口中最具竞争力。使用隐式协方差算法的方法计算性能最佳，墙钟时间约5分钟，而其他方法需要3至30小时。

三、主要技术方法

研究人员主要采用了以下关键技术方法：一是概率模型构建，将观测数据建模为矩阵正态分布，并对协方差矩阵施加逆威沙特先验，通过后验期望得到正则化协方差估计；二是经验贝叶斯收缩参数估计，边际化协方差矩阵后最大化边际似然，数据仅通过样本奇异值进入似然函数；三是基于奇异值分解的隐式协方差算法，利用白化变换仅修改奇异值而保留奇异向量的性质，避免直接构造和求逆协方差矩阵；四是低秩近似策略，通过秩k部分SVD将计算复杂度降至样本数和特征数的线性；五是将白化变换应用于z统计量插补，基于条件多元高斯分布推导插补公式。在应用研究中，使用了1000 Genomes Project参考面板（欧洲血统407个无关个体）以及精神分裂症GWAS汇总统计量（58K病例和77K对照）。

四、研究结果

性能模拟研究：在四种相关性场景（恒定相关、自相关、分块恒定相关、分块自相关）中，GIW-EB方法在所有p值范围内均给出最低样本外均方根误差。无收缩方法（λ = 0）即使在小p值下表现不佳，伪逆方法在p > n时有所改进但无法与收缩方法竞争。GIW-EB方法误差随p增大而减小，避免了其他收缩方法在p接近n时误差增大的问题。

1000 Genomes Project遗传数据性能：对欧洲、亚洲和非洲血统个体的独立LD区块进行样本外白化评估，GIW-EB与数据驱动估计收缩参数的其他方法相当或更具竞争力。低秩GIW-EB（k = 50）在少于1500个变异的窗口中最具竞争力。隐式协方差算法方法计算性能最佳，墙钟时间约5分钟，其他方法需3至30小时。

GWAS汇总统计量插补应用：使用精神分裂症GWAS（58K病例和77K对照，欧洲血统）的560万变异汇总统计量，基于1000 Genomes Project 407个欧洲血统无关个体参考面板，采用1 Mb滑动窗口和250 kb侧翼区域插补1500万变异的z统计量。GIW-EB在广泛的z统计量值范围内给出高插补精度，全基因组均方根误差最小。插补误差随次要等位基因频率（minor allele frequency, MAF）减小而增大，GIW-EB在MAF范围内给出最小平方误差。计算上，GIW-EB方法全基因组分析所需时间小于1小时，其他方法需4至126小时。

串联重复序列eQTL信号插补应用：研究人员整合了包含1000 Genomes Project供体的170万串联重复序列群体参考面板与人类大脑基因表达遗传调控研究（近4000个bulk RNA-seq样本）的汇总统计量，对3个疾病相关基因进行原理验证。在C9orf72（肌萎缩侧索硬化关键风险基因）、TSPAN14和TYK2（阿尔茨海默病风险相关基因）中，研究人员识别出比区域中SNP与基因表达关联更强的串联重复序列，证明z统计量插补能够识别解释疾病相关基因遗传调控信号的串联重复序列。

五、讨论与结论

白化变换是许多广泛使用的统计方法中核心但常被忽视的组成部分。标准白化变换受限于计算时间、内存使用，并在中等规模数据集上给出较差的样本外性能。正则化协方差估计改善了统计性能，但对高维数据集仍计算需求极高。研究人员提出的正则化概率白化变换基于数据矩阵的奇异值分解考察其统计性质，概率模型的特殊形式使得隐式协方差方法能够避免直接估计协方差矩阵，将计算时间从特征数量的三次方降至线性，实际应用中可减少约一个数量级的计算时间。在模拟和1000 Genomes Project高维基因型数据应用中，GIW-EB方法在样本外白化基准测试中匹配或超越其他方法。在精神分裂症GWAS z统计量插补应用中，GIW-EB给出最佳插补性能。作为原理验证，研究人员证明基于GIW-EB白化变换的z统计量插补能够识别解释肌萎缩侧索硬化和阿尔茨海默病相关基因遗传调控信号的串联重复序列。随着基因组数据集规模持续增大，概率模型的隐式协方差方法将使白化变换能够应用于现有方法不可行的高维数据，未来应用包括单细胞基因表达和更大的遗传参考面板。