📚 课程目录第3章数据预处理

数据预处理：标准化与对数转换

拿到原始表达矩阵后，第一步不是跑差异分析，更不是扔给批次校正工具。你得先问问自己：这些数据能不能直接比？

我早期做RNA-seq项目时，接过一份同事的counts矩阵，二话不说就丢进ComBat跑了。结果下游验证怎么也重复不上。后来才发现，样本间总reads数差了3倍，library size都没校正，校正结果全是假的。所以今天这章，我带你把三大核心操作讲透：标准化（CPM/TPM/FPKM）、log1p转换、质量控制过滤。这三步做好了，后面批次检测才不会出幺蛾子。

为什么需要标准化？

说白了，原始counts受两个因素影响：测序深度和基因长度。不同样本总reads数不一样，长的基因更容易被检测到。我们要消除这些技术变异，保留真正的生物学信号。

你想想看，两个样本同一个基因多50%的reads，到底是表达真的高了，还是这组测得深？标准化就是干这个的。

标准化方法	消除因素	适用场景
CPM (Counts Per Million)	测序深度	样本间library size差异不大时
TPM (Transcripts Per Million)	测序深度 + 基因长度	跨基因、跨样本比较推荐
FPKM/RPKM	测序深度 + 基因长度	经典方法，已逐渐被TPM取代

CPM：最简单，也最容易埋坑

CPM = (基因counts / 总counts) × 1e6。就一行公式。我最开始觉得这个太粗暴了，结果在项目中反而常被要求用，因为很多下游工具（比如edgeR的cpm函数）默认就是这个。

代码实现很简单：

import pandas as pd
import numpy as np

def cpm_normalize(counts_matrix):
    # counts_matrix: 基因×样本的DataFrame
    col_sums = counts_matrix.sum(axis=0)
    cpm = counts_matrix.div(col_sums, axis=1) * 1e6
    return cpm

避坑警告： 如果某个样本的总counts非常低（比如< 1e6），CPM会放大噪声。我曾经在单细胞数据上用过CPM，结果导致一批低质量细胞被高表达基因淹没……后来加了QC过滤才稳住。

TPM：更科学的长度校正

TPM是现在的主流。它先用基因长度标准化counts（变成RPK），再除总RPK乘以1e6。好处是：所有样本的TPM之和恒为1e6，跨样本比较更公平。

需要基因长度信息，一般从GTF文件提取。我这里假设你已经有了一个gene_lengths Series：

def tpm_normalize(counts_matrix, gene_lengths):
    # gene_lengths: 基因名对应长度（bp）
    # 第一步：RPK = counts / (length/1000)
    rpk = counts_matrix.div(gene_lengths / 1000, axis=0)
    # 第二步：TPM = RPK / sum(RPK per sample) * 1e6
    tpm = rpk.div(rpk.sum(axis=0), axis=1) * 1e6
    return tpm

记得检查基因长度列是否与counts矩阵行名对齐。我踩过这个坑：GTF文件里的基因ID是ENSEMBL，而矩阵装的是gene_symbol，一乘就全NaN了。

FPKM vs TPM：差在哪？

FPKM是RNA-seq早期经典方法。公式是 RPK / (总 reads / 1e6)。TPM是 RPK / (总 RPK / 1e6)，区别只在分母：FPKM用总reads数，TPM用总RPK值。TPM的均值在不同样本中更稳定，更适合批次校正前的标准化。我个人建议：能用TPM就别用FPKM**。

log1p转换：让数据更像正态

标准化之后，表达值还是严重右偏：大部分基因低表达，少数高表达。直接做线性模型会出问题。log1p = log2(expression + 1)，目的是处理零值（log2(0)是-inf）。

我记得第一次做PCA时没做log1p，前两个主成分全是library size差异，批次效应完全被掩盖了。log转换后才能真正看到生物学结构。

def log1p_transform(df):
    return np.log2(df + 1)

为什么加1不是加0.5？没有绝对标准。我习惯加1，因为log2(0+1)=0，零值保留为0，解释起来直观。注意：如果是TPM这种已经标准化过的，加0.5或1都可以，差别不大。

质量控制过滤：别让垃圾数据拖后腿

过滤是预处理中最容易被跳过的步骤。实际上，原始counts矩阵里可能有大量未表达基因（所有样本counts都为零），或者低表达基因（平均counts < 1）。这些基因没有统计能力，留着只会增加多重检验负担。

典型过滤策略：

去除所有样本中counts全为零的基因
保留在至少N个样本中counts ≥ X的基因（N一般设为最小分组样本数，X常用5或10）
可选：去除表达量极低的样本（总counts < 1e5）

举个具体的：

def filter_low_expression(counts_matrix, min_count=5, min_samples=3):
    # 对每个基因，统计在多少样本中counts ≥ min_count
    mask = (counts_matrix >= min_count).sum(axis=1) >= min_samples
    return counts_matrix.loc[mask]

实战经验： 我习惯先用 min_count=10, min_samples=3 过滤一遍，然后观察过滤后基因数。如果原始2万个基因只剩不到1万，我会调低min_count到5，保证有足够基因用于下游批次检测。别过度过滤，否则信号都没了。

完整预处理流程（SVG图解）

下面这张流程图画出了本章知识结构：从原始counts矩阵到可用于批次检测的表达矩阵。

集成到一个函数中

下面是我项目中常用的一个预处理函数，把CPM/TPM、log1p和过滤打包在一起。你可以按需替换标准化方法：

def preprocess_expression(counts_matrix, gene_lengths=None,
                         method='cpm', log_transform=True,
                         min_count=5, min_samples=3):
    """
    counts_matrix: pd.DataFrame, 基因×样本
    gene_lengths: pd.Series, 基因长度, TPM/FPKM需要
    method: 'cpm', 'tpm', 'fpkm'
    """
    df = counts_matrix.copy()

    # 1. 标准化
    if method == 'cpm':
        col_sums = df.sum(axis=0)
        df = df.div(col_sums, axis=1) * 1e6
    elif method == 'tpm':
        if gene_lengths is None:
            raise ValueError("TPM需要提供gene_lengths")
        rpk = df.div(gene_lengths / 1000, axis=0)
        df = rpk.div(rpk.sum(axis=0), axis=1) * 1e6
    elif method == 'fpkm':
        if gene_lengths is None:
            raise ValueError("FPKM需要提供gene_lengths")
        rpk = df.div(gene_lengths / 1000, axis=0)
        col_sums = df.sum(axis=0)
        df = rpk.div(col_sums, axis=1) * 1e9
    else:
        raise ValueError("method must be 'cpm','tpm','fpkm'")

    # 2. log1p转换
    if log_transform:
        df = np.log2(df + 1)

    # 3. 质量控制过滤：先过滤低表达基因
    mask = (df >= (1 if log_transform else min_count)).sum(axis=1) >= min_samples
    df = df.loc[mask]

    # 如果没做log转换，仍然基于原始counts过滤
    # 这里简化：统一用标准化后的数值判断，实践中可以提前过滤
    return df

核心原则总结：
计数型数据（counts）：必须标准化，推荐TPM>CPM>FPKM。
连续型数据（如微阵列）：通常直接log2转换，不需要长度校正。
过滤：先过滤再跑批次检测，否则噪声会干扰校正。
顺序铁律：原始counts → 标准化 → log转换 → 过滤 → 输出。别改顺序。

好，数据已经洗干净了。接下来就可以放心地在干净的矩阵上做批次检测了。下一章？不，你先把这章代码跑通，检验一下自己处理的数据是什么分布。动手试试。