大家好，我是在生信技能書練習(xí)了2月半的學(xué)徒，愛(ài)好是唱，跳，rap和生信。之前在復(fù)現(xiàn)論文的時(shí)候，有遇到過(guò)PPI網(wǎng)絡(luò)，老師說(shuō)可以使用STRING這個(gè)網(wǎng)頁(yè)工具來(lái)完成。這個(gè)PPI網(wǎng)絡(luò)圖不僅可以表達(dá)蛋白之間的相互作用，而且還能把相應(yīng)基因的功能展現(xiàn)出來(lái)。其實(shí)還可以做KEGG/GO注釋，那么R語(yǔ)言的代碼也可以做，他們二者的區(qū)別是什么呢？就是我們接下來(lái)要講的東西啦。

      首先，在生信學(xué)徒培訓(xùn)的前一個(gè)半月里，主要是攻克了R語(yǔ)言，然后做了一些RNA-seq和GEO數(shù)據(jù)的挖掘。這里展示一下GEO數(shù)據(jù)挖掘的流程。

1.GEO主頁(yè)下載原始數(shù)據(jù)（RAW.tar）

      下載下來(lái)是CEL格式，需要自己進(jìn)行預(yù)處理。hgu 95系列和133系列的芯片常用affy包中ReadAffy函數(shù)進(jìn)行讀取，有些平臺(tái)用affy包處理不了，例如[HuGene-1_1-st] Affymetrix Human Gene 1.1 ST Array這個(gè)平臺(tái)，就需要oligo包read.celfiles函數(shù)進(jìn)行讀取。illumina的芯片用lumi包來(lái)處理。之后可以進(jìn)行rma或mas5進(jìn)行歸一化操作。

2.GEO主頁(yè)下載series_matrix.txt文件

3.GEOquery包直接下載表達(dá)量

          getGEO這個(gè)函數(shù)，輸入不同的參數(shù)，下載的東西不一樣。輸入?yún)?shù)是GDS號(hào)時(shí)，下載soft文件，參數(shù)是GPL號(hào)時(shí)，下載芯片設(shè)計(jì)信息，參數(shù)是GSE號(hào)時(shí)，下載series_matrix.txt.gz文件，返回的是ExpressionSet對(duì)象，需要掌握geneNames，sampleNames，pData，exprs等對(duì)ExpressionSet對(duì)象操作的函數(shù)。

從GEO上下載的表達(dá)譜的行名是probe_id探針名，但是不同的平臺(tái)，探針名不同，我們也無(wú)法直觀地知道某個(gè)樣本在某個(gè)探針上的表達(dá)量是那個(gè)基因的表達(dá)量，于是就需要將探針名轉(zhuǎn)換為大家公認(rèn)的NCBI的entrez ID，或者HUGO組織規(guī)定的gene symbol以便于后續(xù)分析。于是，我們要根據(jù)不同的GPL找到該芯片平臺(tái)有對(duì)應(yīng)的bioconductor注釋包來(lái)找到探針與基因的對(duì)應(yīng)關(guān)系，再進(jìn)行轉(zhuǎn)換。這里會(huì)遇到，“一個(gè)探針對(duì)應(yīng)著多個(gè)基因”或者“一個(gè)基因?qū)?yīng)多個(gè)探針”或者“探針沒(méi)有對(duì)應(yīng)基因”的情況，這就需要過(guò)濾整合表達(dá)矩陣，處理方法不盡相同。

        表達(dá)矩陣描述的就是各個(gè)基因在各個(gè)樣本上的表達(dá)量。這講主要是表達(dá)矩陣的可視化。無(wú)論是芯片表達(dá)數(shù)據(jù)或是轉(zhuǎn)錄組高通量測(cè)序數(shù)據(jù)，下載完表達(dá)譜需要根據(jù)生物學(xué)背景驗(yàn)證一下表達(dá)譜是不是正確的。只有確定了所得表達(dá)譜是正確的，之后的差異分析等一系列分析手段才是有意義的。這里提到的方法是看看管家基因的表達(dá)量是不是在表達(dá)譜中處于高表達(dá)。也可以用boxplot看看每個(gè)樣本的表達(dá)量分布圖，看看是否有批次效應(yīng)等等。這里就需要去了解一些畫圖函數(shù)的使用方法。

