昨天發(fā)了一個差點捅了馬蜂窩的推文：費九牛二虎之力也無法重現(xiàn)的GSEA圖 ，當(dāng)時的確花了兩三個小時也沒能成功重復(fù)文章的GSEA圖，覺得很有趣，就跟粉絲分享了，純粹是交流數(shù)據(jù)處理技巧。
后來跟文章作者充分交流了，也理解到了他們研究的亮點，并不是在這個簡單的芯片數(shù)據(jù)上面，而且不同的數(shù)據(jù)處理方法的選擇，本身就會有不同的結(jié)果，但并不會影響他們文章的重大發(fā)現(xiàn)。
正是因為跟第一作者充分的溝通，才了解到為什么當(dāng)時他們的芯片數(shù)據(jù)沒有生物學(xué)重復(fù)，作者的解釋是：樣本數(shù)量的選擇根據(jù)研究的內(nèi)容、目的和芯片對于文章的重要性來看。他們是從小鼠淋巴結(jié)中分離的基質(zhì)細(xì)胞，從小鼠的淋巴結(jié)中獲取足夠量的基質(zhì)細(xì)胞來做芯片，每組用了50只小鼠，是個混合樣本。此芯片的目的是對功能學(xué)行為進(jìn)行佐證，沒有挑選關(guān)鍵基因進(jìn)行后續(xù)的機制研究，故而沒有進(jìn)行多組樣本的檢測。
那我們現(xiàn)在就繼續(xù)討論一下，文中的GSEA的顯著性結(jié)果圖（下圖）到底該如何來重現(xiàn)。
請大家再注意看圖中的基因集，應(yīng)該是50個基因左右的數(shù)量，我昨天主要是就糾結(jié)在這里，因為我在不同版本的GO數(shù)據(jù)庫資源，都找不到50個基因左右的數(shù)量的這兩個通路，大家也可以自行去數(shù)據(jù)庫搜索了解這兩個GO基因集：
# GO:0002690 Positive regulation of leukocyte chemotaxis.
# http://amigo.geneontology.org/amigo/term/GO:0030595
# Accession GO:0030595
# Name leukocyte chemotaxis
# Ontology biological_process
# Synonyms immune cell chemotaxis, leucocyte chemotaxis
# http://amigo.geneontology.org/amigo/term/GO:0006959
# Accession GO:0006959
# Name humoral immune response
# Ontology biological_process
不懂就要問，所以我直接和文章的第一作者顧博士聊了一下這個數(shù)據(jù)問題，了解到了他們完成研究課題的艱辛和付出。只是術(shù)業(yè)有專攻，他們專注濕實驗以及科研想法，在數(shù)據(jù)分析方面略微薄弱，也忽略了其中的細(xì)節(jié)描述。然后我們就干脆直接聯(lián)系了完成芯片服務(wù)的公司，拿到了全部數(shù)據(jù)分析資料。
其中一張圖文并茂的指示引起了我的注意：
圖中強調(diào)了是上調(diào)基因集合，作為背景基因，然后我剎那間相通了為什么他們圖中是50個基因左右的數(shù)量的這兩個通路，因為背景基因由兩萬多個被挑選為上調(diào)基因的4000多個，所以通路基因集必然也會縮水很多。
再次回顧他們的原始圖：
的確，p值和NES值都對，就是基因只有2700左右，很明顯，他們是選擇了上調(diào)探針的4000多個，但是只有其中的2700左右個探針是可以注釋到基因的（關(guān)于這個探針注釋到基因，是另外一個坑，跟本文無關(guān)）
既然知道了他們的策略，那么我就很容易寫代碼重復(fù)這個分析：
rm(list = ls()) ## 魔幻操作，一鍵清空~
options(stringsAsFactors = F)
library(GEOquery)
library(ggplot2)
library(clusterProfiler)
library(GSEABase)
# 這里載入前面下載到的GEO數(shù)據(jù)集，看前面推文即可
load(file = 'step1-output.Rdata')
colnames(dat)
dat[1:4,1:4]
m=dat[,2]+dat[,1]/2
logFC1=log2(dat[,2]/dat[,1]);fivenum(logFC1)
logFC2=dat[,2]-dat[,1];fivenum(logFC2)
fivenum(logFC2)
## 這里我們也挑選基因
table(logFC2>0.7)
logFC2=logFC2[logFC2>0.7]
geneList=scale(logFC2)
names(geneList)=rownames(geneList)
geneList=sort(geneList,decreasing = T)
head(geneList)
# 這里需要注意的是 要 以 entrez ID來進(jìn)行后續(xù)GSEA分析
# 如果是symbol，需要使用其它包，比如 GSEABase
library(org.Mm.eg.db)
s2g=select(org.Mm.eg.db,names(geneList),
'ENTREZID','SYMBOL')
head(s2g)
s2g=s2g[!is.na(s2g$ENTREZID),]
s2g=s2g[s2g$SYMBOL %in% names(geneList),]
geneList=geneList[s2g$SYMBOL]
names(geneList)=s2g$ENTREZID
head(geneList)
tail(geneList)
## 這里比較耗時：
if(F){
go_bp_gsea <- gseGO(geneList = geneList,
OrgDb='org.Mm.eg.db',
ont = 'BP',
nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.9,
verbose = FALSE)
gseaplot(go_bp_gsea, geneSetID = rownames(go_bp_gsea[1,]))
gseaplot(go_bp_gsea, geneSetID = "GO:0006959")
gseaplot(go_bp_gsea, geneSetID = "GO:0030595")
tmp=go_bp_gsea@result
table(tmp$pvalue<0.01)
}
看其中一個結(jié)果，可以看到幾乎再現(xiàn)了作者的結(jié)果，現(xiàn)在的差異就來自于我們選取的基因數(shù)量不一樣，但都是根據(jù)logFC值來挑選了上調(diào)基因后做的GSEA分析，我們需要討論這樣的合理性。
雖然我這邊重復(fù)出來了作者的GSEA顯著結(jié)果，但是并不意味著萬事大吉，我們?nèi)匀灰懻撘幌卵芯空咚械墓臼褂眠@樣的策略來做GSEA分析是否合理。
首先，在GSEA官網(wǎng)文檔，很清楚的描述：
The GSEA algorithm does not filter the expression dataset and does not benefit from your filtering of the expression dataset. During the analysis, genes that are poorly expressed or that have low variance across the dataset populate the middle of the ranked gene list and the use of a weighted statistic ensures that they do not contribute to a positive enrichment score. By removing such genes from your dataset, you may actually reduce the power of the statistic.
也就是說，不建議自行過濾基因，更別說是挑選上調(diào)表達(dá)基因了。
這個公司的策略是挑選上調(diào)基因，再去做gsea看某些通路是否富集，的確是合理性不足！
但是顧博士給我一個鏈接：
Although GSEA does not require that you preprocess the expression dataset, it can be used effectively on preprocessed datasets. For example, Monti et al used a filtered dataset to further analyze genes consistently expressed across two datasets
參考文獻(xiàn)是：“Molecular profiling of diffuse large B-cell lymphoma identifies robust subtypes including one characterized by host inflammatory response” (http://www.bloodjournal.org/cgi/content/full/bloodjournal;105/5/1851).
里面提到了作者收集整理了相關(guān)領(lǐng)域的4個研究的281個基因集，但是他們的GSEA方法，就是挑選基因后的GSEA，Enrichment was assessed by:
(1) ranking the 2118 genes in the top 5% pool with respect to the phenotype “cluster X versus not cluster X”;
(2) locating the represented members of a given gene set within the ranked 2118 genes;
關(guān)于公司這個GSEA分析策略的統(tǒng)計學(xué)合理性，在我們生信技能樹VIP群有著激烈的討論：
藥企東哥發(fā)言：
很多研究領(lǐng)域已經(jīng)有了確定的發(fā)現(xiàn)，這些背景知識可以作為positive control用于判斷以及合理的調(diào)整分析策略。跟研究目的也有關(guān)吧。不一定非得為了“totally unbiased result”閉著眼跑“標(biāo)準(zhǔn)流程”，畢竟大多數(shù)生物信息方法都是不完美的。分析結(jié)果可解釋很重要
其他人發(fā)言：
略
