|
Pathview是一個用于整合表達譜數(shù)據(jù)并用于可視化KEGG通路的一個R包,其會先下載KEGG官網(wǎng)上的通路圖,然后整合輸入數(shù)據(jù)對通路圖進行再次渲染,從而對KEGG通路圖進行一定程度上的個性化處理,并且豐富其信息展示。(KEGG在線數(shù)據(jù)庫使用攻略)
Pathview的安裝一種方法是通過Bioconductor安裝,需要Bioconductor版本3.9,R的版本3.6 (推薦)。 if (!requireNamespace('BiocManager', quietly=TRUE)) install.packages('BiocManager') BiocManager::install('pathview')
另一種方法是去網(wǎng)上下載包的壓縮包,https:///packages/release/bioc/html/pathview.html。 下載完壓縮包之后,進入Rstudio,選擇Tools——Install Packages——Browse,找到下載壓縮包的位置,安裝即可。 

什么是KEGG pathway?KEGG (Kyoto Encyclopedia of Genes and Genomes)是系統(tǒng)分析基因功能、基因組信息的數(shù)據(jù)庫。它有助于研究者把基因及表達信息作為一個整體網(wǎng)絡(luò)進行研究。KEGG將基因組信息和基因功能信息有機地結(jié)合起來,通過對細(xì)胞內(nèi)已知生物學(xué)過程進行計算機化處理和將現(xiàn)有的基因功能解釋標(biāo)準(zhǔn)化,對基因的功能進行系統(tǒng)化的分析。KEGG的另一個用途是將基因組中的一系列基因用一個細(xì)胞內(nèi)的分子相互作用網(wǎng)絡(luò)連接起來,如一個通路或是一個復(fù)合物,通過它們來展現(xiàn)更高一級的生物學(xué)功能。
為什么要用KEGG的代謝通路?KEGG提供的整合代謝途徑(pathway)查詢十分出色,包括碳水化合物、核苷、氨基酸等的代謝及有機物的生物降解,不僅提供了所有可能的代謝途徑,而且對催化各步反應(yīng)的酶進行了全面的注解,包含有氨基酸序列、PDB庫的鏈接等等。(沒錢買KEGG怎么辦?REACTOME開源通路更強大) KEGG是進行生物體內(nèi)代謝分析、代謝網(wǎng)絡(luò)研究的強有力工具。與其他數(shù)據(jù)庫相比,KEGG 的一個顯著特點就是具有強大的圖形功能,它利用圖形而不是繁縟的文字來介紹眾多的代謝途徑以及各途徑之間的關(guān)系。
在KEGG中有兩種代謝圖。 參考代謝通路圖 reference pathway,是根據(jù)已有的知識繪制的具有一般參考意義的代謝圖;這種圖上所有小框都是無色的,不會有綠色的小框,并且都可以點擊查看更詳細(xì)的信息; 特定物種的代謝圖 species-specific pathway,會用綠色來標(biāo)出這個物種所有的基因或酶,只有這些綠色的框點擊以后才會給出更詳細(xì)的信息。
這兩種圖很好區(qū)分,reference pathway 在KEGG中的名字是以map開頭的,比如map00010,就是糖酵解途徑的參考圖;而特定物種的代謝通路圖開頭三個字符不是map而是種屬英文單詞的縮寫(一個屬的首字母+2個種的首字母)比如酵母的糖酵解通路圖,就是sce00010,大腸桿菌的糖酵解通路圖就應(yīng)該是eco00010吧。 對差異基因進行pathway分析 (代謝通路),就是把基因表達變化映射到具體的代謝網(wǎng)絡(luò)中,可以研究某個實驗條件下顯著改變的代謝通路,在機制研究中顯得尤為重要。 研究pathway的原因是:生物學(xué)問題中設(shè)定一個“蝴蝶效應(yīng)”假設(shè),1個Pathway上游基因的改變,會導(dǎo)致下游相關(guān)基因改變,從而改變通路中大量基因的表達。 準(zhǔn)備開始Pathview主要用于可視化pathway圖上的數(shù)據(jù)。pathview可以生成KEGG視圖和Graphviz視圖,前者將用戶數(shù)據(jù)呈現(xiàn)在原生KEGG pathway圖上,更自然,更易于閱讀。后者使用Graphviz引擎對pathway圖進行布局,可以更好地控制節(jié)點或邊緣屬性和pathway拓?fù)洹?/p>
首先我們在R里面調(diào)用該包,使用該包自帶的數(shù)據(jù)集。這是一個乳腺癌數(shù)據(jù)集,可以查看下演示數(shù)據(jù)是什么格式的。 列名是每個樣本名,行名是每個基因的entrez id。自己準(zhǔn)備的數(shù)據(jù)要符合這個格式,因為entrez id是行名字,而entrez id都是數(shù)字,讀取需設(shè)置check.names=F。其它類型的ID也支持,需要做一些參數(shù)設(shè)置或轉(zhuǎn)換,具體文后有介紹。 library(pathview) data(gse16873.d) head(gse16873.d)
# 讀取自己的文件可以使用類似下面的命令 # gse16873.d <- read.table('filename', sep='\t', row.names=1, check.names=F, header=T)

