FishTaco: Functional Shifts’ Taxonomic Contributors

Other Input Data Options

FishTaco relates any taxonomic abundances to any community-level gene abundances in the context of a case-control study. You may have prior information on the functions encoded by each taxon, or not. The taxonomic and functional abundances can be generated by any method, meaning that many different data types and processing methods may be appropriate for generating input data for FishTaco, each with different pros and cons. A few possibilities include:

• 16S rRNA sequencing with PICRUSt functional predictions
• Paired 16S rRNA sequencing and metagenomics datasets
• Metagenomic data that has been annotated with both taxonomic and functional information (e.g. using Metaphlan2 and Humann2)

Other analysis options

You can examine the other possible options to provide for a FishTaco analysis by using the “-h” flag:

run_fishtaco.py -h

Some of these are briefly described below.

• Inferring and/or providing genomic content (-gc, -inf): A file specifying the genomic KO content for each taxon can be provided with the “-gc” flag. If the “-inf” flag is also provided, FishTaco will still infer a genomic profile for each taxon, but will use the provided content as a prior for this inference.

• Analyzing different levels of functions (-map_function_level): By default, FishTaco analyzes taxonomic contributors to the differential abundance of KEGG pathways. To decompose a more specific functional category, FishTaco can also analyze KEGG modules, or custom groups of functions (see the documentation).

• Single-taxa versus multi-taxa permutations (-assessment): By default, FishTaco analyzes the contribution of each taxon to each function by permuting the observed abundances of varying subsets of taxa across samples. Supplying “-assessment single_taxa” will cause FishTaco to instead permute the abundances of each single taxon independently, which will speed up the runtime but produce less accurate results. See the publication for more details.

• Analyzing only a specific subset of functions (-single_function_filter, -multi_function_filter_list): Specify a single function (KEGG pathway or module) or list of functions to be analyzed with FishTaco.

Using the FishTaco web server

To use the web server, navigate to http://elbo-spice.gs.washington.edu/shiny/FishTacoPlot/, and select the following options on the right-hand menu:

• Upload FishTaco results file -> Upload results file -> Choose file -> fishtaco_out_main_output_SCORE_wilcoxon_ASSESSMENT_single_taxa.tab
• Select taxonomy file -> Upload custom file -> Choose file -> taxonomy_vaginal_fishtaco.txt
• Experiment with the various options for functions to display and plot settings.

Using the FishTacoPlot R package:

0. Install RStudio (here) if you haven’t done so yet, and open a new session.

1. First, install the FishTacoPlot package from GitHub, using the devtools package:

if(!requireNamespace("devtools", quietly = T)){ # Install devtools if you haven't already
  install.package("devtools")
}
devtools::install_github("borenstein-lab/fishtaco-plot")

2. Next, run the following code, which will plot the taxonomic contributors for the top 5 most differentially abundant functions. The code below also uses the KEGGREST package to label each bar plot with the names of each pathway (rather than the ID numbers). If you do not have the KEGGREST package installed, you can install it from the Bioconductor repository (see also the pre-workshop instructions). This code also assumes all data and results are in the same directory - change file paths as needed.

library(FishTacoPlot)
library(KEGGREST)

##Number of functions to include in plot
n_functions = 5

## Obtain list of the n most differentially abundant functions
top_functions = fread("fishtaco_out_STAT_DA_function_SCORE_wilcoxon_ASSESSMENT_single_taxa.tab")[order(abs(StatValue), decreasing = T)][1:n_functions, Function]

## Obtain the names of the top pathways from KEGG
top_functions_names = sapply(top_functions, function(x){
  return(keggGet(x)[[1]]$NAME)
})

# Make the plot
p = MultiFunctionTaxaContributionPlots(input_dir=getwd(), input_prefix="fishtaco_out",
input_taxa_taxonomy="taxonomy_vaginal_fishtaco.txt", sort_by="predicted_da", plot_type="bars", add_predicted_da_markers=TRUE, add_original_da_markers=TRUE, add_case_control_line = T, 
add_names_in_bars = T, input_function_filter_list = top_functions, add_facet_labels = T)

## Adjust the plot formatting
p = p  + scale_x_continuous(breaks=seq(1:n_functions), labels = top_functions_names) +
guides(fill=guide_legend(nrow=7)) + ylab("Wilcoxon test statistic (W)") +
theme(plot.title=element_blank(), axis.title.x=element_text(size=12,colour="black",face="plain"),
axis.text.x=element_text(size=10,colour="black",face="plain"), axis.title.y=element_blank(),
axis.text.y=element_text(size=9,colour="black",face="plain"), axis.ticks.y=element_blank(),
axis.ticks.x=element_blank(), panel.grid.major.x = element_line(colour="light gray"), panel.grid.major.y = element_line(colour="light gray"),
panel.grid.minor.x = element_line(colour="light gray"), panel.grid.minor.y = element_line(colour="light gray"), panel.background = element_rect(fill="transparent",colour=NA), panel.border = element_rect(fill="transparent",colour="black"), legend.background=element_rect(colour="black"), legend.title=element_text(size=10), legend.text=element_text(size=8,face="plain"),
legend.key.size=unit(0.8,"line"), legend.spacing=unit(0.1,"line"), legend.position="bottom")

## Display the plot
p

As you examine the FishTaco plots, recall that the upper bar represents case-associated taxa, while the lower bar represents control-associated taxa. Some questions to consider:

• Which taxa are associated with bacterial vaginosis in the example data?
• Which functions are determined more by shifts in the case-associated taxa, and which are determined more by the control-associated?
• What are the possible implications of these differences in driver taxa?