表達(dá)矩陣的提取（這里的'GSE42872’是個(gè)例子）

分組信息的話，就通過(guò)下面的方法得到。而且，分組信息的個(gè)數(shù)和樣本是一一對(duì)應(yīng)的。

在R中看到的是這樣子的：

差異分析呢，就是把表達(dá)量特別高和表達(dá)量特別低的基因給篩選出來(lái)，因?yàn)槔碚撋希挥羞@種不平凡的基因，才會(huì)對(duì)你想要研究的東西影響最大。提取出來(lái)了之后，用圖形和表格直觀地展示出來(lái)，就是所謂的數(shù)據(jù)可視化。

下面的代碼，就是在R中，設(shè)置條件，篩選出差異基因DEGs（differentially expressed genes）.一般來(lái)說(shuō)，火山圖，MA圖和熱圖都是我們DEGs可視化的選擇。

###不知道怎么作差異分析和可視化，不知道怎么用R。無(wú)所謂，大神已經(jīng)把代碼post到這里https://github.com/jmzeng1314/GEO，操作順序放了在這里。https://www.bilibili.com/video/av25643438

   介紹基因的注釋及富集分析。差異分析通過(guò)自定義的閾值挑選了有統(tǒng)計(jì)學(xué)顯著的基因列表后我們其實(shí)是需要對(duì)它們進(jìn)行注釋才能了解其功能，最常見(jiàn)的就是GO/KEGG數(shù)據(jù)庫(kù)注釋，當(dāng)然也可以使用Reactome和Msigdb數(shù)據(jù)庫(kù)來(lái)進(jìn)行注釋。而最常見(jiàn)的注釋方法就是超幾何分布檢驗(yàn)。

當(dāng)然還有其他的注釋方法。超幾何分布檢驗(yàn)，運(yùn)用到通路的富集概念就是“總共有多少基因（這個(gè)地方值得注意，主流認(rèn)為只考慮那些在KEGG等數(shù)據(jù)庫(kù)注釋的背景基因），你的通路有多少基因，你的通路被抽中了多少基因（在差異基因里面屬于你的通路的基因）?！?/span>目的就是知道，哪些通路中的哪些基因的表達(dá)因?yàn)樗幬锘蛘吣承┎僮鞯淖饔冒l(fā)生了較大的變化，導(dǎo)致通路有較大改變。

       KEGG輸入的基因是EntrezID，在此之前需要進(jìn)行轉(zhuǎn)換。當(dāng)然，上面的ID轉(zhuǎn)換已經(jīng)包括在里面了，其實(shí)蠻多人是會(huì)嫌麻煩漏掉這一步的。

       在Ｒ中如何進(jìn)行注釋，這里就不在多說(shuō)，不知道如何運(yùn)用R或者還沒(méi)有試過(guò)在R中進(jìn)行GO/KEGG注釋的小伙伴們，可以到JM大神的b站觀看視頻。

       我們得到了篩選出來(lái)的DEGs，還可以通過(guò)包來(lái)做ID轉(zhuǎn)換，把symbol轉(zhuǎn)換成ENSEMBL的蛋白ID。但是，之前本人轉(zhuǎn)換過(guò)了，發(fā)現(xiàn)ENSEMBL的prot ID post上去匹配不了，后來(lái)的某天早晨，由于看了cxk的籃球視頻，我直接把symbol list放上了STRING,發(fā)現(xiàn)居然可以識(shí)別，而且自動(dòng)匹配成對(duì)應(yīng)的蛋白ID！

只要把你的基因粘貼到右邊的大方框，下面選好物種就OK了。當(dāng)然，記得左邊選擇第三個(gè)，除非你是有蛋白ID或者是AA序列。

就會(huì)導(dǎo)出一個(gè)PPI圖。有很多圓球和連線。這些又大又圓的球代表的是基因，也可以是蛋白質(zhì)。在圖中，用Node表示。而且那些又細(xì)又長(zhǎng)的是連接Node的線，叫做edge。edge不僅連接node,而且還有表示interaction的功能。

點(diǎn)擊Node，還有有相關(guān)的信息和域的顯示。

不僅如此，下面的Analysis，還有整個(gè)PPI基因/蛋白的GO，KEGG注釋。

