Heuristic Mining of Hierarchical Genotypes and Accessory Genome Loci in Bacterial Populations

Podsumowanie

This analytical computational platform provides practical guidance for microbiologists, ecologists, and epidemiologists interested in bacterial population genomics. Specifically, the work presented here demonstrated how to perform: i) phylogeny-guided mapping of hierarchical genotypes; ii) frequency-based analysis of genotypes; iii) kinship and clonality analyses; iv) identification of lineage differentiating accessory loci.

Streszczenie

Routine and systematic use of bacterial whole-genome sequencing (WGS) is enhancing the accuracy and resolution of epidemiological investigations carried out by Public Health laboratories and regulatory agencies. Large volumes of publicly available WGS data can be used to study pathogenic populations at a large scale. Recently, a freely available computational platform called ProkEvo was published to enable reproducible, automated, and scalable hierarchical-based population genomic analyses using bacterial WGS data. This implementation of ProkEvo demonstrated the importance of combining standard genotypic mapping of populations with mining of accessory genomic content for ecological inference. In particular, the work highlighted here used ProkEvo-derived outputs for population-scaled hierarchical analyses using the R programming language. The main objective was to provide a practical guide for microbiologists, ecologists, and epidemiologists by showing how to: i) use a phylogeny-guided mapping of hierarchical genotypes; ii) assess frequency distributions of genotypes as a proxy for ecological fitness; iii) determine kinship relationships and genetic diversity using specific genotypic classifications; and iv) map lineage differentiating accessory loci. To enhance reproducibility and portability, R markdown files were used to demonstrate the entire analytical approach. The example dataset contained genomic data from 2,365 isolates of the zoonotic foodborne pathogen Salmonella Newport. Phylogeny-anchored mapping of hierarchical genotypes (Serovar -> BAPS1 -> ST -> cgMLST) revealed the population genetic structure, highlighting sequence types (STs) as the keystone differentiating genotype. Across the three most dominant lineages, ST5 and ST118 shared a common ancestor more recently than with the highly clonal ST45 phylotype. ST-based differences were further highlighted by the distribution of accessory antimicrobial resistance (AMR) loci. Lastly, a phylogeny-anchored visualization was used to combine hierarchical genotypes and AMR content to reveal the kinship structure and lineage-specific genomic signatures. Combined, this analytical approach provides some guidelines for conducting heuristic bacterial population genomic analyses using pan-genomic information.

Wprowadzenie

The increasing use of bacterial whole-genome sequencing (WGS) as a basis for routine surveillance and epidemiological inquiry by Public Health laboratories and regulatory agencies has substantially enhanced pathogen outbreak investigations^1,2,3,4. As a consequence, large volumes of de-identified WGS data are now publicly available and can be used to study aspects of the population biology of pathogenic species at an unprecedented scale, including studies based on: population structures, genotype frequencies, and gene/allele frequencies across multiple reservoirs, geographical regions, and types of environments⁵. The most commonly used WGS-guided epidemiological inquiries are based on analyses using only the shared core-genomic content, where the shared (conserved) content alone is used for genotypic classification (e.g., variant calling), and these variants become the basis for epidemiological analysis and tracing^1,2,6,7. Typically, bacterial core-genome-based genotyping is carried out with multi-locus sequence typing (MLST) approaches using seven to a few thousand loci^8,9,10. These MLST-based strategies encompass mapping of pre-assembled or assembled genomic sequences onto highly curated databases, thereby combining allelic information into reproducible genotypic units for epidemiological and ecological analysis^11,12. For instance, this MLST-based classification can generate genotypic information at two levels of resolution: lower-level sequence types (STs) or ST lineages (7 loci), and higher-level core-genome MLST (cgMLST) variants (~ 300-3,000 loci)¹⁰.

