介绍

方差稳定变换（Variance Stabilizing Transformation，VST）是一种统计方法，用于将一个具有异方差性的随机变量（即方差随着均值的变化而变化的变量）转换为方差相对稳定的变量。这种转换在数据分析和建模中非常有用，因为许多统计方法（如线性回归）假设数据的方差是恒定的，即同方差性（Homoscedasticity）。通过使用VST，可以满足这些假设，从而提高模型的有效性和准确性

方差稳定变换的背景和动机

在许多实际数据集中，响应变量的方差可能会随其均值变化。例如，在Poisson分布中，方差等于均值；在二项分布中，方差与均值的函数关系为
𝑛𝑝(1−𝑝)其中 𝑝 是成功的概率。因此，当数据存在异方差性时，直接进行统计分析可能会导致错误的结论。

方差稳定变换的主要目标是找到一个合适的变换，使得变换后的数据具有稳定的方差，这样就能更好地满足统计模型的假设。

常见的方差稳定变换

平方根变换（Square Root Transformation）：
对于具有Poisson分布的计数数据，方差与均值成正比关系。通过平方根变换，可以使得方差趋于稳定。
变换公式：$Y^{\prime }$=$\sqrt{Y}$
这种变换通常用于计数数据，当数据的方差接近均值时特别有效
​
对数变换（Logarithmic Transformation）
对于数据方差随均值呈指数关系的数据，可以通过对数变换来稳定方差.
变换公式：$Y^{\prime }$=log($Y$+$c$)
其中 𝑐 是一个常数，通常用于避免对零值进行对数变换而导致的数学问题。对数变换在处理经济数据和生物数据时非常常见。
反正弦平方根变换（Arcsine Square Root Transformation）
对于比例数据或百分比数据（通常在0到1之间），可以使用反正弦平方根变换。
变换公式: $Y^{\prime }$=arcsin($\sqrt{Y}$)
这种变换在生态学和遗传学中处理比例数据时经常使用
Box-Cox变换
Box-Cox变换是一种广义的变换方法，可以通过调整参数来找到最适合数据的方差稳定变换。
变换公式：$Y^{\prime }$=$\frac{Y^{\lambda }-1}{\lambda }$($\lambda$≠0)
当 𝜆=0 时，Box-Cox变换退化为对数变换。
通过估计最佳的 𝜆 值，Box-Cox变换能够灵活适应不同类型的数据

方差稳定变换的应用场景

生物统计学
在处理生物数据（如基因表达数据、计数数据）时，经常需要应用方差稳定变换，以使得数据满足模型的假设。
例如，在RNA-seq数据分析中，通常使用对数变换来稳定基因表达的方差
经济学和金融学
对数变换在经济和金融时间序列数据中应用广泛，如股票价格、GDP等。这些数据的方差往往随着均值的增加而增加
生态学
在生态学中，处理物种丰度、种群密度等数据时，常用平方根变换或反正弦平方根变换来处理比例或计数数据
图像处理
在图像处理领域，方差稳定变换可以用于增强图像的对比度，减少噪声对图像的影响

方差稳定变换的优缺点

优点：

使数据满足模型假设：通过稳定方差，方差稳定变换帮助数据更好地满足统计模型（如线性回归）的同方差性假设。
增强数据解释性：变换后的数据通常更适合于解释和预测，特别是在非线性关系显著的情况下。
通用性：多种变换方法可以适用于不同类型的数据和分布，具有较强的适应性

缺点：

信息丢失：某些情况下，变换可能导致数据的原始信息或直观意义丢失，特别是当变换过度时。
复杂性增加：对数据进行变换后，可能需要额外的步骤将结果转换回原始尺度进行解释，这增加了分析的复杂性。
模型误差：不适当的变换可能导致模型误差的引入，特别是如果变换后的数据依然不满足模型假设时。

本文示例

我们将使用方差稳定变换技术对这些数据进行处理，以便更好地进行下游分析（如差异表达分析或聚类分析）
假设我们有RNA-seq基因表达数据，其中每个基因的表达水平是通过测序得到的原始计数数据。由于RNA-seq数据通常具有Poisson或负二项分布的特性，这意味着这些数据的方差随着均值的增加而增加。因此，直接使用这些数据进行分析可能会导致不准确的结果。为了稳定方差，我们将使用对数变换、方差稳定化变换（如VST或rlog变换）来处理这些数据

核心代码

% 清理工作区
clear;
clc;
close all;

%% Step 1: 生成或导入RNA-seq基因表达计数数据
% 假设我们有500个基因和100个样本的表达数据
numGenes = 500;
numSamples = 100;

% 生成Poisson分布的模拟数据，基因表达水平不同，方差随均值变化
rng(0); % 固定随机种子
trueMean = randi([10, 1000], numGenes, 1); % 每个基因的真实表达均值
countData = poissrnd(repmat(trueMean, 1, numSamples)); % 生成计数数据