在它輸出的文檔里面，前面三個(gè)download都屬于圖形文件，下面的三個(gè)屬于文本文件，可以用來(lái)導(dǎo)入cytoscape.可以從下面的表中看到，PPI各個(gè)node的關(guān)系都已經(jīng)列好，還對(duì)應(yīng)出每個(gè)蛋白ID與注釋信息，連它們間的score都有了。這樣，就可以基本得到了PPI和比較全面的interaction信息了。

好了，下面的就是從STRING上面下載的5個(gè)download文件。可見(jiàn)，下面5個(gè)文本文件分別為：string_interaction.tsv(以tab分割)；XML 總結(jié)；網(wǎng)絡(luò)坐標(biāo)；蛋白序列和蛋白注釋。應(yīng)該都可以用excel打開(kāi)。

在string_interaction中，有15列，上圖為前面6列。第一列為每一個(gè)node的基因名，3,5分別是它們對(duì)應(yīng)的內(nèi)部ID和蛋白ID（這里的蛋白ID還在前面加了物種編號(hào)）。2則是和之前那個(gè)node有關(guān)系的另一個(gè)node,4,6也分別是node2的對(duì)應(yīng)內(nèi)部ID和蛋白ID。

后面9列，分別為染色體上臨近點(diǎn)，基因融合，系統(tǒng)發(fā)育，同源性，共表達(dá)，實(shí)驗(yàn)性相互關(guān)系，數(shù)據(jù)庫(kù)注釋，自動(dòng)文本挖掘，綜合評(píng)分。

network_coordinate記錄的是Node名稱，坐標(biāo)，顏色和注釋。

protein_sequences記錄的是氨基酸的序列。

protein_annotations及記錄了基因名，蛋白ID和結(jié)構(gòu)域的信息來(lái)源。但是由于格式太大，用excel不能完全打開(kāi)。最后一個(gè)XML，是以psimi格式制備的，因此不適宜用excel打開(kāi)。不然看起來(lái)就像cxk打籃球一樣。

       而在Ｒ中，也同樣可以對(duì)基因進(jìn)行KEGG/GO注釋。那到底哪個(gè)更方便，更可信呢。

　　在Ｒ中如何進(jìn)行注釋，這里就不在多說(shuō)，不知道如何運(yùn)用R或者還沒(méi)有試過(guò)在R中進(jìn)行GO/KEGG注釋的小伙伴們，可以到JM大神的b站觀看視頻。        

       那我們就分別來(lái)對(duì)比下同樣的基因，在STRING和R得到的KEGG/GO注釋有啥區(qū)別。這里主要是比較STRING和R中的KEGG/GO的background gene庫(kù)的大小。如果數(shù)據(jù)庫(kù)太久或比較小，有很多基因就沒(méi)有被收錄進(jìn)去，這樣有可能我們的目的基因就不會(huì)被注釋到。（GO注釋包括BP，CC和MF）。

       基因名如何導(dǎo)入，和網(wǎng)站如何使用，JM大神已經(jīng)在視頻里有詳細(xì)說(shuō)明了，而我們就在PPI圖下面的exports下載相應(yīng)的文件就好了。

#注意：在名為'GeneRation’和'BgRatio’兩列的數(shù)據(jù)里，我們只看分子。

在KEGG注釋方面，我們可以看到，各自的區(qū)別不大。

那么下面，我們來(lái)看看GO中的BP、MF和CC。

在GO注釋方面，同樣識(shí)別的基因和background區(qū)別也不大，所以在KEGG/GO功能注釋中，兩種方法大家都可以放心使用。

PS：

雖說(shuō)大多說(shuō)情況如此，既然可以在STRING這種online tools中做出來(lái)的東西，為何我要在R中敲代碼來(lái)實(shí)現(xiàn)呢。

然后，我就發(fā)現(xiàn)了某些功能，STRING是很籠統(tǒng)的歸為一類，而R中，則會(huì)進(jìn)行比較細(xì)致的分類。在這，R中，可以通過(guò)p,q值進(jìn)行cutoff ,而在STRING中，只能通過(guò)調(diào)節(jié)interaction score來(lái)cutoff了。所以STRING幾秒鐘的便捷，和R中細(xì)膩還是有一點(diǎn)區(qū)別的，看大家所需吧。畢竟魚與熊掌不可兼得。（但AJ和鋼絲球可以）