先看下Pathview最常見的一種用法:將某個樣本的表達量(讀入的數(shù)據(jù)需要是歸一化后的表達量);其實也可以任何列數(shù)據(jù),不僅僅是表達量數(shù)據(jù),也可以是fold change, p-value, 組蛋白修飾數(shù)據(jù)等,人為指定的數(shù)值型數(shù)據(jù)也行 (關(guān)鍵是要懂需要展示什么數(shù)據(jù)、說明什么問題,原理最重要,就像GSEA基因集富集分析也是一樣);最后以color bar的形式在KEGG通路圖上的對應(yīng)節(jié)點 (一定注意節(jié)點名字的匹配)展示; 如下例子所示,我們通過指定gene.data和pathway.id來觀察單個樣本在典型信號通路細(xì)胞周期上的表達變化。該基因芯片是在人體組織上進行的,因此species=“hsa”。 p <- pathview(gene.data = gse16873.d[, 1], pathway.id = '04110', species = 'hsa', out.suffix = 'gse16873', kegg.native = T)

如何獲取自己關(guān)注的通路的ID呢? 如下動圖,可以得到Cell cyle的ID04110。 
該例子中的圖只有一個單層,在原始圖層修改節(jié)點顏色,保留原始KEGG節(jié)點標(biāo)簽 (節(jié)點名)。這樣輸出的文件大小與原始的KEGG PNG文件一樣小,但是計算時間相對較長。如果我們想要一個快速的視圖,并且不介意將輸出文件大小,我們可以通過same.layer = F使用兩層圖。通過這種方式,節(jié)點顏色和標(biāo)簽被添加到原始KEGG的額外圖層上。原來的KEGG基因標(biāo)簽(或EC編號)被替換為官方基因符號。 p <- pathview(gene.data = gse16873.d[, 1], pathway.id = '04110', species = 'hsa', out.suffix = 'gse16873.2layer', kegg.native = T, same.layer = F)

上面的兩個例子中,我們查看了KEGG pathway圖的數(shù)據(jù),在KEGG圖上我們可以得到所有注釋和meta-data,因此數(shù)據(jù)更具可讀性和解釋性。然而輸出的是PNG格式的柵格圖像。我們也可以使用Graphviz engine重新繪制pathway圖來查看數(shù)據(jù),這樣我們對節(jié)點和邊緣屬性能有更多的控制,更重要的是可以保存為PDF格式的矢量圖像。 p <- pathview(gene.data = gse16873.d[, 1], pathway.id = '04110', species = 'hsa', out.suffix = 'gse16873', kegg.native = F, sign.pos='bottomleft')

該圖的主圖和圖例都在一個圖層或者說一個頁面中,我們只列出KEGG邊的類型,忽略圖例中的節(jié)點類型,以節(jié)省空間。如果我們想要完整的圖例,我們可以使用兩個層來創(chuàng)建Graphviz視圖: 第1頁是主圖,第2頁是圖例。
注意,對于Graphviz視圖 (PDF文件),“層”的概念與KEGG視圖 (PNG文件)略有不同。在這兩種情況下,我們都為兩層圖設(shè)置參數(shù)same.layer=F。 p <- pathview(gene.data = gse16873.d[, 1], pathway.id = '04110', species = 'hsa', out.suffix = 'gse16873.2layer', kegg.native = F, sign.pos='bottomleft', same.layer = F)