MLST-based genotypic classification is computationally portable and highly reproducible between laboratories, making it widely accepted as an accurate sub-typing approach beneath the bacterial species level^13,14. However, bacterial populations are structured with species-specific varying degrees of clonality (i.e., genotypic homogeneity), complex patterns of hierarchical kinship between genotypes^15,16,17, and a wide range of variation in the distribution of accessory genomic content^18,19. Thus, a more holistic approach goes beyond discrete classifications into MLST genotypes and incorporates the hierarchical relationships of genotypes at different scales of resolution, along with mapping of accessory genomic content onto genotypic classifications, which facilitates population-based inference^18,20,21. Moreover, analyses can also focus on shared patterns of inheritance of accessory genomic loci among even distantly-related genotypes^21,22. Overall, the combined approach enables agnostic interrogation of relationships between population structure and the distribution of specific genomic compositions (e.g., loci) among geospatial or environmental gradients. Such an approach can yield both fundamental and practical information about the ecological characteristics of specific populations that may, in turn, explain their tropism and dispersion patterns across reservoirs, such as food animals or humans.

This systems-based hierarchical population-oriented approach demands large volumes of WGS data for sufficient statistical power to predict distinguishable genomic signatures. Consequently, the approach requires a computational platform capable of processing many thousands of bacterial genomes at once. Recently, ProkEvo was developed and is a freely available, automated, portable, and scalable bioinformatics platform that allows for integrative hierarchical-based bacterial population analyses, including pan-genomic mapping²⁰. ProkEvo allows for the study of moderate-to-large scale bacterial datasets while providing a framework to generate testable and inferable epidemiological and ecological hypotheses and phenotypic predictions that can be customized by the user. This work complements that pipeline in providing a guide on how to utilize ProkEvo-derived output files as input for analyses and interpretation of hierarchical population classifications and accessory genomic mining. The case study presented here utilized the population of Salmonella enterica lineage I zoonotic serovar S. Newport as an example and was specifically aimed at providing practical guidelines for microbiologists, ecologists, and epidemiologists on how to: i) use an automated phylogeny-dependent approach to map hierarchical genotypes; ii) assess the frequency distribution of genotypes as a proxy for evaluating ecological fitness; iii) determine lineage-specific degrees of clonality using independent statistical approaches; and iv) map lineage-differentiating AMR loci as an example of how to mine accessory genomic content in the context of the population structure. More broadly, this analytical approach provides a generalizable framework to perform a population-based genomic analysis at a scale that can be used to infer evolutionary and ecological patterns regardless of the targeted species.

Protokół

1. Prepare input files

NOTE: The protocol is available here - https://github.com/jcgneto/jove_bacterial_population_genomics/tree/main/code. The protocol assumes that the researcher has specifically used ProkEvo (or a comparable pipeline) to get the necessary outputs available in this Figshare repository (https://figshare.com/account/projects/116625/articles/15097503 - login credentials are required - The user must create a free account to have file access!). Of note, ProkEvo automatically downloads genomic sequences from the NCBI-SRA repository and only requires a .txt file containing a list of genome identifications as an input²⁰, and the one used for this work on S. Newport USA isolates is provided here (https://figshare.com/account/projects/116625/articles/15097503?file=29025729).  Detailed information on how to install and use this bacterial genomics platform is available here (https://github.com/npavlovikj/ProkEvo/wiki/2.-Quick-start)²⁰

