In mirnavazquez/RbiMs: R Tools for Reconstruncting Bin Metabolisms

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

First, load the rbims package.

library(rbims)

The second thing to do would be to read the KO_picrust2_data. The meta file is a tab-separated file containing the name of your ids (bins or sites) and any extra information you would like to use for visualization.

Read metadata

The metadata table can be read in various formats (csv, tsv, txt, xlsx); you will need to use the corresponding function to read the type of file you have. meta and try.

The example data is from an article already published link and were sediment core samples that were collected from three zones of mangrove ecosystem: Basin, Fringe and Impaired. The samples were taken during two seasonal conditions, dry and flood seasons, at three sediment depths: surface (5), mid-layer(20), and deeper sediment layers (40), to assess the impact of degradation level or zone, seasonal variation, and sediment depth on microbial communities.

library(readr)
meta<-read_tsv("../inst/extdata/Metadatos-MANGLAR")

head(meta)

Let's chaghe the names of ids for visualization purposes:

library(dplyr)
library(stringr)
meta <- meta %>%
  mutate(
    Bin_name = str_extract(Bin_name, "(?<=zr2502_)(\\d+)(?=_R1)") 
  )

head(meta)

Let's read again the KO_picrust2 data setting profile=F.

KO_picrust2_profile_long <- read_picrust2("../inst/extdata/KO_metagenome_out/KO_pred_metagenome_unstrat.tsv",
                                          database = "KO",
                                          profile = F)

head(KO_picrust2_profile_long)

The same change in Bin names:

KO_picrust2_profile_long <- KO_picrust2_profile_long %>%
  mutate(
    Bin_name = str_extract(Bin_name, "(?<=zr2502_)(\\d+)(?=_R1)")
  )
head(KO_picrust2_profile_long)

Map to the KEGG database

Then map the KO to the rest of the features of the KEGG and rbims database.

ko_bin_mapp<-mapping_ko(KO_picrust2_profile_long)

Metabolism subsetting

In this example, we will use energy metabolism to explore the rest of the functions.

Create a vector with the metabolism of interest.

Other_energy<-c("Fermentation", "Carbon fixation", "Methane metabolism", 
                "Sulfur metabolism", "Nitrogen metabolism")

Use get_subset_pathway to subset the table using the cycles and the information of the energy metabolism.

library(tidyr)
Energy_metabolisms<-ko_bin_mapp %>%
  drop_na(Cycle, Pathway_cycle) %>%
  get_subset_pathway(Cycle, Other_energy)

Plot the information using the plot_heatmap function. The argument order_y will order the rows according to a metabolic feature; in this case, we order the pathways_cycle according to cycle.

plot_heatmap(tibble_ko=Energy_metabolisms, 
             y_axis=Pathway_cycle,
             order_y = Cycle,
             split_y = TRUE,
            data_experiment = meta,
             order_x = "zone",
             analysis = "KEGG",
             calc="Percentage")

Example of exploring a specific KEGG pathway

specific<-c("map05416", "map00363", "map00604", "map00908", "map00941")

Now, let’s extract the profile associated with that metabolic pathway.

library(tidyr)
specific_subset<-ko_bin_mapp%>%
  drop_na(Pathway) %>%
  get_subset_pathway(Pathway, specific)

head(specific_subset)

Now, let’s create a plot of the presence and absence of the different KO associated with that pathway. Besides presence and absence, it is possible to plot abundance or the percentage of genes within certain pathways (See plot_bubble, calc argument).

plot_bubble(tibble_ko = specific_subset,
            x_axis = Bin_name, 
            y_axis = KO,
            analysis="KEGG",
            data_experiment = meta,
            color_character = zone,
            calc="Binary",
            range_size = c(1,10))