圖形布局樣式調(diào)整在Graphviz視圖中,我們對圖形布局有更多的控制,比如可以將節(jié)點組拆分為獨立的節(jié)點,甚至可以將多基因節(jié)點擴展為單個基因。分裂的節(jié)點或擴展的基因可能從未分裂的組或未擴展的節(jié)點繼承邊。這樣我們就可以得到一個基因/蛋白-基因/蛋白相互作用網(wǎng)絡(luò)。 可以更好地查看網(wǎng)絡(luò)特性(模塊化等)和基因方面(而不是節(jié)點方面)的數(shù)據(jù)。注意在KEGG視圖中,一個基因節(jié)點可能代表多個功能相似或重復(fù)的基因/蛋白。成員基因的數(shù)量從1到幾十不等。 為了更好的清晰度和可讀性,一般將它們作為路徑圖上的單個節(jié)點放在一起。因此,默認(rèn)情況下,我們不分割節(jié)點并單獨標(biāo)記每個成員基因。但是,我們可以通過總結(jié)基因數(shù)據(jù)來可視化節(jié)點數(shù)據(jù),用戶可以使用node.sum參數(shù)指定摘要方法。 p <- pathview(gene.data = gse16873.d[, 1], pathway.id = '04110', species = 'hsa', out.suffix = 'gse16873.split', kegg.native = F, sign.pos='bottomleft', split.group = T)

下面的圖更復(fù)雜了,對簡單通路適用,復(fù)雜通路就頭禿了?。?/p> p <- pathview(gene.data = gse16873.d[, 1], pathway.id = '04110', species = 'hsa', out.suffix = 'gse16873.split.expanded', kegg.native = F, sign.pos='bottomleft', split.group = T, expand.node = T)

數(shù)據(jù)整合Pathview為數(shù)據(jù)集成提供了強大的支持。它可以用來整合、分析和可視化各種各樣的生物數(shù)據(jù):基因表達、遺傳關(guān)聯(lián)、代謝產(chǎn)物、基因組數(shù)據(jù)、文獻和其他可映射到通路的數(shù)據(jù)類型。當(dāng)數(shù)據(jù)映射到KEGG ortholog pathways時,它可以直接用于宏基因組、微生物組或未知物種的數(shù)據(jù)。
化合物和基因集同時繪制在通路上在上面的例子中,我們查看了具有典型的信號通路的基因數(shù)據(jù)。有時候我們也想研究代謝通路。除了基因節(jié)點外,這些通路還有復(fù)合節(jié)點。因此,我們可以將基因數(shù)據(jù)和化合物數(shù)據(jù)與代謝途徑進行整合或可視化。這里的基因數(shù)據(jù)是一個廣泛的概念,包括基因、轉(zhuǎn)錄本、蛋白質(zhì)、酶及其表達、修飾和任何可測量的屬性。化合物數(shù)據(jù)也是如此,包括代謝物、藥物、它們的測量值和屬性。這里我們?nèi)匀皇褂萌橄侔┪㈥嚵袛?shù)據(jù)集作為基因數(shù)據(jù)。然后生成模擬的化合物或代謝組數(shù)據(jù),并加載適當(dāng)?shù)幕衔颕D類型(具有足夠數(shù)量的惟一條目)進行演示。 sim.cpd.data=sim.mol.data(mol.type='cpd', nmol=3000) data(cpd.simtypes) head(sim.cpd.data)
我們查看該數(shù)據(jù)部分內(nèi)容如下: C00232 C01881 C02424 C07418 C13756 C07378 -1.09759772 0.12891537 2.07851027 0.93299245 -0.00786048 -0.09023300
我們生成了一個包含基因數(shù)據(jù)和化合物數(shù)據(jù)的KEGG視圖。pathview生成的代謝通路圖與原始KEGG圖相同,只是為了更好地查看顏色,將復(fù)合節(jié)點放大。 p <- pathview(gene.data = gse16873.d[, 1], cpd.data = sim.cpd.data, pathway.id ='00640', species = 'hsa', out.suffix = 'gse16873.cpd', keys.align = 'y', kegg.native = T, key.pos = 'topright')

