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

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

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

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

2.GEO主頁下載series_matrix.txt文件

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

那么下面，我們來看看GO中的BP、MF和CC。

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

PS：

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

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