芯片表达数据预处理：从探针RMA-based value 到基因水平表达量
引入：提出问题
GEO 里海量的芯片表达数据为许多validation 分析提供了丰富的素材。有时，直接从GEO页面得到的series matrix file 是 normalization 之后的数据用起来不太友好，就需要我们从原始的.cel文件获取基因的表达矩阵来满足我们自己的分析需要。
我们都知道芯片数据是基于一个个探针信号来对基因进行定量的，并且经常出现设计多个探针对应同一个基因进行检测，因此必须将探针的表达量映射回基因水平才能进行分析或结果解释。因此这里涉及一个比较基础的问题：
为了得到基因水平的表达矩阵，我们在处理芯片原始数据时，应该走下面那条路：
先把同一基因的不同探针的原始信号整合起来，再进行定量处理（例如RMA），得到基因水平的表达；先直接按照探针水平进行表达信号的预处理（例如RMA），再对同一个基因的多个探针表达量进行aggregate（例如取均值或者中值）芯片探针水平分析
芯片表达定量过程分为background correct、quantile normalize、summarize 三个主要步骤，可以通过R包 affy 或 oligo 中的一体化算法 RMA（Robust Multiarray Average）直接完成上述四步的数据处理，这个网上有较多的R代码教程，这里不进行赘述。需要注意的是，RMA是针对探针水平的原始数据进行处理，最终得到的表达矩阵仅仅是将荧光强度值从探针水平汇总到探针组（probeset）水平，并非我们平时分析时使用的基因水平的表达谱。
探针组水平的表达矩阵
查了很多关于芯片预处理的教程，发现大多都是说明如何得到探针组的表达矩阵之后点到即止，立马转到下游的差异表达分析了。
设计同一个基因的不同探针是为了对应检测基因不同的转录片段，并且有一些探针可以识别非特异性转录本，因此可以检测多个基因(promiscuous probes)。所以通常大家在最终从具体基因总结分析结果之前，会尽可能长时间地使用探针水平数据进行分析。但是这样一来得到的差异表达“基因”其实是一些差异表达的“探针组”，在结果解释以及可视化时仍然需要对探针ID进行注释，也就是对应到基因ID。
Durham, United States 的高级计算生物学家Michael B Black 称，如果有两个或以上的probes达到的分析的显著性阈值，为了最终的结果解释，可以对同一个基因不同探针组的log2FC取平均。但他同时也强调，在所有的统计分析都完成之前，绝不能将probe整合到基因水平。使用探针水平的数据进行所有的显著性测试分析，基因水平用于最终的结果展示。
But you should never summarize probes to genes before all your statistical anaylses are complete. Use probe level data for all significance testing, and only summarize by gene ID for final presentation.
因此，第一条路貌似走不通，会混淆该基因不同探针的信号，导致定量计算不准确。
但如果我的分析目的并不局限于差异表达，比如想做GSEA这一类必须基于基因表达才能够进行的分析时，必须拿到一个以基因为单位的表达矩阵，那似乎只剩第二条路，可是这么做有点不踏实，看看有没有支持的证据。
芯片探针表达整合到基因水平
1、RMA函数可直接获取基因水平的表达
经过查证，对于基因和外显子芯片（名字中带有'Gene' 或 'Exon'，还可能有 'ST'），函数rma()可通过设置参数target = "core"直接指定获取基因水平的表达。这是因为对于这些芯片，Affymetrix提供了额外的一个注释文件（注释meta-probesets (MPSs)），用于将探针表达整合到基因水平。
# target: Level of summarization (only for Exon/Gene arrays)
geneSummaries <- oligo::rma(affyExonFS, target="core")
对于外显子芯片，rma可以直接定量到外显子水平的表达，此时可以设置target="probeset"来获得探针水平的表达。
probesetSummaries <- rma(affyExonFS, target="probeset")
2、其他芯片
其他芯片缺少这个注释文件只能采取一般的方法，即先进行探针定量，然后再对同一个基因的探针表达量进行aggregate（取均值或者中值）。这个方法也在网上技术大佬们给出的回答中发现。
If you wish to summarise over transcripts after you normalise with rma(), you can do something rudimentary like using aggregate() (base R) or avereps() (limma)
关键步骤代码：
# 读取cel文件
data.raw <- read.celfiles(filenames = cel.files)
# 使用RMA方法进行预处理：依次进行背景处理，log2转化，分位数标化和探针表达量计算
eset.rma <- rma(data.raw)
# 提取探针水平表达矩阵
ex.mat.rma <- exprs(eset.rma)
# 多个探针的表达均值作为基因表达，probe2gene为探针和基因对应信息数据框
gene.ex.mat <- aggregate(ex.mat.rma, by = probe2gene$SYMBOL, FUN = mean)
rownames(gene.ex.mat) <- probe2gene$SYMBOL
PS. 探针ID的注释也有很多人讨论过怎么做，重点在于如何获取探针和基因对应信息。
最后，这个小查证解决了本人在对不同平台芯片数据上手时一个棘手的问题，希望也可以以帮到同样有问题的小伙伴。希望对于GEO里的数据资源大家都可以手到擒来，玩的飞起~
参考资料
https://www.biostars.org/p/271379/#271588https://bioconductor.org/packages/release/bioc/vignettes/oligo/inst/doc/oug.pdfhttps://www.researchgate.net/post/When-I-download-microarray-processed-normalized-data-from-GEO-it-contains-multiple-probe-IDs-against-one-gene-Which-probe-value-I-should-select了解更多生信分析干货，欢迎关注灵魂工具人
芯片表达数据预处理：从探针RMA-based value 到基因水平表达量
相关文章