探针水平的表达矩阵转换为基因水平的表达矩阵是芯片数据分析中关键的一步
将探针水平的表达矩阵转换为基因水平的表达矩阵是芯片数据分析中关键的一步,因为多个探针可能对应同一个基因,需要整合为一个值。以下是主要的转换方法和步骤。
方法 基本原理 特点/常用函数 主要考虑
保留最大值 对同一基因的多个探针,保留表达量(行和或行均值最大)的探针。 能代表该基因最高表达活性。 可能忽略该基因其他亚型或转录本的信息。
取平均值 对同一基因的多个探针取平均表达量。 limma::avereps() 操作简单,得到基因的整体表达趋势。但可能平滑掉有生物学意义的差异。
随机去重 当多个探针对应一个基因时,随机保留其中一个。 仅在快速处理或探针间差异不大时考虑。 缺乏生物学依据,一般不推荐。
🧬 转换前的重要准备
在进行转换前,有两项准备工作至关重要:
-
获取探针注释信息:你需要一个数据框(通常命名为ids或probe2gene),其中至少包含两列:一列是探针ID(与表达矩阵行名匹配),另一列是对应的基因符号(Gene Symbol) 或 Gene ID。
◦ 来源:这通常来自GEO数据库的平台文件(GPL)。你可以通过getGEO()函数下载,或从GEO官网下载并读取注释文件。◦ 处理复杂注释:有些平台的Gene Symbol列可能包含多个基因(如"DDR1 /// MIR4640"),需要进行清理。常见做法是提取第一个基因符号或进行拆分。
示例:从GPL的gene_assignment列提取Symbol
symbol <- trimws(tstrsplit(probe$gene_assignment, “//”, fixed=TRUE)[[2]])
-
确保探针顺序一致:在合并注释与表达矩阵前,确保注释文件中的探针ID顺序与表达矩阵的行名顺序一致,或者通过merge()函数进行匹配。
确保顺序一致并取行和最大的探针
ids <- ids[order(rowSums(exp), decreasing = T), ]
exp_sorted <- exp[idsprobeid,]identical(idsprobe_id, ] identical(idsprobeid,]identical(idsprobe_id, rownames(exp_sorted)) # 检查顺序是否一致
🧪 主要转换方法
-
使用 tinyarray 包中的 trans_array 函数
该函数默认保留每个基因的第一个探针(随机去重)。
library(tinyarray)exp: 探针表达矩阵; ids: 包含probe_id和symbol两列的注释数据框
gene_exp <- trans_array(exp, ids)
-
使用 limma 包中的 avereps 函数取平均值
这是非常常用且稳健的方法。
library(limma)首先将表达矩阵的行名替换为基因名
exp_with_gene <- exp
rownames(exp_with_gene) <- ids$symbol # 确保ids与exp行顺序匹配对重复基因取平均值
gene_exp <- avereps(exp_with_gene)
-
手动处理(如保留最大值)
若想保留表达量最高的探针,可先对注释框按表达量排序后再用trans_array。计算每个探针的行和,并从大到小排序
ids_ordered <- ids[order(rowSums(exp), decreasing = T), ]
使用排序后的ids进行转换,trans_array会保留第一个出现的(即行和最大的)
gene_exp_max <- trans_array(exp, ids_ordered)
📌 注意事项
• 转换时机:建议在完成探针水平的标准化和质量控制后,再进行转换。差异表达分析最好在探针水平进行,将结果转换到基因水平用于解释和展示。
• 缺失值处理:转换前常需移除未注释到基因的探针(NA值)。
• 方法选择:取平均值(avereps) 是最常用且通用的方法。若研究关注特定转录本或亚型,保留最大值可能更合适。
• 特殊平台:对于Affymetrix Gene/Exon ST芯片,使用oligo包的rma()函数时,可设置参数target = "core"直接得到基因水平表达矩阵。
💎 总结工作流程
一个典型的工作流程如下:
- 获取数据:表达矩阵(如从GSE的series_matrix.txt.gz)和探针注释信息(如从GPL)。
- 预处理:清洗注释信息(处理多个基因名、去除NA),确保与表达矩阵探针匹配。
- 转换:选择合适方法(如limma::avereps)将探针矩阵转换为基因矩阵。
- 后续分析:基于基因水平矩阵进行后续分析,如差异表达、聚类、GSEA等。
希望这些信息能帮助你顺利完成转换!