DIABLO多组学联合建模工程 - 第3节

3. 多组学数据类型与常见格式

3. 多组学数据类型与常见格式：转录组、蛋白组、代谢组等数据的读取与整理

做多组学联合分析，第一步永远不是跑模型，而是把手头的原始数据变成干净、对齐的矩阵。我在不同项目里见过太多“模型跑不动”的案例，最后发现都是数据格式没处理好。这一章咱们就集中精力，把转录组、蛋白组、代谢组最常见的数据格式捋一遍，然后看看怎么用 R 把它们读进来、整理好，直接喂给后面的 DIABLO。

3.1 转录组数据：从 FASTQ 到表达矩阵

转录组数据，说白了最开始都是一堆 FASTQ 文件。但你想想看，我们做多组学联合时很少直接用 FASTQ——那玩意儿太大、太原始。一般会先比对、定量，得到一个“表达矩阵”，也就是基因 × 样本的计数表。常见的格式有：

我个人习惯，在做联合分析前，只保留基因名在行、样本名在列的表达矩阵。别在 R 里读 FASTQ——那是工具链的事。我们在 DIABLO 环节需要用到的数据，通常是 counts、FPKM 或者 TPM 格式。

读取示例

用 data.table::fread 读大矩阵非常快。我之前接过一个项目，转录组矩阵 60 万行 × 200 样本，read.table 读不动，换成 fread 十几秒就搞定了。

# 最基础的方式，适合小文件
expr_mat <- read.csv("rnaseq_counts.csv", row.names = 1)

# 推荐方式：data.table 读取
library(data.table)
expr_mat <- fread("rnaseq_counts.tsv", header = TRUE, sep = "\t")
expr_mat <- as.data.frame(expr_mat)
rownames(expr_mat) <- expr_mat[[1]]
expr_mat <- expr_mat[, -1]  # 删掉基因名列，交给行名

# 检查维度
dim(expr_mat)

3.2 蛋白组数据：质谱原始峰到定量矩阵

蛋白组数据和转录组不太一样。最常见的来源是 LC-MS/MS 质谱分析，经过搜库（比如 MaxQuant、Proteome Discoverer、MSFragger）后，会输出一个蛋白定量结果表。常见的格式包括：

• 蛋白质组学搜库结果（如 MaxQuant 的 proteinGroups.txt）：包含各种定量值 iBAQ、LFQ intensity、MS/MS counts 等。
• 肽段谱图数据（如 mzML、mzXML）：一般在预处理环节，不会直接进联合分析。
• 标准化后的蛋白表达矩阵：蛋白 ID × 样本，数值通常是 log2 转换后的 Intensity。

避坑指南：我曾经有一次从合作伙伴那拿到一个蛋白组结果文件，里面既有“peptide sequences”又有“protein accessions”，行长 20 万行。我只想拿来跑 DIABLO，结果发现同一蛋白对应多行肽段——得先做汇总（取最大值或中位数）才能拿到蛋白级别的表达。

# 读取 MaxQuant 经典输出
prot_raw <- fread("proteinGroups.txt", sep = "\t", header = TRUE)

# 选择需要的列：蛋白 ID（'Protein IDs'） 和样本定量列（如 'LFQ intensity Sample1' 等）
# 注意：同一蛋白可能对应多个肽段，建议保留唯一蛋白
# 这里演示简单过滤
prot_df <- prot_raw[, .(Protein.IDs, Intensity_S1, Intensity_S2)]
rownames(prot_df) <- prot_df$Protein.IDs
prot_mat <- as.matrix(prot_df[, -1])

# 处理缺失值：蛋白组常出现 0 或 NA，建议先过滤掉 50% 以上样本为缺失的蛋白
prot_mat[prot_mat == 0] <- NA
keep <- rowSums(is.na(prot_mat)) / ncol(prot_mat) < 0.5
prot_mat <- prot_mat[keep, ]

3.3 代谢组数据：谱图、峰表和浓度矩阵

代谢组复杂度又高了一截。数据来源可能是 GC-MS、LC-MS 或 NMR，原始文件格式五花八门（.cdf、.mzML、.raw、.d 等）。经过峰检测、对齐、去卷积后，通常会得到一个“峰表”或“特征表”：

我在处理代谢组时有个习惯：拿到峰表后，先确认缺失值比例。代谢组峰表的缺失率经常很高（比如 60%），因为很多代谢物只在部分样本中检出。直接把这些特征删掉会损失很多信息，所以我会用 KNN 或最小值填充（视下游分析要求而定）。

# 读取 XCMS 输出的峰表
peak_table <- fread("xcms_peaks.csv", header = TRUE)
# 通常第一列是特征 ID（如 "M123T456"），其余列为样本峰面积
rownames(peak_table) <- peak_table[[1]]
peak_mat <- as.matrix(peak_table[, -1])

# 缺失值填充：用每列最小值的 1/2 填充
for (j in 1:ncol(peak_mat)) {
    col_min <- min(peak_mat[, j], na.rm = TRUE)
    peak_mat[is.na(peak_mat[, j]), j] <- col_min / 2
}

3.4 三组数据统一整合：最终对齐

DIABLO 要求所有组学数据矩阵样本顺序一致（相同的样本 ID 在行，变量在列）。所以我一般在分别处理后，做一个样本交集筛选。

这里给个经典流程：

1. 分别从每个组学提取样本名（列名或行名）
2. 找到共有样本（一般用 intersect）
3. 按共有样本顺序重排每个矩阵的行
4. 将每个矩阵转置为 n（样本）× p（特征）形式

# 假设已有 expr_mat, prot_mat, peak_mat
# 样本在列，所以先取列名
common <- Reduce(intersect, list(colnames(expr_mat), 
                                 colnames(prot_mat), 
                                 colnames(peak_mat)))
# 对齐
expr_aligned <- expr_mat[, common, drop = FALSE]
prot_aligned <- prot_mat[, common, drop = FALSE]
peak_aligned <- peak_mat[, common, drop = FALSE]

# 转置为样本×变量（DIABLO 期望的格式）
expr_final <- t(expr_aligned)
prot_final <- t(prot_aligned)
peak_final <- t(peak_aligned)

# 用 list 打包，后面直接送 DIABLO
data_list <- list(RNA = expr_final, Protein = prot_final, Metabolite = peak_final)

3.5 注意事项与“我曾经”系列

我在前面挖的坑，再补几个：

• 基因/蛋白 ID 的重复问题： 如果读取后发现有重复的行名（比如一个基因对应多个探针），建议按均值或中位数去重，不要直接删。
• 计数数据的处理： 转录组 counts 是整数，DIABLO 内部默认做标准化，但如果 counts 值范围太大（比如 0～1e6），可以先做 log2(counts+1) 转换，避免极端值主导。
• 代谢组峰表的保留： 不建议把所有特征都保留。先做缺失率筛选（比如保留至少 80% 样本有值的特征），否则后续降维时噪音太多。

这一章的内容就到这里。你如果前面几步都做扎实了，后面用 DIABLO 就是水到渠成的事。好，咱们继续。