% 可视化原始计数数据的分布
figure;
subplot(1,2,1);
histogram(log10(countData(:) + 1), 50);
title('Log10 of Raw Counts');
xlabel('Log10(Count + 1)');
ylabel('Frequency');
grid on;

% 基因表达的均值-方差关系
meanExpression = mean(countData, 2);
varExpression = var(countData, 0, 2);
subplot(1,2,2);
scatter(meanExpression, varExpression);
title('Mean-Variance Relationship');
xlabel('Mean Expression');
ylabel('Variance');
grid on;

%% Step 2: 方差稳定变换

% 对数变换
logData = log2(countData + 1);

% 方差稳定化变换（使用伽马分布逼近）
% 近似的方差稳定化变换可以考虑 rlog 转换
% 这里我们用 log(countData + pseudo-count) 的变体来近似


% 使用DESeq2等软件包实现的rlog变换通常是推荐的，这里展示简单近似
% 实际中可以考虑使用第三方包进行处理

% 可视化变换后的数据分布
figure;
subplot(1,2,1);
histogram(logData(:), 50);
title('Log2 Transformed Data');
xlabel('Log2(Count + 1)');
ylabel('Frequency');
grid on;

subplot(1,2,2);
histogram(vstData(:), 50);
title('Variance Stabilized Data');
xlabel('Transformed Expression');
ylabel('Frequency');
grid on;

% 变换后的均值-方差关系
meanVstExpression = mean(vstData, 2);
varVstExpression = var(vstData, 0, 2);
figure;
scatter(meanVstExpression, varVstExpression);
title('Mean-Variance Relationship (After VST)');
xlabel('Mean Expression');
ylabel('Variance');
grid on;

%% Step 3: 差异表达分析
% 生成两个条件下的样本（50个样本每个条件）
condition1 = vstData(:, 1:50); % 条件1的样本
condition2 = vstData(:, 51:100); % 条件2的样本

% 简单t检验来评估差异表达（实际中更常用的是DESeq2, edgeR等专用方法）


% 使用Benjamini-Hochberg方法进行p值调整
adjPValues = mafdr(pValues, 'BHFDR', true);

% 显示显著差异表达的基因数量
significantGenes = sum(adjPValues < 0.05);


% 可视化火山图（Volcano Plot）
logFoldChange = mean(condition2, 2) - mean(condition1, 2);
figure;
scatter(logFoldChange, -log10(pValues));
hold on;
scatter(logFoldChange(adjPValues < 0.05), -log10(pValues(adjPValues < 0.05)), 'r');
title('Volcano Plot');
xlabel('Log2 Fold Change');
ylabel('-Log10(p-value)');
grid on;

%% Step 4: 主成分分析和聚类分析


% 可视化前两个主成分
figure;
scatter(score(1:50, 1), score(1:50, 2), 'r', 'DisplayName', 'Condition 1');
hold on;
scatter(score(51:100, 1), score(51:100, 2), 'b', 'DisplayName', 'Condition 2');
xlabel(['PC1 (' num2str(explained(1), '%.2f') '%)']);
ylabel(['PC2 (' num2str(explained(2), '%.2f') '%)']);
legend('show');
title('PCA of Variance Stabilized Data');
grid on;

% 使用k-means聚类
k = 2; % 假设我们知道有两个群体
clusterIdx = kmeans(vstData', k);

% 可视化聚类结果
figure;
gscatter(score(:, 1), score(:, 2), clusterIdx);
xlabel(['PC1 (' num2str(explained(1), '%.2f') '%)']);
ylabel(['PC2 (' num2str(explained(2), '%.2f') '%)']);
title('K-means Clustering on PCA Results');
grid on;

代码说明

数据生成或导入：
生成了500个基因在100个样本中的RNA-seq计数数据，使用Poisson分布模拟表达水平不同的基因表达数据。

方差稳定变换：
对数变换：应用 log2(count + 1) 变换以降低方差。
方差稳定化变换：使用简单的伽马分布逼近方法来实现VST（方差稳定化变换），这是一种更复杂的对数变换。

差异表达分析：
对两个条件下的样本进行差异表达分析，使用 t-test 和Benjamini-Hochberg方法进行p值调整，并生成火山图（Volcano Plot）来可视化显著差异表达基因。

主成分分析和聚类分析：
使用PCA将数据降维，并可视化前两个主成分的分布。
使用 k-means 聚类方法对数据进行聚类分析，并将聚类结果与PCA结果结合进行可视化。

效果图示

创新扩展

rlog 变换：在实际应用中，推荐使用DESeq2等工具包实现的rlog变换，这些工具会根据数据的特性自动进行复杂的方差稳定化处理。
更复杂的统计分析：使用更复杂的差异表达分析工具，如 DESeq2、edgeR 或 limma-voom，它们专为处理RNA-seq数据而设计。
多组数据比较：可以扩展分析到多个条件下的数据比较，并使用更高级的多变量统计方法

完整代码获取

关注下列卡片公众号，回复“VST”获取完整代码