03_sesion2.Rmd

# Inferencia metabólica

Dra. Mirna Vazquez Rosas Landa

03 de agosto de 2022


Vamos a utilizar [prodigal](https://github.com/hyattpd/Prodigal) para predecir las proteínas:

```{bash, eval=FALSE}
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

```{bash, eval=FALSE}
grep ">" *.faa
```

## KEEG

Okay, ahora vamos a utilizar [kofam_scan](https://github.com/takaram/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:

```{bash, eval=FALSE}
/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í: 

```{bash, eval=FALSE}
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.

```{bash, eval=F}
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 :)

## Explorando el metabolismo con rbims.

Vamos a hacer una exploración rápida del metabolismo con [rbims](https://mirnavazquez.github.io/RbiMs/index.html).

Okay iniciamos con la librería de Rbims

```{r, eval=F}
library(rbims)
library(tidyverse)
```

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

```{r, eval=F}
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. 

```{r, eval=F}
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.

```{r, eval=F}
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. 

```{r, eval=F}
Metadatos<-read_delim("11.GTDBTK/Metadatos.txt", delim="\t")
```

Hagamos un plot

```{r, eval=F}
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

```{r, eval=F}
Secretion_system_htn<-htn_mapp %>%
  drop_na(Cycle) %>%
  get_subset_pathway(Cycle, "Secretion system")
```

Y hagamos un heatmap

```{r, eval=F}
plot_heatmap(tibble_ko=Secretion_system_htn, 
             y_axis=Genes,
             analysis = "KEGG",
             calc="Binary")
```

Ahora agreguemos metadatos

```{r, eval=F}
plot_heatmap(tibble_ko=Secretion_system_htn, 
             y_axis=Genes,
             data_experiment = Metadatos,
             order_x = Phylum,
             analysis = "KEGG",
             calc="Binary")
```

```{r, eval=F}
plot_heatmap(tibble_ko=Secretion_system_htn, 
             y_axis=Genes,
             data_experiment = Metadatos,
             order_y = Pathway_cycle,
             order_x = Phylum,
             analysis = "KEGG",
             calc="Binary")
```

## Anotación con [InterproScan](https://interproscan-docs.readthedocs.io/en/latest/)

```{bash, eval=FALSE}
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.

```{bash, eval=FALSE}
cat *.tsv > htn_interpro.tsv
```

Vamos a R

```{r, eval=FALSE}
library(rbims)
library(tidyverse)
```

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

```{r, eval=FALSE}
important_INTERPRO<-get_subset_pca(tibble_rbims=interpro_Interpro_profile, 
                                cos2_val=0.95,
                                analysis="INTERPRO")
```

```{r, eval=FALSE}
plot_heatmap(important_INTERPRO, y_axis=INTERPRO, analysis = "INTERPRO", distance = T)

```

```{r, eval=FALSE}
plot_heatmap(important_INTERPRO, y_axis=INTERPRO, analysis = "INTERPRO", distance = F)
```

### Finalmente podemos escribir esto en una tabla.

```{r, eval=FALSE}
write_metabolism("09.Interpro/01.Proteomas/htn_interpro.tsv", 
                 "08.Kofamscan/02.KO_results/")
```