春節(jié)后的第一次寫學習筆記震庭，還是要提醒大家沒事別聚集鼻由、別聚集、別聚集橱赠！重要的事情說三遍尤仍！
前幾天看到公眾號里寫了GSEA的文章，雖然我之前跟著很多網上的教程走過一遍GSEA分析狭姨，但是對于GSEA的結果圖并不知道怎么看宰啦，所以也想學習一下這方面的知識。
在跟著教程走代碼之前饼拍，想先了解一下什么是GSEA赡模，以及這個東西是干嘛用的。
（一）什么是GSEA?
GSEA: Gene Set Enrichment Analysis
翻譯應該不難吧： 基因集富集分析
（二）上面說的集是什么意思惕耕？
參考文章：GSEA入門------原理
在以前的實驗中發(fā)表的數據或表達譜上共表達的基因信息數據集合（瞧這專業(yè)的說法纺裁，我看了好幾遍才看明白）。說的通俗一點：就是某一個通路（相關的所有基因的總和）。
（三）GSEA分析的是什么東東欺缘？
我們在利用DESeq2栋豫，或者edgeR，或者limma進行差異分析之后谚殊，會得到一個列表丧鸯，里面有你的所有樣品的差異基因。我們對這個列表里的基因表達水平進行一個排序嫩絮，從高到低丛肢。GSEA分析就是確定一個基因集S（或者一個通路）的成員是否傾向于出現在列表L的頂部(或底部)，如果聚集在頂部/底部剿干，該基因集與表型分類區(qū)分相關蜂怎。GSEA是一種無閾值方法，可根據其差異表達等級或其他分數對所有基因進行分析，無需事先進行基因過濾（學習基因通路富集分析軟件GSEA）。
（四）GSEA原理是什么鸟顺？
實在是不想看一堆公式。幽歼。。所以就把教程里（劉小澤學習GSEA）的一段話再說的明白點：
假設：某一個通路的全部基因在排序后的差異基因列表中隨機分布谬盐。
但是：如果我們看到它們“意外地”出現在基因列表的某一端（從圖上看是在某一側形成一個峰）甸私，那么就可以計算顯著性來看看富集程度如何。如果富集結果顯著飞傀，那么就拒絕原假設皇型，認為這個通路的基因在我們的基因列表中富集，并且看到富集分數（是個發(fā)paper的契機爸觥）犀被。
（五）上面提到的富集分數是什么意思？
富集分數：Enrichment Score外冀，簡稱ES。
從我們排序的差異基因列表里掀泳，從高到低一個一個來看：遇到一個基因雪隧，屬于我們要的基因集S時，ES就會加分员舵；如果遇到的這個基因不屬于我們要的基因集脑沿，ES就會減分。來一張圖马僻，更直觀一些：
我們通常看到的GSEA的圖是長這樣的：
（六）實際操作：方法一：R語言
這里我用的是DESeq2得到的差異基因列表瞭郑，這個列表是怎么得來的辜御，請參考：TCGA的差異基因分析。如果你有自己的差異基因列表更好屈张。這個差異基因的列表可以是DESeq2擒权，也可以是edgeR或者limma分析得到的。
#GSEA analysis steps
#得到差異分析結果
> geneList <- DESeq2_DEG$log2FoldChange
#調用需要的R包
> library(stringr)
> library(clusterProfiler)
> library(org.Hs.eg.db)
先看一下我的這個DESeq2差異基因列表阁谆，可以看到基因名稱都是Ensemble ID：
我們需要把基因名稱先轉換成Entrez ID:
#如果是Ensemble ID碳抄，并且如果還帶著版本號，需要先去除版本號
> names(geneList) <- str_split(rownames(DESeq2_DEG) ,
pattern = '\\.',simplify = T)[,1]
#下面用clusterProfiler包進行基因名稱的轉換：
> geneList_tr <- bitr(names(geneList),
fromType = "ENSEMBL",
toType = c("ENTREZID","SYMBOL"),
OrgDb = org.Hs.eg.db)
> new_list <- data.frame(ENSEMBL=names(geneList),
logFC = as.numeric(geneList))
> new_list <- merge(new_list, geneList_tr, by = "ENSEMBL")
> geneList <- new_list$logFC
> names(geneList) <- geneList_tr$ENTREZID
# 最后從大到小排序场绿，得到一個字符串
> geneList <- sort(geneList,decreasing = T)
#進行GO分析
> go_result <- gseGO(geneList = geneList,
ont = "BP",
OrgDb = org.Hs.eg.db,
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 1)
#將其中的基因名變成symbol ID
> go_result <- setReadable(go_result, OrgDb = org.Hs.eg.db)
#還可以直接點擊查看纳鼎，只需要轉成數據框
> go_result_df <- as.data.frame(go_result)
#如果對結果的某個通路感興趣，可以作圖
> gseaplot(go_result, geneSetID = c("GO:0010638"))
#進行KEGG分析
> kk <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 1,
verbose = FALSE)
> go_result_kk <- as.data.frame(kk)
> gseaplot(kk, geneSetID = "hsa05202")
（七）實際操作：GSEA軟件(這個方法麻煩得很啊裳凸，還是R語言方便快捷贱鄙，但做出來的圖漂亮，有閑工夫的可以試試）
（1）軟件下載
GSEA的官網：https://www.gsea-msigdb.org/gsea/index.jsp
這里以windows系統(tǒng)為例姨谷。安裝軟件就不說了逗宁，so easy。
運行一下吧梦湘，界面長這樣：
（2）準備需要的文件
先別著急用瞎颗，為什么我說這個軟件很煩躁呢，因為它不是啥格式的文件都識別的捌议。
首先哼拔，是你的基因列表。這里的基因列表請注意瓣颅！不是差異基因列表倦逐，是你所有的基因！９埂Ｃ世选（參考文章：Question: GSEA error ）可以先打開看一下我們的表（這個表是怎么來的，可以參照利用R包TCGAbiolinks進行各種數據下載去下載RNA-seq的FPKM結果來進行練習）：
然而這樣的表是不能夠被GSEA軟件識別的贫贝。根據官網的說明：Dataformats:
The TXT format is a tab delimited file format that describes an expression dataset. It is organized as follows:
The first line contains the labels Name and Description followed by the identifiers for each sample in the dataset. NOTE: The Description column is intended to be optional, but there is currently a bug such that it is treated as required. We hope to fix this in a future release. If you have no descriptions available, a value of NA will suffice.
注意注意注意１獭ｒ让铡！這里的示范列表里比我們多出了一列DESCRIPTION崇堵。這一列是必須有的型诚。你要自己添加進去。如果沒有什么可描述的筑辨，就寫上na俺驶。如果這里你不注意，之后軟件會一直報錯棍辕，報到你懷疑人生Ｄ合帧！楚昭！
你的基因列表還需要符合下面的格式：
第一行: Name(tab)Description(tab)(sample 1 name)(tab)(sample 2 name) (tab) ... (sample N name)
NOTE：必須是制表符分隔栖袋！別的符號統(tǒng)統(tǒng)不可以！
那么怎么讓它變成制表符分隔呢抚太？請看：
在excel里點擊“導出”塘幅，再點擊“更改文件類型”：
選擇“文本文件（制表符分隔）”，最后點擊“另存為”就好了尿贫。要確保是.txt后綴的文件格式电媳。
接下來還需要一個文件，就是你的分組信息庆亡。在前面的文章里匾乓，提到了如何把TCGA數據庫里的樣品分組（TCGA的差異基因分析）。那么這個分組信息的文件格式也是有講究的又谋，先看一下官網上給出的示范：
那么根據官網的例子丰介，我也整理了一下我的分組信息：
NOTE:這里需要注意的是哮幢，這個文件你需要保存成.cls后綴的文件带膀。
（3）萬事俱備，開始分析
有了必須的文件橙垢，就可以進行分析了垛叨。
首先，上傳你的文件（method1選擇你的文件）：
如果你的文件格式沒問題，這里會提示你成功上傳了2個文件裁奇。
然后點擊左邊一欄的"Run GSEA":
然后就是一系列的選擇參數了：
然后選擇Gene sets database桐猬，點右邊的...會出來一個選擇框：
是不是覺得好多？選擇困難癥要犯了框喳。课幕。。我好南五垮。乍惊。。和我一樣的童鞋可以看這個官網：https://www.gsea-msigdb.org/gsea/msigdb放仗，這里對每一個database都進行了說明：
但看完感覺還是不知道選哪個怎么辦润绎？我問了一下實驗室的Steven，他說比如說你的樣品是control和knockout诞挨，你想看看這兩者的區(qū)別莉撇，那么就選擇H，因為這里面是一些經典的生物過程和通路的基因惶傻，比較完全棍郎。或者是每一種里的"all"银室。但是你如果做免疫涂佃，需要專業(yè)性更強的database励翼，那么就選C7」架總之汽抚，根據你的實驗需求來選。這里我選擇H伯病。
接下來造烁，“number of permutation",置換次數。1000是默認值午笛。內存夠大CPU夠強惭蟋，就選1000，這里選100季研。越大越好敞葛，越大結果越精確。電腦不行的請選擇100与涡，當然了值越大惹谐，出來的結果會越精確。參考：學習基因通路富集分析軟件GSEA
更多的置換次數需要更長的計算時間驼卖。 為了計算每個gene set的FDR Q值氨肌，通過置換每個基因組中的基因并重新計算隨機組的P值來隨機化數據集,此參數指定完成此隨機化的次數。執(zhí)行的隨機化越多酌畜，FDR Q值估計就越精確（達到極限怎囚，因為最終FDR Q值將穩(wěn)定在實際值）。
然后桥胞，phenotype label:
我選擇的是tumor vs normal恳守。
下面的collapse/Remap to gene symbols的選擇，如果你的列表里基因名稱和我的一樣贩虾，都是ENSG開頭的催烘，那么你需要選擇collapse；如果你的基因名稱是gene symbol缎罢，那就選No_collapse伊群。
Permutaion type的選擇：樣品數量少，選gene_set策精；樣品數量多舰始，選phenotype。
這一部分的最后一個參數：chip platform咽袜。點開也有很多的選擇丸卷，有的文章里說這一步不用選，但是對我來說询刹，如果我空著它及老，運行時會一直報錯抽莱！這里我選擇的是human_ensembl_gene_ID這一項：
下面是第二部分的參數范抓，我是都用的默認值骄恶，只是修改一下輸出文件的儲存路徑：
最后點擊最下方的"Run"，如果你的文件格式沒有錯匕垫，選擇的參數也對應你的文件僧鲁，比如你的基因名稱和collapse/no_collapse是否對應。那么在左下角就會出現running的字樣：
等一段時間（取決于你的電腦性能象泵，以及你選擇的參數）寞秃，運行成功會顯示：
NOTE:如果點擊運行，左下角出現紅色字Error偶惠，那就是你的文件格式有問題春寿，要么就是你collapse選擇錯了，要么就是你空著chip platform了忽孽。
（4）生成的文件怎么看
運行后绑改，生成了一堆的文件：
這些文件那些比較重要呢？首先你可以找一下這個文件：
點開它兄一，看一下summary：
然后我們點開看看詳細情況厘线，會告訴你這24個gene set都是什么：
點擊details,會出現對于這個gene set的分析結果:
重點來了：這個圖怎么看？
先看官網放的圖出革，有點糊：
根據官網的圖造壮，我得到的這個圖分為三部分：
（1）第一部分：綠色的折線圖。就是我們上面說到的加分減分得到的曲線骂束。我的圖這里的峰是向下的耳璧，它的富集分數ES是-0.733。
（2）第二部分：中間黑色的像條形碼一樣的圖展箱。是這個gene set里所有的基因旨枯，這個gene set里有200個基因，所以這些黑線理論上有200條析藕≌偻ⅲ可以看到大部分的黑線都集中在右邊，也就是說這個gene set里的大部分基因在tumor vs normal里是下調的账胧。
（3）第三部分：灰色部分竞慢。是所有基因的rank 值分布圖。
你還會在下面看到這個gene set里的所有基因的情況：
關于上面的有些參數治泥，我就不詳細說明了筹煮，請參考文章，這兩篇文章里有相當詳細的說明居夹，我就不做搬運工了败潦，copy也沒啥意思：
1.學習基因通路富集分析軟件GSEA
2.GSEA教程
總結
學習生信一定要自己親自實踐一遍本冲！光看是沒有用的！因為網上的教程很多劫扒，說的步驟也都是五花八門檬洞，不同文章里的不同點就是你在操作的時候容易出錯的地方。我也跟著不同的教程走流程沟饥，發(fā)現不是每一個教程都完全正確添怔。在操作過程中會遇到各種各樣的問題，報錯報到頭大贤旷。有時一個報錯要各種的搜索文章广料，各種百度google，折騰上好幾天也是常有的事幼驶。但是每次靠著自己解決問題的那種成就感是無法言喻的艾杏。所以要相信自己！只要沒氣到把電腦砸了盅藻，就一定會找到解決問題的方法购桑！
加油！
