掃碼聯系客服
引入：提出問題
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/#271588
https://bioconductor.org/packages/release/bioc/vignettes/oligo/inst/doc/oug.pdf
