|
你不得不知道的GEO數(shù)據(jù)庫(kù)芯片數(shù)據(jù)分析方法 大家好，我是阿琛。今天來(lái)給大家介紹一個(gè)新的專題內(nèi)容---GEO數(shù)據(jù)庫(kù)的使用方法。烹飪需要食材，分析需要數(shù)據(jù)。數(shù)據(jù)出發(fā)，整個(gè)研究的第一步就是數(shù)據(jù)的下載。對(duì)于大部分的研究者而言，拿公開(kāi)的高通量數(shù)據(jù)，進(jìn)行二次分析，是最佳的選擇途徑。 作為與TCGA數(shù)據(jù)庫(kù)齊名的一個(gè)大型數(shù)據(jù)庫(kù)，GEO數(shù)據(jù)庫(kù)包羅萬(wàn)象，對(duì)于每個(gè)領(lǐng)域的科研工作者很有幫助。GEO數(shù)據(jù)庫(kù)是一個(gè)儲(chǔ)存芯片、二代測(cè)序以及其他高通量測(cè)序數(shù)據(jù)的一個(gè)數(shù)據(jù)庫(kù)。利用這個(gè)數(shù)據(jù)庫(kù)，我們可以檢索到其他一些人上傳的一些實(shí)驗(yàn)測(cè)序數(shù)據(jù)。 下面，我們來(lái)看一下如何使用GEO數(shù)據(jù)庫(kù)中的芯片數(shù)據(jù)進(jìn)行后續(xù)分析。 GEO數(shù)據(jù)庫(kù)的簡(jiǎn)介 GEO數(shù)據(jù)庫(kù)，GENE EXPRESSION OMNIBUS，是由美國(guó)國(guó)立生物技術(shù)信息中心NCBI創(chuàng)建并維護(hù)的基因表達(dá)數(shù)據(jù)庫(kù)（官網(wǎng)：https://www.ncbi.nlm./geo/）。 它最初創(chuàng)建于2000年，主要用于收錄各國(guó)研究機(jī)構(gòu)提交的高通量基因表達(dá)數(shù)據(jù)，也就是說(shuō)只要是在目前已經(jīng)發(fā)表的絕大部分論文中，其涉及到的基因表達(dá)檢測(cè)的數(shù)據(jù)，包括芯片數(shù)據(jù)，二代測(cè)序結(jié)果，以及其他形式的高通量檢測(cè)結(jié)果，都可以通過(guò)這個(gè)數(shù)據(jù)庫(kù)中找到。  image
首先，我們進(jìn)入GEO數(shù)據(jù)庫(kù)中，根據(jù)GSE編號(hào)，查看一下該數(shù)據(jù)的一些相關(guān)信息。在搜索欄中輸入編號(hào)，“GSE39582”，然后點(diǎn)擊“Search”按鈕，進(jìn)行檢索。 當(dāng)然，熟悉了以后，我們也可以直接輸入網(wǎng)址進(jìn)行快速檢索，https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE39582；對(duì)于檢索頁(yè)面網(wǎng)址而言，其前面是一樣的，唯一的區(qū)別是acc=后面的GSE編號(hào)，修改編號(hào)，可以快速進(jìn)入對(duì)應(yīng)的結(jié)果頁(yè)面，這樣也可以在一定程度上減少由于網(wǎng)速等原因所帶來(lái)的阻礙。  image
在結(jié)果頁(yè)面中，我們可以初步看一下該結(jié)果對(duì)應(yīng)的發(fā)布時(shí)間，題目，物種等信息。同時(shí)，Experiment type中Expression profiling by array表示該結(jié)果是通過(guò)芯片獲得的表達(dá)譜；Overall design中簡(jiǎn)單介紹了整個(gè)研究的設(shè)計(jì)方案和分組信息。  image
GEO數(shù)據(jù)的下載 當(dāng)獲得了芯片的GSE編號(hào)后，我們接下來(lái)需要將其對(duì)應(yīng)的數(shù)據(jù)進(jìn)行下載，從而根據(jù)自己的需要進(jìn)一步分析。關(guān)于數(shù)據(jù)的下載，我們這里主要介紹三種不同的方法。 方法一：下載芯片的原始數(shù)據(jù) 在檢索頁(yè)面中，一路下拉；在Supplementary file中點(diǎn)擊Download中的custom，展開(kāi)原始數(shù)據(jù)對(duì)應(yīng)列表；點(diǎn)擊“Select all“，然后點(diǎn)擊Download，即可將所有樣本的原始數(shù)據(jù)RAW Data文件下載；  image
雖然這是最直接的方法，但是RAW Data文件相對(duì)較大，對(duì)下載的網(wǎng)速要求相對(duì)較高，而且不同的芯片來(lái)源，有不同的處理方法，甚至有些芯片沒(méi)有處理方法，因?yàn)槠涫菍?duì)應(yīng)定制的。所以，一般情況下，不推薦大家下載原始數(shù)據(jù)。 方法二：下載表達(dá)矩陣（series matrix） 在Download family中點(diǎn)擊Serier Matrix File(s)，進(jìn)入下載頁(yè)面； 待下載完成后，可以直接使用read.table()函數(shù)讀取進(jìn)來(lái)。 可以看到，在芯片中，包含了54675個(gè)基因探針，586個(gè)患者。 方法三：使用R的GEOquery包里面的getGEO()函數(shù)直接讀取進(jìn)來(lái)（推薦） 當(dāng)然，考慮到網(wǎng)速問(wèn)題，我們可以對(duì)參數(shù)進(jìn)行設(shè)置，選擇不下載平臺(tái)的注釋文件，因?yàn)橐话銇?lái)講注釋文件是相對(duì)比較大的。 如果把之前下載的series Matrix文件放在當(dāng)前目錄下，getGEO()函數(shù)會(huì)直接檢測(cè)到該文件，并進(jìn)而直接將其進(jìn)行讀取；  image
我們可以直接查看一下下載結(jié)果gset的變量類型。 可以看到，變量gset是一個(gè)列表的形式。 為什么是list格式呢？因?yàn)橐粋€(gè)GEO芯片項(xiàng)目，是可以對(duì)應(yīng)多個(gè)芯片平臺(tái)的，那么每個(gè)平臺(tái)的數(shù)據(jù)結(jié)果會(huì)對(duì)應(yīng)list里面的一個(gè)元素。 既然是列表，自然可以提取其中的第一個(gè)元素出來(lái)查看一下?？梢钥吹剑渲姓故玖税?個(gè)樣本，33297個(gè)特征，以及相關(guān)的臨床信息，PMID號(hào)，以及注釋平臺(tái)信息。  image
提取表達(dá)和臨床信息 3.1 通過(guò)pData函數(shù)獲取分組信息 通過(guò)pData()函數(shù)，即可提取表達(dá)數(shù)據(jù)中的臨床信息；同時(shí)，點(diǎn)擊Environment中的pdata查看，我們可以查看里面的相關(guān)信息。可以看到，其中，非腫瘤的樣品19例。 因此，根據(jù)臨床信息，我們可以對(duì)樣品進(jìn)行分組，分為腫瘤組和正常組； 最終，得到正常組19例樣品，腫瘤組566例樣品。 3.2 通過(guò)exprs()函數(shù)獲取表達(dá)矩陣并校正 整理完臨床信息后，我們需要提取對(duì)應(yīng)的表達(dá)數(shù)據(jù)。對(duì)于表達(dá)數(shù)據(jù)，除了下載Series Matrix后直接使用read.table()函數(shù)進(jìn)行讀取外，我們也可以直接從GEOquery下載得到的變量gset中進(jìn)行提取。 使用exprs()函數(shù)可以從gset[[1]]提取表達(dá)信息；同時(shí)，我們可以使用boxplot()函數(shù)先簡(jiǎn)單看一下整體樣本的表達(dá)情況。 由于每一次技術(shù)重復(fù)的時(shí)候，都會(huì)有誤差，芯片的原始數(shù)據(jù)是由儀器讀取的，不同的讀取時(shí)間，或者掃描儀光線的強(qiáng)弱都會(huì)導(dǎo)致同一類型的樣本出現(xiàn)誤差。正式分析前，我們需要對(duì)其進(jìn)行人工校正一下。這里我們用limma包內(nèi)置的一個(gè)函數(shù)， normalizeBetweenArrays()函數(shù)。 可以看到，經(jīng)過(guò)校正，整個(gè)表達(dá)水平基本趨于一致。 此外，使用range()函數(shù)查看一下表達(dá)數(shù)據(jù)exp的取值范圍；一般而言，范圍在20以內(nèi)的表達(dá)值基本已經(jīng)經(jīng)過(guò)了log對(duì)數(shù)轉(zhuǎn)換。 ID變換 整理好了表達(dá)矩陣以后，我們需要將探針的id轉(zhuǎn)換成為基因的Gene symbol。對(duì)于探針id的轉(zhuǎn)換過(guò)程，目前主要是通過(guò)R包來(lái)進(jìn)行轉(zhuǎn)換。接下來(lái)，我們來(lái)看一下如何進(jìn)行芯片探針id的轉(zhuǎn)換過(guò)程。 方法一：使用R包轉(zhuǎn)換 隨著芯片平臺(tái)的普遍使用，其基因的注釋信息也被整理成了不同的R包；因此，通常情況下我們使用R包來(lái)注釋。不同的平臺(tái)，對(duì)應(yīng)著不同的R包。首先，我們來(lái)看一下這個(gè)數(shù)據(jù)集使用的平臺(tái)類型。 通過(guò)提取列表gset[[1]]中的注釋信息，可以看到，該芯片使用的是我們最常見(jiàn)的平臺(tái)，GPL570。  image
對(duì)于GPL570，其對(duì)應(yīng)的R包是hgu133plus2.db包；查找顯示，其儲(chǔ)存在Bioconductor中，下載并進(jìn)行安裝R包。 首先，我們來(lái)看看，在hgu133plus2.db包中，包含了哪些信息； 可以看到，除了Symbol信息外，在其中還包含了Ensemble id，Entrez id等信息，可以需要進(jìn)行提取。  image
提取其中的Symbol信息，可以看到，最終獲得了probe id和Gene symbol的對(duì)應(yīng)信息。  image
其中，經(jīng)過(guò)去重復(fù)，一共存在20174個(gè)不同的Gene symbol，且部分基因存在多條探針的對(duì)應(yīng)關(guān)系。 接下來(lái)，我們需要將其進(jìn)行一一對(duì)應(yīng)匹配。 經(jīng)過(guò)id匹配，并去重復(fù)，最終得到了20174個(gè)基因的表達(dá)結(jié)果；  image
同時(shí)，我們可以查看一下前3行前3列的表達(dá)情況。 當(dāng)然，除了R包注釋外，還有其他的注釋方法，比如使用網(wǎng)頁(yè)下載的soft文件進(jìn)行注釋，或者有些特殊的芯片內(nèi)容，需要自己手工比對(duì)注釋。 方法二：使用soft文件注釋 方法三：手工注釋 PCA分析 表達(dá)矩陣到此基本整理完成。接下來(lái)，在正式的差異分析之前，我們首先可以做一個(gè)PCA分析，整體水平查看正常和腫瘤兩組樣品直接是否存在顯著的差異。 結(jié)果顯示：在該芯片中，癌和癌旁組織的表達(dá)水平存在一定的差異。  image
差異分析 對(duì)于芯片數(shù)據(jù)的差異分析，我們一般使用limma包來(lái)進(jìn)行。關(guān)于差異分析的輸入文件，主要是兩個(gè)，第一是整理好的表達(dá)矩陣，其中行名為基因名，列名為樣本名；第二是分組信息（group list）。 最終，使用topTable()函數(shù)提取所有基因的差異分析。 可以看到，在結(jié)果表格中，包含了6塊的內(nèi)容，包括我們常見(jiàn)的logFC值，以及P.Value，adj.P.Val等等。 接下來(lái)，我們需要根據(jù)設(shè)定的閾值，｜logFC｜>1.5和P值 可視化分析 1、火山圖 對(duì)于差異分析結(jié)果，火山圖和熱圖是兩種最常見(jiàn)的展示方式。首先，我們來(lái)看一下火山圖的繪制方法。對(duì)于火山圖的繪制，大家可以參考之前的推文（生信最重要的圖之一，十分鐘幫你搞定！建議收藏！?。?/p> 方法一：基于ggpubr包繪制火山圖 結(jié)果顯示： image
接下來(lái)，我們可以進(jìn)一步給火山圖添加標(biāo)簽，把顯著上調(diào)和顯著下調(diào)基因中前5名的基因名進(jìn)行標(biāo)注； 結(jié)果顯示： image
方法二：基于ggplot2包繪制火山圖 結(jié)果顯示： image
當(dāng)然，我們也可以對(duì)其進(jìn)行添加標(biāo)簽的操作； 結(jié)果顯示： image
2、熱圖 提取差異表達(dá)基因的表達(dá)情況； 結(jié)果顯示： image
到此，GEO數(shù)據(jù)庫(kù)芯片數(shù)據(jù)的下載，probe id的轉(zhuǎn)換，差異分析已經(jīng)基本完成了，整個(gè)文章中最難的80%內(nèi)容也已經(jīng)基本解決。接下來(lái)，就是針對(duì)這些差異基因的常規(guī)分析，包括GO分析，KEGG分析，GSEA分析，蛋白-蛋白互作網(wǎng)絡(luò)（Protein-protein interaction，PPI）。 |
|
|
來(lái)自： imtravelinghah > 《待分類》