1. Generate core-genome phylogeny using FastTree²³ as previously described²⁰, which is not part of the bioinformatics platform²⁰. FastTree requires the Roary²⁴ core-genome alignment as an input file. The phylogeny file is named newport_phylogeny.tree  (https://figshare.com/account/projects/116625/articles/15097503?file=29025690).
2. Generate SISTR²⁵ output containing the information regarding serovars classifications for Salmonella and cgMLST variant calling data (sistr_output.csv - https://figshare.com/account/projects/116625/articles/15097503?file=29025699).
3. Generate BAPS file by fastbaps^26,27 containing the BAPS levels 1-6 classification of genomes into sub-groups or haplotypes (fastbaps_partition_baps_prior_l6.csv - https://figshare.com/account/projects/116625/articles/15097503?file=29025684).
4. Generate MLST-based classification of genomes into STs using the MLST program (https://github.com/tseemann/mlst)²⁸ (salmonellast_output.csv - https://figshare.com/account/projects/116625/articles/15097503?file=29025696).
5. Generate ABRicate (https://github.com/tseemann/abricate)²⁹ output as a .csv file containing AMR loci mapped per genome (sabricate_resfinder_output.csv - https://figshare.com/account/projects/116625/articles/15097503?file=29025693).
   NOTE: The user can turn off specific parts of the ProkEvo bioinformatics pipeline (check here for more information - https://github.com/npavlovikj/ProkEvo/wiki/4.2.-Remove-existing-bioinformatics-tool-from-ProkEvo). The analytical approach presented here provides guidelines for how to conduct a population-based analysis after the bioinformatics pipeline has been run. 

2. Download and install the statistical software and integrated development environment (IDE) application

1. Download the most up-to-date freely available version of the R software for Linux, Mac, or PC³⁰. Follow the default installation steps.
2. Download the most up-to-date freely available version of the RStudio desktop IDE here³¹. Follow the default steps for installation.
   ​NOTE: The next steps are included in the available script, including detailed information of code utilization, and should be run sequentially to generate the outputs and figures presented in this work (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/data_analysis_R_code.Rmd). The user may decide to use another programming language to conduct this analytical/statistical analysis such as Python. In that case, use the steps in the scripts as a framework to carry out the analysis. 

3. Install and activate data science libraries

1. Install all data science libraries at once as a first step in the analysis. Avoid installing the libraries every time the script needs to be re-run. Use the function install.packages() for library installation. Alternatively, the user may click on the Packages tab inside of the IDE and automatically install the packages. The code used to install all needed libraries is presented here:
   # Install Tidyverse
   install.packages("tidyverse")
   # Install skimr
   install.packages("skimr")
   # Install vegan
   install.packages("vegan")
   # Install forcats
   install.packages("forcats")
   # Install naniar
   install.packages("naniar")
   # Install ggpubr
   install.packages("ggpubr")
   # Install ggrepel
   install.packages("ggrepel")
   # Install reshape2
   install.packages("reshape2")
   # Install RColorBrewer
   install.packages("RColorBrewer")
   # Install ggtree
   if (!requireNamespace("BiocManager", quietly = TRUE))
     install.packages("BiocManager")
   BiocManager::install("ggtree")
   # Installation of ggtree will prompt a question about installation - answer is "a" to install/update all dependencies
2. Activate all the libraries or packages using the library() function at the beginning of the script, right after installation. Here is a demonstration on how to activate all necessary packages:
   # Activate the libraries and packages
   library(tidyverse)
   library(skimr)
   library(vegan)
   library(forcats)
   library(naniar)
   library(ggtree)
   library(ggpubr)
   library(ggrepel)
   library(reshape2)
   library(RColorBrewer)
3. Suppress outputting the code used for library and package installation and activation by using {r, include = FALSE} in the code chuck, as follows:
   ``` {r, include = FALSE}
   # Install Tidyverse
   install.packages("tidyverse")
   ```
   NOTE: This step is optional but avoids showing chunks of unnecessary code in the final html, doc, or pdf report.
4. For a brief description of the specific functions of all libraries along with some useful links to gather further information, refer to steps 3.4.1-3.4.11.
   1. Tidyverse - use this collection of packages used for data science, including data entry, visualization, parsing and aggregation, and statistical modeling. Typically, ggplot2 (data visualization) and dplyr (data wrangling and modeling) are practical packages present in this library³².
   2. skimr - use this package for generating summary statistics of data frames, including identification of missing values³³.
   3. vegan - use this package for community ecology statistical analyses, such as calculating diversity-based statistics (e.g., alpha and beta-diversity)³⁴.
   4. forcats - use this package to work with categorical variables such as re-ordering classifications. This package is part of the Tidyverse library³².
   5. naniar - use this package to visualize the distribution of missing values across variables in a data frame, by using the viss_miss() function³⁵.
   6. ggtree - use this package for the visualization of phylogenetic trees³⁶.
   7. ggpubr - use this package to improve the quality of ggplot2-based visualizations³⁷.
   8. ggrepel - use this package for text labeling inside of graphs³⁸.
   9. reshape2 - use the melt() function from this package for the transformation of data frames from wide to long format³⁹.
   10. RColorBrewer - use this package to manage colors in ggplot2-based visualizations⁴⁰.
   11. Use the following basic functions for exploratory data analysis: head() to check the first observations in a data frame, tail() to check the last observations of a data frame, is.na() to count the number of rows with missing values across a data frame, dim() to check the number of rows and columns in a dataset, table() to count observations across a variable, and sum() to count the total number of observations or instances.

4. Data entry and analysis

NOTE: A detailed information on each step of this analysis can be found in the available script (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/data_analysis_R_code.Rmd). However, here are some important points to be considered:

1. Do all genomic data entry, including all genotypic classifications (serovar, BAPS, ST, and cgMLST) using the read_csv() function.
2. Rename, create new variables, and select columns of interest from each dataset before multi-dataset aggregation.
3. Don't remove missing values from any independent dataset. Wait until all datasets are aggregated to modify or exclude missing values. If new variables are created for each dataset, then missing values are by default categorized into one of the newly generated classifications.
4. Check for erroneous characters such as hyphens or interrogations marks and replace them with NA (Not applicable). Do the same for missing values.
5. Aggregate data based on the hierarchical order of genotypes (serovar -> BAPS1 -> ST -> cgMLST), and by grouping based on the individual genome identifications.
6. Check for missing values using multiple strategies and deal with such inconsistencies explicitly. Only remove a genome or isolate from the data if the classification is unreliable. Otherwise, consider the analysis being done and remove NAs on a case-by-case basis.
   NOTE: It is highly recommended to establish a strategy to deal with such values a priori. Avoid removing all genomes or isolates with missing values across any variables. For instance, a genome may have ST classification without having cgMLST variant number. In that case, the genome can still be used for the ST-based analysis.
7. Once all datasets are aggregated, assign them to a data frame name or object that can be used in multiple locations in the follow-up analysis, to avoid having to generate the same metadata file for every figure in the paper.

5. Conduct analyses and generate visualizations

NOTE: A detailed description of each step needed to produce all the analysis and visualizations can be found in the markdown file for this paper (https://github.com/jcgneto/jove_bacterial_population_genomics/tree/main/code). Code for each figure is separated in chunks and the entire script should be run sequentially. Additionally, the code for each main and supplementary figure is provided as a separate file (see Supplementary File 1 and Supplementary File 2). Here are some essential points (with snippets of code) to be considered while generating each main and supplementary figures.

1. Use ggtree to plot a phylogenetic tree along with genotypic information (Figure 1).
   1. Optimize the ggtree figure size, including diameter and width of rings, by changing the numerical values inside of the xlim() and gheatmap(width = ) functions, respectively (see example code below).
      tree_plot <- ggtree(tree, layout = "circular") + xlim(-250, NA)
      figure_1 <- gheatmap(tree_plot, d4, offset=.0, width=20, colnames = FALSE)
      NOTE: For a more detailed comparison of programs that can be used for phylogenetic plotting, check this work²⁰. The work highlighted an attempt made to identify strategies to improve ggtree-based visualizations such as decreasing the dataset size, but branch lengths and tree topology were not as clearly discriminating as compared to phandango⁴¹.
   2. Aggregate all metadata into as few categories as possible to facilitate the choice of coloring panel when plotting multiple layers of data with the phylogenetic tree (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/figure_1.Rmd). Conduct the data aggregation based on the question of interest and domain knowledge.
2. Use a bar plot to assess relative frequencies (Figure 2).
   1. Aggregate data for both ST lineages and cgMLST variants to facilitate visualizations. Choose an empirical or statistical threshold used for data aggregation, while considering the question being asked.
   2. For an example code that can be used to inspect the frequency distribution of ST lineages to determine the cut-off see below:
      st_dist <- d2 %>% group_by(ST) %>% # group by the ST column
      count() %>% # count the number of observations
      arrange(desc(n)) # arrange the counts in decreasing order
   3. For an example code showing how minor (low-frequency) STs can be aggregated refer below. As demonstrated below, STs that are not numbered as 5, 31, 45, 46, 118, 132, or 350, are grouped together as "Other STs". Use a similar code for cgMLST variants (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/figure_2.Rmd).
      d2$st <- ifelse(d2$ST == 5, "ST5", # create a new ST column for which minor S Ts are aggregated as Others
       ifelse(d2$ST == 31, "ST31",
        ifelse(d2$ST == 45, "ST45",
         ifelse(d2$ST == 46, "ST46",
          ifelse(d2$ST == 118, "ST118",
      ifelse(d2$ST == 132, "ST132", ifelse(d2$ST == 350, "ST350", "Other STs")))))))
3. Use a nested approach to calculate the proportion of each ST lineage within each BAPS1 sub-group to identify STs that are ancestrally related (belong to the same BAPS1 sub-group) (Figure 3). The code below exemplifies how the ST-based proportion can be calculated across BAPS1 sub-groups (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/figure_3.Rmd):
   baps <- d2b %>% filter(serovar == "Newport") %>% # filter Newport serovars
   select(baps_1, ST) %>% # select baps_1 and ST columns
   mutate(ST = as.numeric(ST)) %>% # change ST column to numeric
   drop_na(baps_1, ST) %>% # drop NAs
   group_by(baps_1, ST) %>% # group by baps_1 and ST
   summarise(n = n()) %>% # count observations
   mutate(prop = n/sum(n)*100) # calculate proportions
4. Plot the distribution of AMR loci across ST lineages using the Resfinder-based gene annotation results (Figure 4).
   NOTE: Resfinder has been widely used in ecological and epidemiological studies⁴². Annotation of protein-coding genes can vary depending on how often databases are curated and updated. If using the suggested bioinformatics pipeline, the researcher can compare AMR-based loci classifications across different databases²⁰. Be sure to check which databases are continually being updated. Do not use out-of-date or poorly curated databases, in order to avoid miscalls.
   1. Use an empirical or statistical threshold to filter out the most important AMR loci to facilitate visualizations. Provide a raw .csv file containing the calculated proportions of all AMR loci across all ST lineages, such as shown here (https://figshare.com/account/projects/116625/articles/15097503?file=29025687).
   2. Calculate the AMR proportion for each ST using the following code (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/figure_4.Rmd):
      # Calculations for ST45
      d2c <- data6 %>% filter(st == "ST45") # filter ST45 data first
      # for ST45, calculate the proportion of AMR loci and only keep proportion greater than 10%
      d3c <- d2c %>% select(id, gene) %>% # select columns
      group_by(id, gene) %>% # group by id and gene
      summarize(count = n()) %>% # count observations
      mutate(count = replace(count, count == 2, 1)) %>% # replace counts equal to 2 with 1 to only consider one copy of each gene (duplications may not be reliable), but the researcher can decide to exclude or keep them. If the researcher wants to exclude them, then use the filter(count != 2) function or else leave as is
      filter(count <= 1) # filter counts below or equal to 1
      d4c <- d3c %>% group_by(gene) %>% # group by gene
      summarize(value = n()) %>% # count observations
      mutate(total = table(data1$st)[6]) %>% # get the total counts of st mutate(prop = (value/total)*100) # calculate proportions
      ​d5c <- d4c %>% mutate(st = "ST45") # create a st column and add ST information
   3. After calculations are done for all STs, combine datasets as one data frame, using the following code:
      # Combine datasets
      d6 <- rbind(d5a, d5b, d5c, d5d, d5e, d5f, d5g, d5h) # row bind datasets
   4. To export the .csv file containing the calculated proportions, use the code:
      # Export data table containing ST and AMR loci information
      abx_newport_st <- d6 write.csv(abx_newport_st,"abx_newport_st.csv", row.names = FALSE)
   5. Before plotting the AMR-based distribution across ST lineages, filter the data based on a threshold to facilitate visualizations, as shown below:
      # Filter AMR loci with proportion higher than or equal to 10%
      d7 <- d6 %>% filter(prop >= 10) # determine the threshold empirically or statistically
5. Plot the core-genome phylogeny along with the hierarchical genotypic classifications and AMR data in a single plot using ggtree (Figure 5).
   1. Optimize the figure size inside ggtree using the abovementioned parameters (see step 5.1.1.).
   2. Optimize visualizations by aggregating variables, or using binary classification such as gene presence or absence. The more features are added to the plot, the harder the coloring selection process becomes (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/figure_5.Rmd).
      NOTE: Supplementary figures - detailed description of the entire code can be found here (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/data_analysis_R_code.Rmd).
6. Use a scatter plot in ggplot2, without data aggregation, to display the distribution of ST lineages or cgMLST variants while highlighting the most frequent genotypes (Supplementary Figure 1) (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/supplementary_figure_s1.Rmd). 
7. Do a nested analysis to assess the composition of ST lineages through the proportion of cgMLST variants in order to get a glimpse of the ST-based genetic diversity, while identifying the most frequent variants and their genetic relationships (i.e., cgMLST variants that belong to the same ST shared an ancestor more recently than those belonging to distinct STs) (Supplementary Figure 2) (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/supplementary_figure_s2.Rmd). 
8. Use community ecology metric, namely Simpson's D index of diversity, to measure the degree of clonality or genotypic diversity of each of the major ST lineages⁴³ (Supplementary Figure 3).
   1. Calculate the index of diversity across ST lineages at different levels of genotypic resolution including BAPS level 1 through 6 and cgMLST. Below is the code example on how to do this calculation at the BAPS level 1 (BAPS1) of genotypic resolution:
      # BAPS level 1 (BAPS1)
      # drop the STs and BAPS1 with NAs, group by ST and BAPS1 and then calculate Simpson's index
      baps1 <- data6 %>%
      select(st, BAPS1) %>% # select columns
      drop_na(st, BAPS1) %>% # drop NAs
      group_by(st, BAPS1) %>% # group by columns
      summarise(n = n()) %>% # count observations
      mutate(simpson = diversity(n, "simpson")) %>% # calculate diversity
      group_by(st) %>% # group by column
      summarise(simpson = mean(simpson)) %>% # calculate the mean of the index
      melt(id.vars=c("st"), measure.vars="simpson",
      variable.name="index", value.name="value") %>% # covert into long format
      mutate(strat = "BAPS1") # create a strat column
      NOTE: A more genetically diverse population (i.e., more variants at different layers of genotypic resolution) has a higher index at the cgMLST level and produces an increasing index-based values going from BAPS level 2 to 6 (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/supplementary_figure_s3.Rmd). 
9. Examine the degree of genotypic diversity of ST lineages by plotting the relative frequency of BAPS sub-groups at all levels of resolution (BAPS1-6) (Supplementary Figure 4). The more diverse the population is, the sparser the distribution of BAPS sub-groups (haplotypes) becomes going from BAPS1 (lower level of resolution) to BAPS6 (higher level of resolution) (https://github.com/jcgneto/jove_bacterial_population_genomics/blob/main/code/supplementary_figure_s4.Rmd). 

Wyniki

By utilizing the computational platform ProkEvo for population genomics analyses, the first step in bacterial WGS data mining is comprised of examining the hierarchical population structure in the context of a core-genome phylogeny (Figure 1). In the case of S. enterica lineage I, as exemplified by the S. Newport dataset, the population is hierarchically structured as follows: serovar (lowest level of resolution), BAPS1 sub-groups or haplotypes, ST lineages, and cg...

Dyskusje

The utilization of a systems-based heuristic and hierarchical population structure analysis provides a framework to identify novel genomic signatures in bacterial datasets that have the potential to explain unique ecological and epidemiological patterns²⁰. Additionally, the mapping of accessory genome data onto the population structure can be used to infer ancestrally-acquired and/or recently-derived traits that facilitate the spread of ST lineages or cgMLST variants across reservoirs

Materiały

Name	Company	Catalog Number	Comments
amr_data_filtered			https://figshare.com/account/projects/116625/articles/14829225?file=28758762
amr_data_raw			https://figshare.com/account/projects/116625/articles/14829225?file=28547994
baps_output			https://figshare.com/account/projects/116625/articles/14829225?file=28548003
Core-genome phylogeny			https://figshare.com/account/projects/116625/articles/14829225?file=28548006
genome_sra			https://figshare.com/account/projects/116625/articles/14829225?file=28639209
Linux, Mac, or PC			any high-performance platform
mlst_output			https://figshare.com/account/projects/116625/articles/14829225?file=28547997
sistr_output			https://figshare.com/account/projects/116625/articles/14829225?file=28548000
figshare credentials are required for login and have access to the files