我們還生成了相同pathway和數(shù)據(jù)的Graphviz視圖。Graphviz視圖更好地顯示了層次結(jié)構(gòu)。對于代謝通路,解析xml文件中的反應(yīng)條目,并將其轉(zhuǎn)換為基因和復(fù)合節(jié)點之間的關(guān)系。對復(fù)合節(jié)點使用省略號。標(biāo)簽是從CHEMBL數(shù)據(jù)庫中檢索到的標(biāo)準(zhǔn)化合物名稱 (KEGG在pathway數(shù)據(jù)庫文件中沒有提供它)。化學(xué)名稱是長字符串,我們需要對它們進行換行,以使其符合圖上指定的寬度。 p <- pathview(gene.data = gse16873.d[, 1], cpd.data = sim.cpd.data, pathway.id ='00640', species = 'hsa', out.suffix = 'gse16873.cpd', keys.align = 'y', kegg.native = F, key.pos = 'bottomright', sign.pos ='topright', cpd.lab.offset =-0.8)

多狀態(tài)或樣本同時或分開繪制在前面的所有示例中,我們查看了單個樣本數(shù)據(jù),這些數(shù)據(jù)要么是向量,要么是單列矩陣。Pathview還可以處理多個樣本數(shù)據(jù),用于為每個樣本生成圖形。從1.1.6版開始,Pathview就可以整合并繪制多狀態(tài)或樣本到一個圖中。 set.seed(10) sim.cpd.data2 = matrix(sample(sim.cpd.data, 18000, replace = T), ncol = 6) rownames(sim.cpd.data2) = names(sim.cpd.data) colnames(sim.cpd.data2) = paste('exp', 1:6, sep = '') head(sim.cpd.data2, 3)

在下面的例子中,gene.data有三個樣本,而cpd.data有兩個。我們可以把所有這些樣品放在一張圖里,繪制KEGG視圖或Graphviz視圖。在這些圖中,我們看到基因節(jié)點和化合物節(jié)點被分割成多個對應(yīng)于不同狀態(tài)或樣本的片段 (注意顏色塊,之前是一個節(jié)點一個顏色,現(xiàn)在一個節(jié)點是有多個顏色,每個對應(yīng)一個樣本,基因有3個,化合物有2個,注意下面代碼中的1:3和1:2)。如果兩種數(shù)據(jù)類型中的樣本實際配對,需要選擇匹配數(shù)據(jù),即gene.data和cpd.data的第一列來自同一個實驗/樣本,等等。 # KEGG view p <- pathview(gene.data = gse16873.d[, 1:3], cpd.data = sim.cpd.data2[, 1:2], pathway.id = '00640', species = 'hsa', out.suffix = 'gse16873.cpd.3-2s', keys.align = 'y', kegg.native = T, match.data = F, multi.state = T, same.layer = T)

查看下繪圖的數(shù)據(jù) head(p$plot.data.cpd)

基因表達和化合物都展示前3個樣品,一一對應(yīng)。 #KEGG view with data match p <- pathview(gene.data = gse16873.d[, 1:3], cpd.data = sim.cpd.data2[, 1:2], pathway.id = '00640', species = 'hsa', out.suffix = 'gse16873.cpd.3-2s.match', keys.align = 'y', kegg.native = T, match.data = T, multi.state = T, same.layer = T)

同樣,也可以使用graphviz view。 #graphviz view p <- pathview(gene.data = gse16873.d[, 1:3], cpd.data = sim.cpd.data2[, 1:2], pathway.id ='00640', species = 'hsa', out.suffix = 'gse16873.cpd.3-2s', keys.align = 'y', kegg.native = F, match.data = F, multi.state = T, same.layer = T, key.pos = 'bottomright', sign.pos = 'topright')

