优雅使用SciPy进行分位数归一化：基因表达数据分析实战

引言

在基因表达数据分析中，数据标准化是一个至关重要的预处理步骤。本文将介绍如何使用NumPy和SciPy实现分位数归一化(Quantile Normalization)，这是一种强大的标准化技术，能够确保不同样本的测量值具有相同的统计分布。

分位数归一化原理

分位数归一化是一种使不同样本具有相同分布的技术，主要包含三个步骤：

1. 对每列数据进行排序
2. 计算排序后每行的平均值
3. 用平均列的对应分位数替换每列的分位数

这种方法的优势在于：

• 消除技术变异带来的分布差异
• 保留样本间的相对排序关系
• 适用于各种高通量数据，如微阵列和RNA-seq

实现分位数归一化

下面是用NumPy和SciPy实现分位数归一化的完整代码：

import numpy as np
from scipy import stats

def quantile_norm(X):
    """将X的每列归一化为相同分布
    
    参数:
    X : 2D数组, 形状(M, N)
        输入数据，M行(基因/特征)和N列(样本)
    
    返回:
    Xn : 2D数组, 形状(M, N)
        归一化后的数据
    """
    # 计算分位数
    quantiles = np.mean(np.sort(X, axis=0), axis=1)
    
    # 计算每列的秩
    ranks = np.apply_along_axis(stats.rankdata, 0, X)
    
    # 将秩转换为0到M-1的整数索引
    rank_indices = ranks.astype(int) - 1
    
    # 使用秩矩阵索引每个秩的分位数
    Xn = quantiles[rank_indices]
    
    return Xn

def quantile_norm_log(X):
    """对数转换后的分位数归一化"""
    logX = np.log(X + 1)  # 加1避免log(0)
    logXn = quantile_norm(logX)
    return logXn

关键技术点解析

这段代码展示了NumPy的多个强大特性：

1. 多维数组操作：处理2D矩阵形式的数据
2. 向量化运算：使用单行代码完成整个数组的对数变换
3. 轴操作：通过指定axis参数沿特定维度排序和计算均值
4. 科学计算生态整合：直接使用SciPy的统计函数
5. 沿轴应用函数：使用apply_along_axis函数
6. 高级索引：使用整数数组索引(花式索引)

实际应用：TCGA皮肤癌数据分析

我们使用癌症基因组图谱(TCGA)的皮肤癌RNA-seq数据来演示分位数归一化的效果。

数据加载与预处理

import bz2
import pandas as pd

# 读取压缩的计数数据
with bz2.open('data/counts.txt.bz2', mode='rt') as f:
    data_table = pd.read_csv(f, index_col=0)

# 转换为NumPy数组
counts = data_table.to_numpy()

归一化前后分布对比

我们定义一个函数来绘制每列的密度分布：

from scipy import stats
import matplotlib.pyplot as plt

def plot_col_density(data):
    """绘制每列的密度曲线"""
    density_per_col = [stats.gaussian_kde(col) for col in data.T]
    x = np.linspace(np.min(data), np.max(data), 100)
    
    fig, ax = plt.subplots()
    for density in density_per_col:
        ax.plot(x, density(x))
    ax.set_xlabel('数据值(每列)')
    ax.set_ylabel('密度')

归一化前分布

log_counts = np.log(counts + 1)
plot_col_density(log_counts)

归一化前，不同样本的基因表达分布差异明显，峰值位置可能相差一个数量级。

归一化后分布

log_counts_normalized = quantile_norm_log(counts)
plot_col_density(log_counts_normalized)

归一化后，所有样本的分布几乎完全相同，消除了技术变异带来的系统性差异。

双聚类分析

完成归一化后，我们可以对基因和样本进行双聚类分析：

1. 基因聚类：发现共表达的基因模块
2. 样本聚类：识别具有相似表达谱的样本群体

为提高效率，我们首先选择变异最大的1500个基因：

def most_variable_rows(data, *, n=1500):
    """选择变异最大的n行"""
    rowvar = np.var(data, axis=1)
    sort_indices = np.argsort(rowvar)[-n:]
    return data[sort_indices, :]

结语

分位数归一化是基因表达数据分析中不可或缺的步骤。通过本文的实现，我们不仅标准化了数据，还展示了如何高效使用NumPy和SciPy进行科学计算。这种归一化方法为后续的差异表达分析、聚类分析等提供了可靠的数据基础。

实际应用中，分位数归一化后的数据可以用于：

• 样本分类和预后预测
• 基因共表达网络构建
• 生物标志物发现
• 多种组学数据的整合分析

理解并掌握这一技术将大大提升您在生物信息学数据分析中的能力。