10 Inferencia metabólica

Dra. Mirna Vazquez Rosas Landa

03 de agosto de 2022

Vamos a utilizar prodigal para predecir las proteínas:

for i in /home/mirna/07.Bins/Genoma/01.Bins_named/*.fa ; do prodigal -i $i -o $i.txt -a $i.faa ; done

Veamos un poco la salida

grep ">" *.faa

10.1 KEEG

Okay, ahora vamos a utilizar kofam_scan para anotar las proteínas.

Vamos a dividirnos en equipos para hacer la anotación de los bins!

Esto es como lo haríamos para uno:

/home/programs/DB/kofam/kofam_scan-1.3.0/exec_annotation -o /home/mirna/08.Kofamscan/htn_bins_63.fa.faa.txt /home/mirna/07.Bins/Proteoma/htn_bins_63.fa.faa  --report-unannotated  --cpu 4

En un loop lo podríamos hacer así:

for i in *.faa ; do  /home/programs/DB/kofam/kofam_scan-1.3.0/exec_annotation -o /home/mirna/08.Kofamscan/$i.txt $i  --report-unannotated  --cpu 4; done

Okay antes de correrlo, vamos a crear dos carpetas nuevas y unas ligas simbólicas.

mkdir -p 08.Kofamscan/{01.Proteomas,02.KO_results}
cd 08.Kofamscan/
ln -s /home/mirna/07.Bins/Proteoma/*.faa 01.Proteomas 

Ajusta el loop para correr KofamScan con los genomas que te tocan :)

10.2 Explorando el metabolismo con rbims.

Vamos a hacer una exploración rápida del metabolismo con rbims.

Okay iniciamos con la librería de Rbims

library(rbims)
library(tidyverse)

Ahora, vamos a leer los resultados de KEEG y mapearlos con el resto de la base de datos de KEEG

htn_mapp<-read_ko("08.Kofamscan/02.KO_results/") %>%
    mapping_ko()

Okay, vamos a enfocarnos en los metabolismos encargados de la obtención de energía.

Overview<-c("Central Metabolism", "Carbon Fixation", 
            "Nitrogen Metabolism", "Sulfur Metabolism", "Fermentation", 
            "Methane Metabolism")
Energy_metabolisms_htn<-htn_mapp %>%
  drop_na(Cycle) %>%
  get_subset_pathway(rbims_pathway, Overview) 

Vamos a visualizar los datos.

plot_bubble(tibble_ko = Energy_metabolisms_htn,
            x_axis = Bin_name, 
            y_axis = Pathway_cycle,
            analysis="KEGG",
            calc="Percentage",
            range_size = c(1,10),
            y_labs=FALSE,
            x_labs=FALSE)  

Okay, incorporemos metadatos, por ejemplo la taxonomía.

Metadatos<-read_delim("11.GTDBTK/Metadatos.txt", delim="\t")

Hagamos un plot

plot_bubble(tibble_ko = Energy_metabolisms_htn,
            x_axis = Bin_name, 
            y_axis = Pathway_cycle,
            analysis="KEGG",
            data_experiment = Metadatos,
            calc="Percentage",
            color_character = Class,
            range_size = c(1,10),
            y_labs=FALSE,
            x_labs=FALSE) 

Ahora, vamos a explorar una sola vía

Secretion_system_htn<-htn_mapp %>%
  drop_na(Cycle) %>%
  get_subset_pathway(Cycle, "Secretion system")

Y hagamos un heatmap

plot_heatmap(tibble_ko=Secretion_system_htn, 
             y_axis=Genes,
             analysis = "KEGG",
             calc="Binary")

Ahora agreguemos metadatos

plot_heatmap(tibble_ko=Secretion_system_htn, 
             y_axis=Genes,
             data_experiment = Metadatos,
             order_x = Phylum,
             analysis = "KEGG",
             calc="Binary")

plot_heatmap(tibble_ko=Secretion_system_htn, 
             y_axis=Genes,
             data_experiment = Metadatos,
             order_y = Pathway_cycle,
             order_x = Phylum,
             analysis = "KEGG",
             calc="Binary")

10.3 Anotación con InterproScan

for i in $(ls *.faa); do sed -i "s/\*//g" $i; done
for i in $(ls *.faa); do interproscan.sh -cpu 4 -goterms -pa -i $i > Log_Interpro_Scan_$i.txt; done 

Okay… vamos a juntar todo en un solo.

cat *.tsv > htn_interpro.tsv

Vamos a R

library(rbims)
library(tidyverse)

interpro_Interpro_profile<-read_interpro(
  data_interpro = "09.Interpro/01.Proteomas/htn_interpro.tsv", 
  database="INTERPRO", profile = T) %>%
  filter(!str_detect(INTERPRO, "-"))

important_INTERPRO<-get_subset_pca(tibble_rbims=interpro_Interpro_profile, 
                                cos2_val=0.95,
                                analysis="INTERPRO")

plot_heatmap(important_INTERPRO, y_axis=INTERPRO, analysis = "INTERPRO", distance = T)

plot_heatmap(important_INTERPRO, y_axis=INTERPRO, analysis = "INTERPRO", distance = F)

10.3.1 Finalmente podemos escribir esto en una tabla.

write_metabolism("09.Interpro/01.Proteomas/htn_interpro.tsv", 
                 "08.Kofamscan/02.KO_results/")