同樣,我們可以選擇分別繪制樣本,即每個圖形一個樣本。請注意,在這種情況下,必須匹配gene.data和cpd.data中的樣本,以便將其分配給同一圖表。因此,參數(shù)match.data在這里并不是很有效。(圖3就沒有化合物的映射了) #plot samples/states separately p <- pathview(gene.data = gse16873.d[, 1:3], cpd.data = sim.cpd.data2[, 1:2], pathway.id = '00640', species = 'hsa', out.suffix = 'gse16873.cpd.3-2s', keys.align = 'y', kegg.native = T, match.data = F, multi.state = F, same.layer = T)



如上所述,同一層上的KEGG視圖可能需要一些時間。同樣,如果我們不介意丟失原始的KEGG基因標(biāo)簽(或EC編號),我們可以選擇使用兩層的KEGG視圖來加快這個過程。 p <- pathview(gene.data = gse16873.d[, 1:3], cpd.data = sim.cpd.data2[, 1:2], pathway.id = '00640', species = 'hsa', out.suffix = 'gse16873.cpd.3-2s.2layer', keys.align = 'y', kegg.native = T, match.data = F, multi.state = T, same.layer = F)

離散數(shù)據(jù)標(biāo)記上下調(diào)或是否存在到目前為止,我們一直在處理連續(xù)數(shù)據(jù)。但我們也經(jīng)常處理離散數(shù)據(jù)。例如,我們根據(jù)一些統(tǒng)計數(shù)據(jù)(p值、折疊變化等)選擇顯著基因或化合物列表。輸入數(shù)據(jù)可以命名為兩個層次的向量,1或0(顯著或不顯著),也可以是一個更短的顯著基因/化合物名稱列表。在接下來的兩個例子中,我們只使gene.data和cpd.data或gene.data離散。 require(org.Hs.eg.db) gse16873.t <- apply(gse16873.d, 1, function(x) t.test(x, alternative = 'two.sided')$p.value) sel.genes <- names(gse16873.t)[gse16873.t < 0.1] sel.cpds <- names(sim.cpd.data)[abs(sim.cpd.data) > 0.5]
我們分別查看下sel.genes和sel.cpds的數(shù)據(jù)結(jié)構(gòu): 選擇高亮的基因 head(sel.genes)
[1] '10000' '10003' [3] '10007' '100128414' [5] '100129271' '100129762'
選擇高亮的化合物 head(sel.cpds)
[1] 'C00232' 'C02424' 'C07418' [4] 'C06633' 'C00251' 'C01059'
p <- pathview(gene.data = sel.genes, cpd.data = sel.cpds, pathway.id ='00640', species = 'hsa', out.suffix = 'sel.genes.sel.cpd', keys.align = 'y', kegg.native = T, key.pos = 'topright', limit = list(gene = 1, cpd = 1), bins = list(gene = 1, cpd = 1), na.col = 'gray', discrete = list(gene = T, cpd = T))

p <- pathview(gene.data = sel.genes, cpd.data = sim.cpd.data, pathway.id = '00640', species = 'hsa', out.suffix = 'sel.genes.cpd', keys.align = 'y', kegg.native = T, key.pos = 'topright', limit = list(gene = 1, cpd = 1), bins = list(gene = 1, cpd = 10), na.col = 'gray', discrete = list(gene = T, cpd = F))

