Enriquecimento de vias KEGG com clusterProfiler - lmigueel/Bioinformatica GitHub Wiki
1. Sobre
O serviço KEGG FTP não está disponível gratuitamente para acadêmicos desde 2012, e há muitos pacotes de software usando dados de anotação KEGG desatualizados. O pacote clusterProfiler suporta o download da última versão online dos dados KEGG usando o site KEGG, que está disponível gratuitamente para usuários acadêmicos. O caminho e o módulo KEGG são suportados no clusterProfiler.
2. Instalação
Instalação e uso do pacote é feito no R.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
3. Uso geral
Nesta parte vou comentar coisas bem legais que podemos fazer com o clusterProfiler. Existem milhares de opções, mas vou destacar algumas, principalmente o enriquecimento KEGG.
3.1 Organismos suportados pelo clusterProfiler
O pacote clusterProfiler oferece suporte a todos os organismos que possuem dados de anotação KEGG disponíveis no banco de dados KEGG. O usuário deve passar a abreviatura do nome para o parâmetro organismo. A lista completa de organismos com suporte do KEGG pode ser acessada em http://www.genome.jp/kegg/catalog/org_list.html. O banco de dados de ortologia KEGG (KO) também é suportado especificando organism = "ko".
Quando você faz análises de metagonômica (até em metatranscriptômica) e identifica os KEGG ONTHOLOGY (KO) para seus genes, você pode usar qualquer função abaixo com a opção mencionada: organism = "ko".
O pacote clusterProfiler fornece a função search_kegg_organism() para ajudar na busca de organismos suportados.
library(clusterProfiler)
search_kegg_organism('ece', by = 'kegg_code')
# kegg_code scientific_name common_name
#366 ece Escherichia coli O157:H7 EDL933 (EHEC) <NA>
Caso você busque pelo nome científico do organismo, ele pode voltar várias espécies que fazem parte da mesma família. Por exemplo em E. coli.
ecoli <- search_kegg_organism('Escherichia coli', by='scientific_name')
dim(ecoli)
## [1] 65 3
head(ecoli)
## kegg_code scientific_name PERGUNTA
## 606 eco Escherichia coli K-12 MG1655 Escherichia coli K-12 MG1655
## 607 ecj Escherichia coli K-12 W3110 Escherichia coli K-12 W3110
## 608 ecd Escherichia coli K-12 DH10B Escherichia coli K-12 DH10B
## 609 ebw Escherichia coli K-12 BW2952 Escherichia coli K-12 BW2952
## 610 ecok Escherichia coli K-12 MDS42 Escherichia coli K-12 MDS42
## 611 ece Escherichia coli O157: H7 EDL933 EHEC
3.2 Análise de enriquecimento KEGG
Para identificar as vias enriquecidas na sua lista de genes, o pacote faz uma análise de enriquecimento hipergeométrico ACESSE AQUI, e identifica quais as vias mais representadas. O comando é bem simples, basta ter o ID específico para o organismo que você desejar. O cluster profile suporta como tipo de IDS: kegg, ncbi-geneid, ncbi-proteinid ou uniprot.
O legal aqui é ter a lista de genes que deseja verificar e até alguma informação extra, como expressão diferencial. Nesse caso podemos pintar os genes das vias segundo uma escala de cor. Vou criar um dataset de genes de levedura (S. cerevisiae) com o valor em TPM de expressão. O ID utilizado é UNIPROT. Como default ele espera o ID KEGG de cada gene, que também terá uma coluna no dataset. Vamos usar a função enrichKEGG().
genes <- c('P32465', 'P11412','P46969','P22091','P23254','P14540',
'C7GTA0','P16861','P00924','P23250','P38858','P08417',
'P39533','P53598','P00950','P00360','Q00711','P28834',
'Q02455','P47116','P01119','P14693')
kegg_id <- c('YHR094C','YNL241C','YJL121C','YLR354C','YPR074C','YKL060C',
'YGR240C','YGR254W','YIL119C','YHR163W','YPL262W','YJL200C',
'YOR142W','YKL152C','YJL052W','YKL148C','YNL037C','YKR095W',
'YJR059W','YOR101W','YHR083W','YDH123C')
exp <- c(110,12,3,0,1.2,4,5,12,13,130,11,20,40,30,44,12,1,2,2,2,1,10)
df <- data.frame(genes = genes, kegg_id = kegg_id, expressao=exp)
kk <- enrichKEGG(gene = df$genes,
organism = 'sce',
keyType = 'uniprot',
pvalueCutoff = 0.05)
head(kk)
# ID Description GeneRatio BgRatio
#sce01200 sce01200 Carbon metabolism 13/17 112/2077
#sce01110 sce01110 Biosynthesis of secondary metabolites 14/17 348/2077
#sce00030 sce00030 Pentose phosphate pathway 6/17 28/2077
#sce01230 sce01230 Biosynthesis of amino acids 9/17 125/2077
#sce00680 sce00680 Methane metabolism 4/17 26/2077
#sce00010 sce00010 Glycolysis / Gluconeogenesis 5/17 55/2077
# pvalue p.adjust qvalue
#sce01200 3.254813e-14 7.160589e-13 4.453955e-13
#sce01110 4.583045e-09 5.041349e-08 3.135767e-08
#sce00030 3.808262e-08 2.792726e-07 1.737102e-07
#sce01230 1.260850e-07 6.934678e-07 4.313436e-07
#sce00680 4.118969e-05 1.812347e-04 1.127297e-04
#sce00010 5.263884e-05 1.930091e-04 1.200535e-04
# geneID
#sce01200 P11412/P46969/P23254/P14540/P16861/P00924/P38858/P08417/P53598/P00950/P00360/Q00711/P28834
#sce01110 P11412/P46969/P23254/P14540/P16861/P00924/P38858/P08417/P39533/P53598/P00950/P00360/Q00711/P28834
#sce00030 P11412/P46969/P23254/P14540/P16861/P38858
#sce01230 P46969/P23254/P14540/P16861/P00924/P39533/P00950/P00360/P28834
#sce00680 P14540/P16861/P00924/P00950
#sce00010 P14540/P16861/P00924/P00950/P00360
# Count
#sce01200 13
#sce01110 14
#sce00030 6
#sce01230 9
#sce00680 4
#sce00010 5
Observamos que identificamos as vias KEGG enriquecidas. É fornecido o número de genes identificados na sua lista e na via. Utilizem
de corte o valor de p.adjust <= 0.05. Algumas das vias enriquecidas pra levedura são:
- sce01200 Carbon metabolismo
- sce01110 Biosynthesis of secondary metabolites
- sce00030 Pentose phosphate pathway
- sce00010 Glycolysis / Gluconeogenesis
3.3 Análise de enriquecimento de sub-módulos do KEGG
O Module KEGG é uma coleção de unidades de função definidas manualmente. Em algumas situações, os Módulos KEGG têm uma interpretação mais direta. São subunidades das vias principais.
mkk <- enrichMKEGG(gene = df$genes,
organism = 'sce',
keyType = 'uniprot',
pvalueCutoff = 1,
qvalueCutoff = 1)
head(mkk)
# ID
#M00001 M00001
#M00004 M00004
#M00003 M00003
#M00009 M00009
#M00002 M00002
#M00011 M00011
# Description
#M00001 Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate
#M00004 Pentose phosphate pathway (Pentose phosphate cycle)
#M00003 Gluconeogenesis, oxaloacetate => fructose-6P
#M00009 Citrate cycle (TCA cycle, Krebs cycle)
#M00002 Glycolysis, core module involving three-carbon compounds
#M00011 Citrate cycle, second carbon oxidation, 2-oxoglutarate => oxaloacetate
# GeneRatio BgRatio pvalue p.adjust qvalue
#M00001 5/14 24/542 0.0001697465 0.0007272584 NA
#M00004 4/14 14/542 0.0002424195 0.0007272584 NA
#M00003 4/14 16/542 0.0004276291 0.0008552582 NA
#M00009 4/14 25/542 0.0025922727 0.0038884090 NA
#M00002 3/14 15/542 0.0052131654 0.0062557985 NA
#M00011 3/14 16/542 0.0063175200 0.0063175200 NA
# geneID Count
#M00001 P14540/P16861/P00924/P00950/P00360 5
#M00004 P11412/P46969/P23254/P38858 4
#M00003 P14540/P00924/P00950/P00360 4
#M00009 P08417/P53598/Q00711/P28834 4
#M00002 P00924/P00950/P00360 3
#M00011 P08417/P53598/Q00711 3
Nesse caso, vários sub-módulos específicos da Glicólise estão enriquecidos. Essas sub-unidades são bem legais pra pensar em vias menores.
3.4 Visualize vias KEGG enriquecidas
O pacote enrichplot implementou vários métodos para visualizar termos enriquecidos. A maioria deles são métodos gerais que podem ser usados para GO, KEGG, MSigDb e outras anotações de conjuntos de genes. Aqui, apresentamos as funções clusterProfiler::browseKEGG() e pathview::pathview() para ajudar os usuários a explorar a via KEGG enriquecida com genes interessados.
Para visualizar o caminho KEGG, o usuário pode usar a função browseKEGG, que irá abrir o navegador da web e destacar genes enriquecidos.
Vamos ver a Pentose Fosfato:
browseKEGG(kk, 'sce00030')
Ele irá abrir a via KEGG no browser. Os retângulos (genes) em verde são os de levedura. Acesse o link aqui para vocês verem: PENTOSE FOSFATO - YEAST. Devo salientar uma coisa importante. Atualmente, o browserKEGG funciona apenas com o Kegg ID padrão (ENTREZID para humanos). Quando você usa o Kegg ID, os genes da sua lista ficam em vermelho.
3.5 Pathview
Podemos também podem usar a função pathview() do pacote pathview (Luo e Brouwer 2013) para visualizar caminhos KEGG enriquecidos identificados pelo pacote clusterProfiler (Yu et al. 2012).
O exemplo a seguir ilustra como visualizar a via 'sce00030', que foi enriquecida em nossa análise anterior. Vamos colorir agora segundo a expressão de cada gene (por isso que disse que talvez agregue as análises).
Caso você precise converter para Kegg ID os seus genes em UNIPROT, use o mapping do próprio Uniprot. ACESSE AQUI.
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("pathview")
library("pathview")
sce00030 <- pathview(gene.data = as.character(df$genes),
pathway.id = "sce00030",
species = "sce",
gene.idtype="UNIPROT",
limit = list(gene=max(abs(df$expressao)), cpd=1)
)
Ele irá salvar um arquivo sce00030.pathview.png no seu diretório atual no R. Caso queira um heatmap da sua via mais elaborado, aconselho você a extrair os genes somente deste pathway e achar o maximo da escala por ele.
# resultado da variavel kk na terceira posicao é a via sce00030.
# Eu pego e monto um dataset a partir da lista do resultado do enrichment
# Altero o nome da coluna pra ficar bonito
# Faço um left_join (PROCV do excel) que vai juntar os datasets em um unico
#install.packages("tidyverse")
library(tidyverse)
genes_ppp <- str_split(kk[3][8](/lmigueel/Bioinformatica/wiki/8), "/")
genes_ppp
df_ppp <- data.frame(ppp = genes_ppp)
colnames(df_ppp) <- c("ppp")
df_ppp
df_final_ppp <- left_join(df_ppp,df,
by=c("ppp" = "genes"))
sce00030 <- pathview(gene.data = as.character(df_final_ppp$ppp),
pathway.id = "sce00030",
species = "sce",
gene.idtype="UNIPROT",
limit = list(gene=max(abs(df_final_ppp$expressao)), cpd=1)
)
Espero que tenham gostado pessoal!
Deve ajudar bem vocês!
Bibliografia
Para mais informações, acesse o manual e capítulo do clusterProfile AQUI.