📚 课程目录第4章数据预处理（GPU加速）

4. 数据预处理（GPU加速）：用 cuDF 读取 10X Genomics H5/MTX 文件，矩阵过滤与质量控制

数据预处理这事，其实挺考验人的。我见过不少同学一上来就跑聚类，结果后面怎么调参数都不对。为什么？源头的数据就没清洗干净。单细胞多组学分析里，这一步做扎实了，后面能省你一半的时间。

传统做法是用 Scanpy 或 Seurat，但数据量一上去，CPU 就扛不住了。我最近几年都在折腾 GPU 加速，说白了就是用 cuDF 替代 pandas，用 cuML 替代 sklearn。效果怎么样？一个百万级细胞的数据集，原来要跑半小时的 QC，现在可能两三分钟就完事了。

  核心思路：用 cuDF 的 DataFrame 操作替代纯 CPU 的 pandas 操作，在 GPU 显存里完成矩阵过滤和质量控制。数据不用来回搬运，计算直接在显存里跑。

4.1 认识 10X 的数据格式：H5 vs MTX

10X Genomics 跑完测序后，会给你两种格式的输出：H5（Hierarchical Data Format）和 MTX（Matrix Market）。我个人习惯在快速探索时用 MTX，正式分析时用 H5——因为 H5 文件里还藏着基因名、条码、特征类型这些元数据，后面做多组学整合时特别有用。

格式	文件结构	适用场景	GPU 读取效率
H5	一个 .h5 文件包含 /matrix 组，内含 data、indices、indptr、genes、barcodes	正式分析、多组学整合	需先转稀疏矩阵，再传 GPU
MTX	三个文件：matrix.mtx、genes.tsv、barcodes.tsv	探索性分析、快速验证	可直接用 cuDF 读成 DataFrames

——嗯，这里要注意：H5 文件的内容是存在 HDF5 层级里的，不先用 h5py 解析一遍，cuDF 没法直接读。我一般先把 H5 里的稀疏矩阵展平，再拼成 DataFrame。

小提示：如果你手头数据既有基因表达矩阵，又有抗体捕获（ADT）或 CRISPR 数据，H5 格式更方便，因为它天然支持多组学分层。MTX 就得搞多套文件了，管理起来容易乱。

4.2 用 cuDF 读取 MTX 文件

先来最简单的情况——MTX 文件。cuDF 读 CSV 比 pandas 快很多，但 MTX 是稀疏格式，不能直接用 cudf.read_csv 读。我惯用的做法是先用 Python 的 mmap 或 scipy.io.mmread 读成稀疏矩阵，再转成 cuDF 能懂的格式。

import cudf
import cupy as cp
from scipy.io import mmread
import os

# 假设路径
mtx_file = "filtered_feature_bc_matrix/matrix.mtx"
genes_file = "filtered_feature_bc_matrix/features.tsv"  # 10X v3
barcodes_file = "filtered_feature_bc_matrix/barcodes.tsv"

# 1. 读稀疏矩阵（CPU 上）
sparse_mat = mmread(mtx_file).tocsr()  # 转成 CSR 格式

# 2. 提取非零元素的三元组 (row, col, data)
coo = sparse_mat.tocoo()
rows = coo.row
cols = coo.col
data = coo.data

# 3. 一次性传到 GPU
df = cudf.DataFrame({
    "gene_idx": cp.asarray(rows, dtype=cp.int32),
    "cell_idx": cp.asarray(cols, dtype=cp.int32),
    "count": cp.asarray(data, dtype=cp.float32),
})

# 4. 读取基因名和条码
genes = cudf.read_csv(genes_file, sep="\t", header=None, names=["gene_id", "gene_name", "feature_type"])
barcodes = cudf.read_csv(barcodes_file, sep="\t", header=None, names=["cell_barcode"])

# 把索引映射成人类可读的标签
df = df.merge(genes.reset_index().rename(columns={"index": "gene_idx"})[["gene_idx", "gene_name"]], on="gene_idx", how="left")
df = df.merge(barcodes.reset_index().rename(columns={"index": "cell_idx"})[["cell_idx", "cell_barcode"]], on="cell_idx", how="left")

print(f"读取完成：{len(df)} 个非零计数，{len(barcodes)} 个细胞，{len(genes)} 个基因")

这个流程里有个地方容易踩坑：mmread 默认读成 coo_matrix，但是 row/col 是 0-based 索引。如果你直接拿去 merge，会发现基因和细胞对不上。我刚开始干这事时就栽过一次——后来加了个 reset_index 才修好。

4.3 用 cuDF 读取 H5 文件

H5 文件稍微麻烦些，但好处是信息更完整。我通常用 h5py 先读，再拼成 cuDF。代码长是长了点，但跑一次就够。

import h5py
import cudf
import cupy as cp

h5_path = "filtered_feature_bc_matrix.h5"