不同來源的ID的轉(zhuǎn)換和映射pathview的一個顯著特點是它強大的ID映射能力。集成的Mapper模塊將10多種基因或蛋白ID、20多種化合物或代謝物ID映射到標(biāo)準(zhǔn)KEGG基因或化合物ID。換句話說,使用這些不同ID類型命名的用戶數(shù)據(jù)可以精確地映射到目標(biāo)KEGG路徑。Pathview適用于大約4800個物種的路徑,物種可以以多種格式指定:KEGG代碼、科學(xué)名稱或常用名稱。
cpd.cas <- sim.mol.data(mol.type = 'cpd', id.type = cpd.simtypes[2], nmol = 10000) gene.ensprot <- sim.mol.data(mol.type = 'gene', id.type = gene.idtype.list[4], nmol = 50000)
同樣查看cpd.simtypes、cpd.cas和gene.ensport的數(shù)據(jù)結(jié)構(gòu): head(cpd.simtypes)
[1] 'Beilstein Registry Number' [2] 'CAS Registry Number' [3] 'ChEMBL COMPOUND' [4] 'KEGG COMPOUND accession' [5] 'KEGG DRUG accession' [6] 'Patent accession'
head(cpd.cas)
是一個named vactor,128470-16-6、465-42-9、25812-30-0 為化合物名字,下面的數(shù)字0.9594129是含量。(不知道為什么會模擬出負(fù)值?可能是歸一化(scale)后的數(shù)據(jù)。具體見R語言 - 熱圖美化)。 128470-16-6 465-42-9 0.9594129 0.6882760 25812-30-0 4790-08-3 -1.5775357 1.2579969 95789-30-3 138-14-7 -1.9112655 -0.1616219
head(gene.ensprot)
是一個named vactor,ENSP00000414006是基因名字,-0.1609490基因表達量 ENSP00000414006 ENSP00000403857 -0.1609490 1.5800692 ENSP00000444905 ENSP00000434312 0.2108686 -0.6047296 ENSP00000493916 ENSP00000338796 -0.9955757 0.5357721
注意參數(shù)中的gene.idtype和cpd.idtype,用來指定基因和化合物的ID類型。 p <- pathview(gene.data = gene.ensprot, cpd.data = cpd.cas, gene.idtype = gene.idtype.list[4], cpd.idtype = cpd.simtypes[2], pathway.id = '00640', species = 'hsa', same.layer = T, out.suffix = 'gene.ensprot.cpd.cas', keys.align = 'y', kegg.native = T, key.pos ='bottomright', sign.pos = 'topright', limit = list(gene = 3, cpd = 3), bins = list(gene = 6, cpd = 6))

對于不在自動映射列表中的外部ID,我們可以使用mol.sum函數(shù) (也是Mapper模塊的一部分)顯式地進行ID和數(shù)據(jù)映射。這里我們需要提供id.map (外部ID和KEGG標(biāo)準(zhǔn)ID之間的映射矩陣)。我們使用ID映射函數(shù)id2eg和cpdidmap等來得到id.map矩陣。注意,這些ID映射函數(shù)可以獨立于pathview 主函數(shù)使用。 id.map.cas <- cpdidmap(in.ids = names(cpd.cas), in.type = cpd.simtypes[2], out.type = 'KEGG COMPOUND accession') cpd.kc <- mol.sum(mol.data = cpd.cas, id.map = id.map.cas) id.map.ensprot <- id2eg(ids = names(gene.ensprot), category = gene.idtype.list[4], org = 'Hs') gene.entrez <- mol.sum(mol.data = gene.ensprot, id.map = id.map.ensprot) p <- pathview(gene.data = gene.entrez, cpd.data = cpd.kc, pathway.id = '00640', species = 'hsa', same.layer = T, out.suffix = 'gene.entrez.cpd.kc', keys.align = 'y', kegg.native = T, key.pos ='bottomright', sign.pos = 'topright', limit = list(gene = 3, cpd = 3), bins = list(gene = 6, cpd = 6))

