第2章：数据加载——读取基因表达矩阵

好了，开始动手写代码前，我们先搞定数据。搞不定数据，后面所有分析都是空中楼阁。我个人习惯把数据加载这一步称为“地基工程”，地基没打好，房子迟早要塌。

基因表达矩阵最常见的形式就是CSV或TSV。说白了，就是一个二维表格，行是基因，列是样本。你看过那种几万行乘几百列的矩阵吗？我第一次处理的时候，Excel直接卡死了……嗯，从那以后我就只用Python读了。

2.1 用pandas读取CSV/TSV

pandas的read_csv函数是主力工具。别看名字里带个CSV，它读TSV也完全没问题。你只需要指定分隔符参数就行。

先看看基本用法：

import pandas as pd

# 读取CSV文件
df = pd.read_csv('expression_data.csv', index_col=0)

# 读取TSV文件
df_tsv = pd.read_csv('expression_data.tsv', sep='\t', index_col=0)

这里有个关键参数index_col=0。这个参数的意思是：把第一列作为行索引。基因表达矩阵的第一列通常是基因名，比如“ENSG00000123456”或者“TP53”。如果不加这个参数，pandas会给你生成一个从0开始的数字索引，那基因名就变成了一列数据，后面操作起来会很别扭。

有时候数据文件很大。几万个基因，几百个样本，文件轻松上几百MB。这时候直接读可能会崩。我建议你加上几个参数：

df_large = pd.read_csv('big_matrix.csv', 
                       index_col=0,
                       dtype='float32',
                       low_memory=False)

dtype='float32' 能把内存砍半。表达值本来就是浮点数，用float32精度足够了。low_memory=False 防止pandas分块读取时类型推断出错。

df_preview = pd.read_csv('big_matrix.csv', nrows=100)
df_preview.head()

2.2 处理缺失值

现实世界的数据没有完美的。表达矩阵里经常出现NaN或0值。为什么会有缺失？可能是某个基因在某些样本里没检测到，或者检测技术本身有局限。

我个人处理缺失值时，会先看缺失比例。如果比例太高，这个基因或样本可以直接扔掉。如果只是零星缺失，可以用填补的方式。

先看看缺失值长什么样：

# 检查缺失值
print(df.isnull().sum())

# 按行（基因）统计缺失比例
missing_ratio = df.isnull().mean(axis=1)
print(missing_ratio.sort_values(ascending=False).head(10))

isnull() 返回布尔矩阵，sum() 按列求和，mean() 按行求比例。这样你就能一眼看出哪些基因缺失严重。

接下来是处理策略，我一般用两种：

注意，填补均值是按行（基因）来算的。每个基因有自己的表达水平，不能用全局均值去填，那会引入偏差。

我曾经处理过一个单细胞RNA-seq数据集，发现有一半的基因缺失率超过80%。那时候我选择直接删除这些低表达基因，因为它们在后面的分析里也提供不了什么信息，反而会引入噪音。

完整代码示例：

# 删除缺失超过50%的基因
df_clean = df.dropna(thresh=len(df.columns) * 0.5, axis=0)

# 对剩余基因填补均值（按行计算）
df_clean = df_clean.T.fillna(df_clean.mean(axis=1)).T

这里我用了两次转置（.T），是因为fillna是按列操作的。转置后，每个基因的表达值变成一列，就可以用该列的均值去填补缺失了。填完再转置回来。

2.3 处理重复基因

重复基因也是个常见问题。为什么会有重复？可能是同一个基因对应多个探针，或者数据整合时从不同平台来的，或者注释版本不一致造成的。

我建议的处理步骤：

• 先找出哪些基因重复了。
• 然后决定保留哪个值。一般用均值或中位数。
• 最后删除重复项。

看代码：

# 找出重复的基因名
duplicated_genes = df_clean.index[df_clean.index.duplicated(keep=False)]
print(f"重复基因数量: {len(duplicated_genes)}")

# 按基因名分组，取均值
df_dedup = df_clean.groupby(df_clean.index).mean()

这里keep=False会标记所有重复项（不只是第一次出现）。groupby().mean() 把同一基因的多行合并成一行，值取均值。

你可能会问：取均值好还是中位数好？我的经验是：如果数据比较对称，用均值；如果可能有离群值，用中位数更稳健。不过一步到位的做法是用np.median：

import numpy as np

# 按基因名分组，取中位数
df_dedup = df_clean.groupby(df_clean.index).agg(np.median)

还有一种情况：同一个基因在不同批次里用不同的ID表示。比如，“TP53”和“p53”其实是同一个东西。这时候需要先统一基因ID，再做去重。这种场景我在整合TCGA和GEO数据时遇到过多次。

2.4 数据加载的完整流程

把上面几步串起来，就是一份完整的数据加载脚本。我通常用一个函数封装：

def load_expression_matrix(filepath, sep=',', missing_thresh=0.5):
    """
    加载表达矩阵，处理缺失值和重复基因。
    
    参数:
        filepath: 文件路径
        sep: 分隔符，CSV用','，TSV用'\t'
        missing_thresh: 基因允许的最大缺失比例
    
    返回:
        清洗后的DataFrame
    """
    # 1. 读取数据
    df = pd.read_csv(filepath, index_col=0, sep=sep)
    print(f"原始数据: {df.shape[0]} 基因, {df.shape[1]} 样本")
    
    # 2. 删除缺失过多的基因
    min_count = int(len(df.columns) * (1 - missing_thresh))
    df = df.dropna(thresh=min_count, axis=0)
    print(f"删除缺失后: {df.shape[0]} 基因")
    
    # 3. 填补剩余缺失值（按行均值）
    df = df.T.fillna(df.mean(axis=1)).T
    
    # 4. 处理重复基因（取均值）
    df = df.groupby(df.index).mean()
    print(f"去重后: {df.shape[0]} 基因")
    
    return df

调用起来很简单：

expr = load_expression_matrix('my_data.csv', sep=',')
print(expr.head())

2.5 知识结构图

下面这张图总结了数据加载的核心流程。从原始文件到干净矩阵，中间要经过缺失值处理和重复基因合并两个关键步骤。

看到这个流程了吗？每一层都在清洗数据。你别嫌麻烦，数据质量决定了分析结果的上限。脏数据进，垃圾结果出。这个道理我在多少个项目里验证过。

好了，数据加载这块就讲到这里。你的表达矩阵现在应该干干净净、整整齐齐，可以拿去喂给后续的批次效应检测算法了。