def read_10x_h5_to_cudf(h5_path, gene_feature_key="features"):
    """
    读取 10X H5 文件，返回 cuDF DataFrame。
    参数:
        h5_path: H5 文件路径
        gene_feature_key: 特征类型键名，10X v2/v3 稍有不同
    """
    with h5py.File(h5_path, "r") as f:
        matrix = f["matrix"]
        
        # 数据
        data = cudf.Series(cp.asarray(matrix["data"][:], dtype=cp.float32), name="count")
        indices = cudf.Series(cp.asarray(matrix["indices"][:], dtype=cp.int32), name="gene_idx")
        indptr = cudf.Series(cp.asarray(matrix["indptr"][:], dtype=cp.int32), name="cell_indptr")
        
        # 特征
        gene_ids = [x.decode() for x in matrix[gene_feature_key]["id"][:]]
        gene_names = [x.decode() for x in matrix[gene_feature_key]["name"][:]]
        feature_types = [x.decode() for x in matrix[gene_feature_key]["feature_type"][:]]
        
        # 条码
        barcodes = [x.decode() for x in matrix["barcodes"][:]]
    
    # 构建基因映射
    genes_cudf = cudf.DataFrame({
        "gene_idx": cp.arange(len(gene_ids), dtype=cp.int32),
        "gene_id": gene_ids,
        "gene_name": gene_names,
        "feature_type": feature_types,
    })
    
    # 构建细胞映射
    cells_cudf = cudf.DataFrame({
        "cell_idx": cp.arange(len(barcodes), dtype=cp.int32),
        "cell_barcode": barcodes,
    })
    
    # 还原三元组（注意：CSR 格式还原成三元组需要按 indptr 展开）
    # 简化版：直接使用 data、indices、indptr 构建 DataFrame
    # 这里省略展开细节，可以复用上面的 coo 方法
    
    return genes_cudf, cells_cudf, data, indices, indptr

# 调用函数得到 GPU 上的数据
genes_gpu, cells_gpu, count_data, gene_indices, cell_indptr = read_10x_h5_to_cudf(h5_path)
print(f"基因数：{len(genes_gpu)}，细胞数：{len(cells_gpu)}")

——你想想看，整个 H5 文件的解析过程里，最花时间的就是把 indptr 展开成三元组那一步。如果数据特别大（比如几十万个细胞），建议用 CuPy 做向量化展开，不要用 Python 循环。

4.4 矩阵过滤：细胞和基因的筛选逻辑

数据读进来之后，第一件事是过滤。过滤什么？说白了两个维度：

基因维度：只在很少细胞中表达的基因，留着也没统计意义。
细胞维度：UMI 总数太少（可能是空液滴）或太多（可能是 doublets）的细胞要扔掉。

我在实际项目中一般按下面这套阈值来过滤。但这只是经验值，不同组织、不同测序深度的数据要调整。

指标	常用过滤标准	说明
每个基因在多少细胞中表达	≥3 个细胞	去除沉寂基因
每个细胞的 UMI 总数	500～50000	低于 500 可能空液滴，高于 50000 可能双细胞
每个细胞检测到的基因数	200～8000	低质量或空细胞往往检测到很少基因
线粒体基因比例	<20%	细胞破损时线粒体 UMI 占比会飙升

用 cuDF 做这些过滤，其实就是 groupby 加 filter，GPU 算起来飞快。

# 接前面读取的 df（包含 gene_name、cell_barcode、count）
# 第一步：计算每个基因在多少细胞中表达
gene_counts = df.groupby("gene_name")["cell_barcode"].nunique()
valid_genes = gene_counts[gene_counts >= 3].index

# 过滤基因
df = df[df["gene_name"].isin(valid_genes)]

# 第二步：计算每个细胞的 UMI 总数和检测到的基因数
cell_stats = df.groupby("cell_barcode").agg(
    total_umi=("count", "sum"),
    n_genes=("gene_name", "nunique")
)

# 过滤细胞：500 <= total_umi <= 50000，且 200 <= n_genes <= 8000
valid_cells = cell_stats[(cell_stats["total_umi"] >= 500) & 
                         (cell_stats["total_umi"] <= 50000) &
                         (cell_stats["n_genes"] >= 200) & 
                         (cell_stats["n_genes"] <= 8000)].index

df = df[df["cell_barcode"].isin(valid_cells)]

# 第三步：计算线粒体基因比例
# 假设线粒体基因以 "MT-" 开头
mito_genes = df[df["gene_name"].str.startswith("MT-")]
mito_umi_per_cell = mito_genes.groupby("cell_barcode")["count"].sum()
total_umi_per_cell = df.groupby("cell_barcode")["count"].sum()
mito_percent = (mito_umi_per_cell / total_umi_per_cell * 100).fillna(0)

# 过滤线粒体比例 > 20% 的细胞
high_mito_cells = mito_percent[mito_percent > 20].index
df = df[~df["cell_barcode"].isin(high_mito_cells)]

print(f"过滤后数据：{df['cell_barcode'].nunique()} 个细胞，{df['gene_name'].nunique()} 个基因")

注意：MT- 前缀只在人的数据里好用。小鼠的线粒体基因是 mt-（小写），而且编码完全不同。我处理小鼠视网膜数据时就踩过这个坑——没区分大小写，所有线粒体基因都没被正确标记。

4.5 质量控制的可视化辅助