不同物種使用時名稱的處理當(dāng)對KEGG處理時,物種是一個棘手但重要的問題。KEGG擁有自己的物種專用命名法和數(shù)據(jù)庫,其中包括大約4800個基因組完整的生物體。換句話說,這些生物體的基因數(shù)據(jù)都可以通過pathview進行映射、可視化和分析。這種全面的物種覆蓋是pathview數(shù)據(jù)集成能力的一個突出特點。然而,KEGG并不以同樣的方式對待所有這些生物體/基因組。KEGG使用NCBI GeneID(或Entrez基因)作為最常見的模型動物的默認(rèn)ID,包括人類、小鼠、大鼠等,以及一些作物,如大豆、葡萄和玉米。另一方面,KEGG對所有其他生物體使用 Locus 標(biāo)記和其他id,包括動物、植物、真菌、原生生物以及細(xì)菌和古生菌。 Pathview帶有一個數(shù)據(jù)矩陣korg,其中包括支持的KEGG物種數(shù)據(jù)和默認(rèn)基因ID的完整列表。讓我們探索korg數(shù)據(jù)矩陣,以便對KEGG物種及其Gene ID的使用有所了解。 data(korg) head(korg)
ktax.id tax.id kegg.code scientific.name common.name entrez.gnodes kegg.geneid [1,] 'T01001' '9606' 'hsa' 'Homo sapiens' 'human' '1' '374659' [2,] 'T01005' '9598' 'ptr' 'Pan troglodytes' 'chimpanzee' '1' '474020' [3,] 'T02283' '9597' 'pps' 'Pan paniscus' 'bonobo' '1' '100989900' [4,] 'T02442' '9595' 'ggo' 'Gorilla gorilla gorilla' 'western lowland gorilla' '1' '101125212' [5,] 'T01416' '9601' 'pon' 'Pongo abelii' 'Sumatran orangutan' '1' '100172878' [6,] 'T03265' '61853' 'nle' 'Nomascus leucogenys' 'northern white-cheeked gibbon' '1' '105739221' ncbi.geneid ncbi.proteinid uniprot [1,] '374659' 'NP_001273380' 'Q8N4P3' [2,] '474020' 'XP_001140087' 'Q1XHV8' [3,] '100989900' 'XP_003811308' NA [4,] '101125212' 'XP_018886437' 'G3QNH0' [5,] '100172878' 'NP_001125944' 'Q5R9G0' [6,] '105739221' 'XP_012359712' 'G1RK33'
#number of species which use Entrez Gene as the default ID sum(korg[,'entrez.gnodes']=='1',na.rm=T) #204
#number of species which use other ID types or none as the default ID sum(korg[,'entrez.gnodes']=='0',na.rm=T) #5041
#new from 2017: most species which do not have Entrez Gene annotation any more na.idx=is.na(korg[,'ncbi.geneid']) sum(na.idx) #4674
從上面的korg的探索中,我們看到4800個KEGG物種中,只有少數(shù)沒有NCBI(Entrez)基因ID或基因ID(注釋)其中的一個。幾乎所有物種都具有默認(rèn)的KEGG基因ID(通常是Locus標(biāo)簽)和Entrez Gene ID注釋。因此,pathview接受所有這些物種的gene.idtype =“kegg”或“Entrez”(不區(qū)分大小寫)。用戶需要確保正確的gene.idtype是特定的,因為除了35種常見物種外,KEGG和Entrez基因ID不同。對于19種,Bioconductor提供基因注釋包。用戶可以自由地輸入gene.data和其他gene.idtype。有關(guān)這些Bioconductor注釋物種和額外基因ID類型的列表允許: data(bods) bods
package species kegg code id.type [1,] 'org.Ag.eg.db' 'Anopheles' 'aga' 'eg' [2,] 'org.At.tair.db' 'Arabidopsis' 'ath' 'tair' [3,] 'org.Bt.eg.db' 'Bovine' 'bta' 'eg' [4,] 'org.Ce.eg.db' 'Worm' 'cel' 'eg' [5,] 'org.Cf.eg.db' 'Canine' 'cfa' 'eg' [6,] 'org.Dm.eg.db' 'Fly' 'dme' 'eg' [7,] 'org.Dr.eg.db' 'Zebrafish' 'dre' 'eg' [8,] 'org.EcK12.eg.db' 'E coli strain K12' 'eco' 'eg' [9,] 'org.EcSakai.eg.db' 'E coli strain Sakai' 'ecs' 'eg' [10,] 'org.Gg.eg.db' 'Chicken' 'gga' 'eg' [11,] 'org.Hs.eg.db' 'Human' 'hsa' 'eg' [12,] 'org.Mm.eg.db' 'Mouse' 'mmu' 'eg' [13,] 'org.Mmu.eg.db' 'Rhesus' 'mcc' 'eg' [14,] 'org.Pf.plasmo.db' 'Malaria' 'pfa' 'orf' [15,] 'org.Pt.eg.db' 'Chimp' 'ptr' 'eg' [16,] 'org.Rn.eg.db' 'Rat' 'rno' 'eg' [17,] 'org.Sc.sgd.db' 'Yeast' 'sce' 'orf' [18,] 'org.Ss.eg.db' 'Pig' 'ssc' 'eg' [19,] 'org.Xl.eg.db' 'Xenopus' 'xla' 'eg'
data(gene.idtype.list) gene.idtype.list
'SYMBOL' 'GENENAME' 'ENSEMBL' 'ENSEMBLPROT' 'UNIGENE' 'UNIPROT' 'ACCNUM' 'ENSEMBLTRANS' 'REFSEQ' 'ENZYME' 'TAIR' 'PROSITE' 'ORF'
所有先前的例子用了人類數(shù)據(jù),其中Entrez Gene是KEGG的默認(rèn)基因ID。Pathview現(xiàn)在(從版本1.1.5開始)顯式處理使用Locus標(biāo)簽或其他基因ID作為KEGG默認(rèn)ID的物種。以下是大腸桿菌菌株K12數(shù)據(jù)的幾個例子。首先,我們使用大腸桿菌K12的默認(rèn)KEGG ID(基因座標(biāo)簽)處理基因數(shù)據(jù)。 eco.dat.kegg <- sim.mol.data(mol.type='gene',id.type='kegg',species='eco',nmol=3000) head(eco.dat.kegg)
查看部分eco.dat.kegg數(shù)據(jù): b3649 b0486 b3312 b3566 b3945 b2936 0.3116743 0.9873035 -1.5734323 1.8890136 0.2399454 0.4556519
p <- pathview(gene.data = eco.dat.kegg, gene.idtype='kegg', pathway.id = '00640', species = 'eco', out.suffix = 'eco.kegg', kegg.native = T, same.layer=T)

