In molecular biology, STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) is a biological database and web resource of known and predicted protein–protein interactions.(from Wkkipedia)
大家好,我是在生信技能書練習了2月半的學徒,愛好是唱,跳,rap和生信。之前在復現(xiàn)論文的時候,有遇到過PPI網(wǎng)絡,老師說可以使用STRING這個網(wǎng)頁工具來完成。這個PPI網(wǎng)絡圖不僅可以表達蛋白之間的相互作用,而且還能把相應基因的功能展現(xiàn)出來。其實還可以做KEGG/GO注釋,那么R語言的代碼也可以做,他們二者的區(qū)別是什么呢?就是我們接下來要講的東西啦。
首先,在生信學徒培訓的前一個半月里,主要是攻克了R語言,然后做了一些RNA-seq和GEO數(shù)據(jù)的挖掘。這里展示一下GEO數(shù)據(jù)挖掘的流程。
下載數(shù)據(jù)(這里提供三種方式)
1.GEO主頁下載原始數(shù)據(jù)(RAW.tar)
下載下來是CEL格式,需要自己進行預處理。hgu 95系列和133系列的芯片常用affy包中ReadAffy函數(shù)進行讀取,有些平臺用affy包處理不了,例如[HuGene-1_1-st] Affymetrix Human Gene 1.1 ST Array這個平臺,就需要oligo包read.celfiles函數(shù)進行讀取。illumina的芯片用lumi包來處理。之后可以進行rma或mas5進行歸一化操作。
2.GEO主頁下載series_matrix.txt文件
exprSet = read.table('GSE42872_series_matrix.txt.gz', sep = '\t', fill = T, comment.char = "!", header = T)3.GEOquery包直接下載表達量
library(GEOquery)getGEO這個函數(shù),輸入不同的參數(shù),下載的東西不一樣。輸入?yún)?shù)是GDS號時,下載soft文件,參數(shù)是GPL號時,下載芯片設計信息,參數(shù)是GSE號時,下載series_matrix.txt.gz文件,返回的是ExpressionSet對象,需要掌握geneNames,sampleNames,pData,exprs等對ExpressionSet對象操作的函數(shù)。
ID轉(zhuǎn)換、表達矩陣
從GEO上下載的表達譜的行名是probe_id探針名,但是不同的平臺,探針名不同,我們也無法直觀地知道某個樣本在某個探針上的表達量是那個基因的表達量,于是就需要將探針名轉(zhuǎn)換為大家公認的NCBI的entrez ID,或者HUGO組織規(guī)定的gene symbol以便于后續(xù)分析。于是,我們要根據(jù)不同的GPL找到該芯片平臺有對應的bioconductor注釋包來找到探針與基因的對應關系,再進行轉(zhuǎn)換。這里會遇到,“一個探針對應著多個基因”或者“一個基因?qū)鄠€探針”或者“探針沒有對應基因”的情況,這就需要過濾整合表達矩陣,處理方法不盡相同。
表達矩陣描述的就是各個基因在各個樣本上的表達量。這講主要是表達矩陣的可視化。無論是芯片表達數(shù)據(jù)或是轉(zhuǎn)錄組高通量測序數(shù)據(jù),下載完表達譜需要根據(jù)生物學背景驗證一下表達譜是不是正確的。只有確定了所得表達譜是正確的,之后的差異分析等一系列分析手段才是有意義的。這里提到的方法是看看管家基因的表達量是不是在表達譜中處于高表達。也可以用boxplot看看每個樣本的表達量分布圖,看看是否有批次效應等等。這里就需要去了解一些畫圖函數(shù)的使用方法。
表達矩陣的提取(這里的'GSE42872’是個例子)
library(GEOquery)分組信息的話,就通過下面的方法得到。而且,分組信息的個數(shù)和樣本是一一對應的。
library(GEOquery)在R中看到的是這樣子的:
差異分析及可視化
差異分析呢,就是把表達量特別高和表達量特別低的基因給篩選出來,因為理論上,只有這種不平凡的基因,才會對你想要研究的東西影響最大。提取出來了之后,用圖形和表格直觀地展示出來,就是所謂的數(shù)據(jù)可視化。
下面的代碼,就是在R中,設置條件,篩選出差異基因DEGs(differentially expressed genes).一般來說,火山圖,MA圖和熱圖都是我們DEGs可視化的選擇。
###不知道怎么作差異分析和可視化,不知道怎么用R。無所謂,大神已經(jīng)把代碼post到這里https://github.com/jmzeng1314/GEO,操作順序放了在這里。https://www.bilibili.com/video/av25643438
富集分析KEGG、GO注釋
介紹基因的注釋及富集分析。差異分析通過自定義的閾值挑選了有統(tǒng)計學顯著的基因列表后我們其實是需要對它們進行注釋才能了解其功能,最常見的就是GO/KEGG數(shù)據(jù)庫注釋,當然也可以使用Reactome和Msigdb數(shù)據(jù)庫來進行注釋。而最常見的注釋方法就是超幾何分布檢驗。
當然還有其他的注釋方法。超幾何分布檢驗,運用到通路的富集概念就是“總共有多少基因(這個地方值得注意,主流認為只考慮那些在KEGG等數(shù)據(jù)庫注釋的背景基因),你的通路有多少基因,你的通路被抽中了多少基因(在差異基因里面屬于你的通路的基因)?!?/span>目的就是知道,哪些通路中的哪些基因的表達因為藥物或者某些操作的作用發(fā)生了較大的變化,導致通路有較大改變。
KEGG輸入的基因是EntrezID,在此之前需要進行轉(zhuǎn)換。當然,上面的ID轉(zhuǎn)換已經(jīng)包括在里面了,其實蠻多人是會嫌麻煩漏掉這一步的。
在R中如何進行注釋,這里就不在多說,不知道如何運用R或者還沒有試過在R中進行GO/KEGG注釋的小伙伴們,可以到JM大神的b站觀看視頻。
https://www.bilibili.com/video/av25643438https://www.bilibili.com/video/av26731585
看完教學視頻,下面的圖表唾手可得?。?!
STRING的基操和文件下載
我們得到了篩選出來的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)居然可以識別,而且自動匹配成對應的蛋白ID!
只要把你的基因粘貼到右邊的大方框,下面選好物種就OK了。當然,記得左邊選擇第三個,除非你是有蛋白ID或者是AA序列。
就會導出一個PPI圖。有很多圓球和連線。這些又大又圓的球代表的是基因,也可以是蛋白質(zhì)。在圖中,用Node表示。而且那些又細又長的是連接Node的線,叫做edge。edge不僅連接node,而且還有表示interaction的功能。
點擊Node,還有有相關的信息和域的顯示。
不僅如此,下面的Analysis,還有整個PPI基因/蛋白的GO,KEGG注釋。
在它輸出的文檔里面,前面三個download都屬于圖形文件,下面的三個屬于文本文件,可以用來導入cytoscape.可以從下面的表中看到,PPI各個node的關系都已經(jīng)列好,還對應出每個蛋白ID與注釋信息,連它們間的score都有了。這樣,就可以基本得到了PPI和比較全面的interaction信息了。
好了,下面的就是從STRING上面下載的5個download文件??梢姡旅?個文本文件分別為:string_interaction.tsv(以tab分割);XML 總結(jié);網(wǎng)絡坐標;蛋白序列和蛋白注釋。應該都可以用excel打開。
在string_interaction中,有15列,上圖為前面6列。第一列為每一個node的基因名,3,5分別是它們對應的內(nèi)部ID和蛋白ID(這里的蛋白ID還在前面加了物種編號)。2則是和之前那個node有關系的另一個node,4,6也分別是node2的對應內(nèi)部ID和蛋白ID。
后面9列,分別為染色體上臨近點,基因融合,系統(tǒng)發(fā)育,同源性,共表達,實驗性相互關系,數(shù)據(jù)庫注釋,自動文本挖掘,綜合評分。
network_coordinate記錄的是Node名稱,坐標,顏色和注釋。
protein_sequences記錄的是氨基酸的序列。
protein_annotations及記錄了基因名,蛋白ID和結(jié)構(gòu)域的信息來源。但是由于格式太大,用excel不能完全打開。最后一個XML,是以psimi格式制備的,因此不適宜用excel打開。不然看起來就像cxk打籃球一樣。
STRING與R的background gene區(qū)別
而在R中,也同樣可以對基因進行KEGG/GO注釋。那到底哪個更方便,更可信呢。
在R中如何進行注釋,這里就不在多說,不知道如何運用R或者還沒有試過在R中進行GO/KEGG注釋的小伙伴們,可以到JM大神的b站觀看視頻。
那我們就分別來對比下同樣的基因,在STRING和R得到的KEGG/GO注釋有啥區(qū)別。這里主要是比較STRING和R中的KEGG/GO的background gene庫的大小。如果數(shù)據(jù)庫太久或比較小,有很多基因就沒有被收錄進去,這樣有可能我們的目的基因就不會被注釋到。(GO注釋包括BP,CC和MF)。
基因名如何導入,和網(wǎng)站如何使用,JM大神已經(jīng)在視頻里有詳細說明了,而我們就在PPI圖下面的exports下載相應的文件就好了。
#注意:在名為'GeneRation’和'BgRatio’兩列的數(shù)據(jù)里,我們只看分子。
在KEGG注釋方面,我們可以看到,各自的區(qū)別不大。
那么下面,我們來看看GO中的BP、MF和CC。
在GO注釋方面,同樣識別的基因和background區(qū)別也不大,所以在KEGG/GO功能注釋中,兩種方法大家都可以放心使用。
PS:
雖說大多說情況如此,既然可以在STRING這種online tools中做出來的東西,為何我要在R中敲代碼來實現(xiàn)呢。
然后,我就發(fā)現(xiàn)了某些功能,STRING是很籠統(tǒng)的歸為一類,而R中,則會進行比較細致的分類。在這,R中,可以通過p,q值進行cutoff ,而在STRING中,只能通過調(diào)節(jié)interaction score來cutoff了。所以STRING幾秒鐘的便捷,和R中細膩還是有一點區(qū)別的,看大家所需吧。畢竟魚與熊掌不可兼得。(但AJ和鋼絲球可以)