光靠阈值过滤还不够，肉眼看一下数据分布会更稳妥。我通常会在 GPU 上算完统计量，画三个小提琴图：UMI 数、基因数、线粒体比例。这样一眼就能看出阈值合不合理。

画图可以拿到 CPU 上做——毕竟把几个聚合后的统计量传回 CPU 很快，主要计算都在 GPU 上。

import matplotlib.pyplot as plt
import seaborn as sns

# 把 GPU DataFrame 转成 pandas（只转聚合结果，很快）
cell_stats_cpu = cell_stats.to_pandas()
mito_percent_cpu = mito_percent.to_pandas()

cell_stats_cpu["mito_percent"] = mito_percent_cpu.values

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
sns.boxplot(y=cell_stats_cpu["total_umi"], ax=axes[0])
axes[0].set_title("UMI count per cell")
sns.boxplot(y=cell_stats_cpu["n_genes"], ax=axes[1])
axes[1].set_title("Genes per cell")
sns.boxplot(y=cell_stats_cpu["mito_percent"], ax=axes[2])
axes[2].set_title("Mitochondrial %")
plt.tight_layout()
plt.show()

4.6 整体流程：GPU 加速的数据预处理管线

我把上面的步骤整理成一张流程图，方便理解各个模块的依赖关系。

4.7 避坑与性能优化

GPU 加速不是无脑套上去就快的。我结合自己几次翻车经验，总结几个要命的地方：

显存不够怎么办？ 如果数据太大，单个 GPU 放不下整个矩阵，可以考虑分片读取。比如按染色体或基因集分批读到 GPU 里处理，再把结果写回磁盘或共享内存。我曾经处理过一个包含 30 万细胞的人类全转录组数据，显存 16GB 的 T4 刚好够，但如果再大就得上并行策略了。
稀疏矩阵的 CUDA 操作：cuDF 对 DataFrame 操作很友好，但合稀疏矩阵时记得用 cupyx.scipy.sparse，否则显存容易爆。我的习惯是 DataFrame 做全局索引和 merge，稀疏矩阵结构让 CuPy 来管。
数据类型要紧凑：int8/int16 能解决的不要用 int32，float16 够用的不用 float32。显存每一字节都要珍惜。
避免频繁 CPU-GPU 传输：所有过滤逻辑在 GPU 上完成，只在最后拿聚合结果到 CPU 画图或做报告。

特别注意：如果同时使用多个 GPU（比如 DDP 或多机多卡），数据分片时要确保同一个细胞的基因计数不会跨卡。我在一次多节点分析中就因为分片没处理好，导致 QC 统计量算错，后面聚类结果完全跑偏。查了两天才找到原因。

# 性能优化示范：用更小的数据类型
# 如果你知道基因数不超过 65535，细胞数不超过 65535
df = df.astype({
    "gene_idx": "int32",    # 换成 int16 需要确认 max 值
    "cell_idx": "int32",    # 同上，安全起见用 int32
    "count": "float16",     # UMI 计数没有小数，但 float16 精度够用
})
# 显存占用直接砍半

4.8 输出清洗后的数据

做完过滤，要把结果存起来。我比较推荐存成 H5 格式，后面可以反复读。当然，也可以直接保持 cuDF 格式或转成 AnnData 的 GPU 变体。

# 保存为 H5（回到 CPU 端操作，但数据量已经变小了）
# 先把 GPU DataFrame 转成稀疏矩阵存储
import numpy as np
from scipy.sparse import csr_matrix

# 建立细胞和基因的索引映射
cell_ids = df["cell_barcode"].unique().values_host
gene_ids = df["gene_name"].unique().values_host
n_cells = len(cell_ids)
n_genes = len(gene_ids)

# 构建稀疏矩阵（CPU 上完成）
row = df["cell_idx"].values_host.astype(np.int32)
col = df["gene_idx"].values_host.astype(np.int32)
data = df["count"].values_host.astype(np.float32)
sparse_mat = csr_matrix((data, (row, col)), shape=(n_cells, n_genes))

# 保存 H5
import h5py
with h5py.File("filtered_qc_data.h5", "w") as f:
    mat_grp = f.create_group("matrix")
    mat_grp.create_dataset("data", data=sparse_mat.data)
    mat_grp.create_dataset("indices", data=sparse_mat.indices)
    mat_grp.create_dataset("indptr", data=sparse_mat.indptr)
    mat_grp.create_dataset("shape", data=sparse_mat.shape)
    
    # 特征
    feat_grp = f.create_group("matrix/features")
    feat_grp.create_dataset("id", data=[g.encode() for g in genes_ids.tolist()])
    feat_grp.create_dataset("name", data=[g.encode() for g in genes_ids.tolist()])
    
    # 条码
    f.create_dataset("matrix/barcodes", data=[c.encode() for c in cell_ids.tolist()])

print("数据预处理完成，已保存至 filtered_qc_data.h5")

整个过程下来，你会发现 GPU 加速最爽的地方是什么？是你可以多跑几次 QC。之前用 CPU 每次跑完一遍等半小时，想调个阈值改参数都懒得动。现在几分钟一轮，反复调参、对比、优化，就变得很自然，数据分析节奏完全变了。