我們也可以用與人類相同的方法對大腸桿菌K12的Entrez Gene ID進行基因數(shù)據(jù)處理。 eco.dat.entrez <- sim.mol.data(mol.type='gene',id.type='entrez',species='eco',nmol=3000) head(eco.dat.entrez)
查看部分eco.dat.entrez數(shù)據(jù): 948629 945263 948265 948536 948973 947881 -0.2516281 0.3116743 0.9873035 -1.5734323 1.8890136 0.2399454
p <- pathview(gene.data = eco.dat.entrez, gene.idtype='entrez', pathway.id = '00640', species = 'eco', out.suffix = 'eco.entrez', kegg.native = T, same.layer=T)

基于上述bods數(shù)據(jù),大腸桿菌K12在Bioconductor中有注釋。因此,我們可以進一步研究其基因數(shù)據(jù)與其他ID類型,例如,官方基因符號。在調(diào)用pathview時,首先需要將這些數(shù)據(jù)映射到Entrez Gene ID(如果還沒有),然后再映射到默認(rèn)的KEGG ID(Locus標(biāo)簽)。因此,它比上一個例子花費更長的時間。 egid.eco=eg2id(names(eco.dat.entrez), category='symbol', pkg='org.EcK12.eg.db') eco.dat.symbol <- eco.dat.entrez names(eco.dat.symbol) <- egid.eco[,2] head(eco.dat.kegg) p <- pathview(gene.data = eco.dat.symbol, gene.idtype='symbol', pathway.id = '00640', species = 'eco', out.suffix = 'eco.symbol.2layer', kegg.native = T, same.layer=F)


未注釋物種的處理 (直接用于宏基因組或微生物組數(shù)據(jù))重要的是,當(dāng)數(shù)據(jù)被映射到KEGG ortholog pathways時,pathview可以直接用于宏基因組或微生物組數(shù)據(jù)。來自KEGG中未注釋和未包含的任何新物種(非KEGG物種)的數(shù)據(jù)也可以通過pathview用同樣的方法映射到KEGG ortholog pathways中進行分析和可視化。在下一個例子中,我們首先模擬映射的KEGG ortholog基因數(shù)據(jù)。然后將數(shù)據(jù)作為gene.data輸入,其中species =“ko”。 ko.data=sim.mol.data(mol.type='gene.ko', nmol=5000) head(ko.data)
查看一下ko.data數(shù)據(jù) K01859 K15874 K04035 K16205 K16093 K09447 -0.5198690 0.5832773 -1.3411454 0.5746940 0.2157879 -1.1339828
p <- pathview(gene.data = ko.data, pathway.id = '04112', species = 'ko', out.suffix = 'ko.data', kegg.native = T)

參考https://academic./nar/article/45/W1/W501/3804420 https:///packages/release/bioc/html/pathview.html
|