你不得不知道的GEO數(shù)據(jù)庫芯片數(shù)據(jù)分析方法
大家好妇穴，我是阿琛。今天來給大家介紹一個新的專題內(nèi)容---GEO數(shù)據(jù)庫的使用方法隶债。烹飪需要食材腾它，分析需要數(shù)據(jù)。數(shù)據(jù)出發(fā)死讹，整個研究的第一步就是數(shù)據(jù)的下載瞒滴。對于大部分的研究者而言，拿公開的高通量數(shù)據(jù)赞警，進行二次分析妓忍，是最佳的選擇途徑。
作為與TCGA數(shù)據(jù)庫齊名的一個大型數(shù)據(jù)庫愧旦，GEO數(shù)據(jù)庫包羅萬象世剖，對于每個領(lǐng)域的科研工作者很有幫助。GEO數(shù)據(jù)庫是一個儲存芯片笤虫、二代測序以及其他高通量測序數(shù)據(jù)的一個數(shù)據(jù)庫旁瘫。利用這個數(shù)據(jù)庫祖凫，我們可以檢索到其他一些人上傳的一些實驗測序數(shù)據(jù)。
下面酬凳，我們來看一下如何使用GEO數(shù)據(jù)庫中的芯片數(shù)據(jù)進行后續(xù)分析惠况。
GEO數(shù)據(jù)庫的簡介
GEO數(shù)據(jù)庫，GENE EXPRESSION OMNIBUS宁仔，是由美國國立生物技術(shù)信息中心NCBI創(chuàng)建并維護的基因表達數(shù)據(jù)庫（官網(wǎng)：https://www.ncbi.nlm.nih.gov/geo/）稠屠。
它最初創(chuàng)建于2000年，主要用于收錄各國研究機構(gòu)提交的高通量基因表達數(shù)據(jù)翎苫，也就是說只要是在目前已經(jīng)發(fā)表的絕大部分論文中完箩，其涉及到的基因表達檢測的數(shù)據(jù)，包括芯片數(shù)據(jù)拉队，二代測序結(jié)果弊知，以及其他形式的高通量檢測結(jié)果，都可以通過這個數(shù)據(jù)庫中找到粱快。
首先秩彤，我們進入GEO數(shù)據(jù)庫中，根據(jù)GSE編號事哭，查看一下該數(shù)據(jù)的一些相關(guān)信息漫雷。在搜索欄中輸入編號，“GSE39582”鳍咱，然后點擊“Search”按鈕降盹，進行檢索。
當然谤辜，熟悉了以后蓄坏，我們也可以直接輸入網(wǎng)址進行快速檢索，https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39582丑念；對于檢索頁面網(wǎng)址而言涡戳，其前面是一樣的，唯一的區(qū)別是acc=后面的GSE編號脯倚，修改編號渔彰，可以快速進入對應(yīng)的結(jié)果頁面，這樣也可以在一定程度上減少由于網(wǎng)速等原因所帶來的阻礙推正。
在結(jié)果頁面中恍涂，我們可以初步看一下該結(jié)果對應(yīng)的發(fā)布時間沪摄，題目康震，物種等信息。同時霎桅，Experiment type中Expression profiling by array表示該結(jié)果是通過芯片獲得的表達譜内贮；Overall design中簡單介紹了整個研究的設(shè)計方案和分組信息产园。
GEO數(shù)據(jù)的下載
當獲得了芯片的GSE編號后汞斧，我們接下來需要將其對應(yīng)的數(shù)據(jù)進行下載，從而根據(jù)自己的需要進一步分析什燕。關(guān)于數(shù)據(jù)的下載粘勒，我們這里主要介紹三種不同的方法。
方法一：下載芯片的原始數(shù)據(jù)
在檢索頁面中屎即，一路下拉庙睡；在Supplementary file中點擊Download中的custom，展開原始數(shù)據(jù)對應(yīng)列表技俐；點擊“Select all“乘陪，然后點擊Download，即可將所有樣本的原始數(shù)據(jù)RAW Data文件下載雕擂；
雖然這是最直接的方法啡邑，但是RAW Data文件相對較大，對下載的網(wǎng)速要求相對較高井赌，而且不同的芯片來源谤逼，有不同的處理方法，甚至有些芯片沒有處理方法仇穗，因為其是對應(yīng)定制的流部。所以，一般情況下纹坐，不推薦大家下載原始數(shù)據(jù)枝冀。
方法二：下載表達矩陣（series matrix）
在Download family中點擊Serier Matrix File(s)，進入下載頁面耘子；
待下載完成后果漾，可以直接使用read.table()函數(shù)讀取進來。
可以看到拴还，在芯片中跨晴，包含了54675個基因探針欧聘，586個患者片林。
方法三：使用R的GEOquery包里面的getGEO()函數(shù)直接讀取進來（推薦）
當然，考慮到網(wǎng)速問題怀骤，我們可以對參數(shù)進行設(shè)置费封，選擇不下載平臺的注釋文件，因為一般來講注釋文件是相對比較大的蒋伦。
如果把之前下載的series Matrix文件放在當前目錄下弓摘，getGEO()函數(shù)會直接檢測到該文件，并進而直接將其進行讀群劢臁韧献；
我們可以直接查看一下下載結(jié)果gset的變量類型末患。
可以看到，變量gset是一個列表的形式锤窑。
為什么是list格式呢璧针？因為一個GEO芯片項目，是可以對應(yīng)多個芯片平臺的渊啰，那么每個平臺的數(shù)據(jù)結(jié)果會對應(yīng)list里面的一個元素探橱。
既然是列表，自然可以提取其中的第一個元素出來查看一下绘证∷砀啵可以看到，其中展示了包含了6個樣本嚷那，33297個特征胞枕，以及相關(guān)的臨床信息，PMID號魏宽，以及注釋平臺信息曲稼。
提取表達和臨床信息
3.1 通過pData函數(shù)獲取分組信息
通過pData()函數(shù)，即可提取表達數(shù)據(jù)中的臨床信息湖员；同時贫悄，點擊Environment中的pdata查看，我們可以查看里面的相關(guān)信息娘摔≌梗可以看到，其中凳寺，非腫瘤的樣品19例鸭津。
因此，根據(jù)臨床信息肠缨，我們可以對樣品進行分組逆趋，分為腫瘤組和正常組；
最終晒奕，得到正常組19例樣品闻书，腫瘤組566例樣品。
3.2 通過exprs()函數(shù)獲取表達矩陣并校正
整理完臨床信息后脑慧，我們需要提取對應(yīng)的表達數(shù)據(jù)魄眉。對于表達數(shù)據(jù)，除了下載Series Matrix后直接使用read.table()函數(shù)進行讀取外闷袒，我們也可以直接從GEOquery下載得到的變量gset中進行提取坑律。
使用exprs()函數(shù)可以從gset[[1]]提取表達信息；同時囊骤，我們可以使用boxplot()函數(shù)先簡單看一下整體樣本的表達情況晃择。
由于每一次技術(shù)重復(fù)的時候冀值，都會有誤差，芯片的原始數(shù)據(jù)是由儀器讀取的宫屠，不同的讀取時間池摧，或者掃描儀光線的強弱都會導(dǎo)致同一類型的樣本出現(xiàn)誤差。正式分析前激况，我們需要對其進行人工校正一下作彤。這里我們用limma包內(nèi)置的一個函數(shù)，
normalizeBetweenArrays()函數(shù)乌逐。
可以看到竭讳，經(jīng)過校正，整個表達水平基本趨于一致浙踢。
此外绢慢，使用range()函數(shù)查看一下表達數(shù)據(jù)exp的取值范圍；一般而言洛波，范圍在20以內(nèi)的表達值基本已經(jīng)經(jīng)過了log對數(shù)轉(zhuǎn)換胰舆。
ID變換
整理好了表達矩陣以后，我們需要將探針的id轉(zhuǎn)換成為基因的Gene symbol蹬挤。對于探針id的轉(zhuǎn)換過程缚窿，目前主要是通過R包來進行轉(zhuǎn)換。接下來焰扳，我們來看一下如何進行芯片探針id的轉(zhuǎn)換過程倦零。
方法一：使用R包轉(zhuǎn)換
隨著芯片平臺的普遍使用，其基因的注釋信息也被整理成了不同的R包吨悍；因此扫茅，通常情況下我們使用R包來注釋。不同的平臺育瓜，對應(yīng)著不同的R包葫隙。首先，我們來看一下這個數(shù)據(jù)集使用的平臺類型躏仇。
通過提取列表gset[[1]]中的注釋信息恋脚，可以看到，該芯片使用的是我們最常見的平臺钙态，GPL570慧起。
對于GPL570，其對應(yīng)的R包是hgu133plus2.db包册倒；查找顯示，其儲存在Bioconductor中磺送，下載并進行安裝R包驻子。
首先灿意，我們來看看，在hgu133plus2.db包中崇呵，包含了哪些信息缤剧；
可以看到，除了Symbol信息外域慷，在其中還包含了Ensemble id荒辕，Entrez id等信息，可以需要進行提取犹褒。
提取其中的Symbol信息抵窒，可以看到，最終獲得了probe id和Gene symbol的對應(yīng)信息叠骑。
其中李皇，經(jīng)過去重復(fù)，一共存在20174個不同的Gene symbol宙枷，且部分基因存在多條探針的對應(yīng)關(guān)系掉房。
接下來，我們需要將其進行一一對應(yīng)匹配慰丛。
經(jīng)過id匹配卓囚，并去重復(fù)，最終得到了20174個基因的表達結(jié)果诅病；
同時捍岳，我們可以查看一下前3行前3列的表達情況。
當然睬隶，除了R包注釋外锣夹，還有其他的注釋方法，比如使用網(wǎng)頁下載的soft文件進行注釋苏潜，或者有些特殊的芯片內(nèi)容银萍，需要自己手工比對注釋。
方法二：使用soft文件注釋
方法三：手工注釋
PCA分析
表達矩陣到此基本整理完成恤左。接下來贴唇，在正式的差異分析之前，我們首先可以做一個PCA分析飞袋，整體水平查看正常和腫瘤兩組樣品直接是否存在顯著的差異戳气。
結(jié)果顯示：在該芯片中，癌和癌旁組織的表達水平存在一定的差異巧鸭。
差異分析
對于芯片數(shù)據(jù)的差異分析瓶您，我們一般使用limma包來進行。關(guān)于差異分析的輸入文件，主要是兩個呀袱，第一是整理好的表達矩陣贸毕，其中行名為基因名，列名為樣本名夜赵；第二是分組信息（group list）明棍。
最終，使用topTable()函數(shù)提取所有基因的差異分析寇僧。
可以看到摊腋，在結(jié)果表格中，包含了6塊的內(nèi)容嘁傀，包括我們常見的logFC值兴蒸，以及P.Value，adj.P.Val等等心包。
接下來类咧，我們需要根據(jù)設(shè)定的閾值，｜logFC｜>1.5和P值
可視化分析
1蟹腾、火山圖
對于差異分析結(jié)果痕惋，火山圖和熱圖是兩種最常見的展示方式。首先娃殖，我們來看一下火山圖的繪制方法值戳。對于火山圖的繪制，大家可以參考之前的推文（生信最重要的圖之一炉爆，十分鐘幫你搞定堕虹！建議收藏！７沂住）赴捞。
方法一：基于ggpubr包繪制火山圖
結(jié)果顯示：
接下來，我們可以進一步給火山圖添加標簽郁稍，把顯著上調(diào)和顯著下調(diào)基因中前5名的基因名進行標注赦政；
結(jié)果顯示：
方法二：基于ggplot2包繪制火山圖
結(jié)果顯示：
當然，我們也可以對其進行添加標簽的操作耀怜；
結(jié)果顯示：
2恢着、熱圖
提取差異表達基因的表達情況；
結(jié)果顯示：
到此财破，GEO數(shù)據(jù)庫芯片數(shù)據(jù)的下載掰派，probe id的轉(zhuǎn)換，差異分析已經(jīng)基本完成了左痢，整個文章中最難的80%內(nèi)容也已經(jīng)基本解決靡羡。接下來系洛，就是針對這些差異基因的常規(guī)分析，包括GO分析亿眠，KEGG分析碎罚，GSEA分析磅废，蛋白-蛋白互作網(wǎng)絡(luò)（Protein-protein interaction纳像，PPI）。
