2. 组学数据基础与工程化读取

组学数据融合建模的第一步，不是建模，而是把四类组学数据“读进来、对齐好、能批量跑”。基因组、转录组、蛋白质组、代谢组在数据形态上差异很大：基因组以区间和变异为核心，转录组以表达矩阵和剪接异构体为核心，蛋白质组以肽段—蛋白映射和丰度为核心，代谢组以质荷比、保留时间和代谢物鉴定为核心。工程化读取的目标，是把这些异构数据统一成“样本 × 特征”的数值矩阵，并保留可追溯的注释与元数据。

四类组学的数据格式与读取要点

基因组数据常见格式包括 FASTA（参考序列）、FASTQ（原始测序读段）、BAM/CRAM（比对结果）、VCF（变异位点）和 BED（区间注释）。工程化读取时，FASTQ 通常先做质控（如 fastp），再用 BWA 或 STAR 比对生成 BAM，最后用 bcftools 或 GATK 得到 VCF。VCF 的关键字段是 CHROM、POS、REF、ALT、QUAL、FILTER、INFO 和 FORMAT，其中 GT 表示基因型，DP 表示深度，AD 表示等位深度。读取 VCF 时要注意：多等位位点需要拆分，缺失基因型（./.）要显式处理，否则后续矩阵会出现非数值列。

转录组数据常见格式包括 FASTQ、BAM、表达矩阵（counts/TPM/FPKM）和 GTF/GFF 注释。工程化读取的核心是“定量—注释—矩阵化”。例如用 featureCounts 或 Salmon 得到 counts 矩阵后，行是基因或转录本，列是样本。读取时需注意：counts 是整数，TPM 是浮点；不同工具对多映射读段的处理不同；基因 ID 与样本 ID 必须与元数据表严格一致。若使用 TPM，需明确其不可直接用于差异表达统计，应使用 counts 配合 DESeq2 或 edgeR。

蛋白质组数据常见格式包括 RAW（质谱原始文件）、MGF（峰列表）、mzML（开放峰列表）、pepXML 和 proteinGroups.txt（MaxQuant 输出）。工程化读取的关键是“肽段—蛋白映射”和“缺失值处理”。proteinGroups.txt 中，Intensity、LFQ intensity、iBAQ 是常用丰度列，Reverse、Potential contaminant、Only identified by site 是过滤列。读取时应先过滤反向库和污染物，再对 LFQ 强度取 log2，最后处理缺失值（如用最小值填充或 KNN 插补）。

代谢组数据常见格式包括 mzML、mzXML、CSV（峰表）和 nmrML。峰表通常包含 m/z、retention time、intensity、sample、metabolite 等列。工程化读取时，需先做峰对齐和保留时间校正，再按代谢物 ID 汇总。代谢组数据缺失值多、批次效应强，读取阶段就应记录批次列，便于后续用 ComBat 或 limma 校正。

公共数据库与注释来源

工程化读取离不开公共数据库。基因组常用 Ensembl、UCSC、NCBI RefSeq；转录组常用 GEO、ArrayExpress、TCGA；蛋白质组常用 PRIDE、ProteomeXchange、MaxQuant 的 Contaminants 库；代谢组常用 HMDB、MetaboLights、KEGG、MassBank。实际读取时，建议把数据库版本号写入元数据，例如 ensembl_release=110、hmdb_version=5.0，否则不同批次注释不一致会导致特征错位。

一个常见的工程化做法是建立“注释映射表”，至少包含：feature_id、source_db、source_id、gene_symbol、uniprot_id、hmdb_id、chrom、start、end。这样四类组学可以按 gene_symbol 或 uniprot_id 做跨组学对齐。

批量加载策略与示例

当样本量达到数百甚至上千时，逐文件读取会非常慢。推荐策略是：先扫描目录生成文件清单，再按样本分块读取，每块读取后立即做类型转换和缺失标记，最后写入 Parquet 或 HDF5。下面是一个可参考的 Python 示例，用于批量读取转录组 counts 矩阵并合并：

import pandas as pd
from pathlib import Path

data_dir = Path("data/rnaseq")
files = sorted(data_dir.glob("*.counts.tsv"))
chunks = []

for f in files:
    sample_id = f.stem.replace(".counts", "")
    df = pd.read_csv(f, sep="\t", index_col=0)
    df.columns = [sample_id]
    chunks.append(df)

expr = pd.concat(chunks, axis=1)
expr = expr.fillna(0).astype("float32")
expr.to_parquet("cache/rnaseq_counts.parquet")
print(expr.shape)  # 例如 (20000, 500)

对于蛋白质组，读取 proteinGroups.txt 时建议先过滤：Reverse != "+"、Potential contaminant != "+"、Only identified by site != "+"，再取 LFQ intensity 列，取 log2 后写入缓存。对于代谢组，峰表读取后先按 metabolite 聚合，再按 sample 透视成矩阵，同时保留 batch 列。

读取质量检查与常见问题

读取完成后，至少检查四项：样本数是否与元数据一致、特征数是否在合理范围、缺失比例是否过高、数值分布是否异常。例如转录组 counts 矩阵中若某样本总 counts 低于其他样本一个数量级，可能是测序深度不足；蛋白质组若某样本 LFQ 缺失率超过 50%，需要标记为低质量样本；代谢组若某批次保留时间整体偏移，说明峰对齐需要重做。

工程化读取的最终交付物，不是一堆原始文件，而是一份可复现的“数据加载清单”：记录每个组学的输入路径、读取函数、过滤条件、注释版本、输出缓存路径和校验值（如行数、列数、缺失率）。这样在融合建模阶段，任何一次矩阵变化都能追溯到读取环节，避免“数据读错但模型跑通”的隐性错误。