The COMBAT-TB Workbench: Making powerful TB bioinformatics accessible

Data management & User interface

The IRIDA platform, a project of the Public Health Agency of Canada (PHAC-NML), is a web application written in Java, that provides a user-friendly web-interface for sequencing data and metadata storage, workflow execution and results visualisation. We adopted IRIDA as the basis of the COMBAT-TB workbench because of its integration with the Galaxy platform, its proven track record in public health bioinformatics (PHAC-NML has analysed over 200,000 pathogen samples using IRIDA (44)) and the fact that it is open source. IRIDA stores sequence samples on disk and sample metadata in a MariaDB database (9). Sequence data is shared between IRIDA and the Galaxy analysis platform (reducing data duplication).

Scientific workflows and IRIDA plugins

The Workbench uses Galaxy for its bioinformatics workflow composition and execution. Two workflows, for M. tuberculosis sample analysis and phylogeny, were implemented in Galaxy. Building workflows in Galaxy involves connecting Galaxy tools to construct an analysis workflow where the Galaxy tools (also known as tool wrappers) themselves connect command line bioinformatics tools to the Galaxy framework. Building our workflows involved updating some Galaxy tools (Trimmomatic (10), snippy (11) and snp-sites (12)) and writing new Galaxy tools (snp-dists (13), TB-Profiler (5), tb_variant_filter (14) and TB Variant Report (15)). In addition to the work on Galaxy tools, the command line tb_variant_filter and TB Variant Report tools were created as part of the COMBAT-TB project.

The irida-wf-ga2xml tool (16) was updated to be able to work with recent versions of Galaxy and it was used to build IRIDA plugin skeletons. These plugin skeletons were further developed into IRIDA plugins for both workflows, in the case of the M. tuberculosis sample analysis involving the addition of an interface between the Galaxy workflow output and the metadata stored for each sequence sample. The M. tuberculosis Sample Report and M. tuberculosis Phylogeny plugins are hosted in Github repositories (https://github.com/COMBAT-TB/irida-plugin-tb-sample-report and https://github.com/COMBAT-TB/irida-plugin-tb-phylogeny respectively) and deployed into IRIDA as part of COMBAT-TB Workbench deployment process.

Deployment

The COMBAT-TB Workbench application (available at: https://github.com/COMBAT-TB/irida-galaxy-deploy) deploys IRIDA, the COMBAT-TB IRIDA plugins, the MariaDB database IRIDA uses for metadata storage and Galaxy using Docker containers (17). The Docker containers are orchestrated using docker-compose (18). This allows the entire Workbench to be installed by users without advanced Linux systems administration knowledge and hides the complexity of the underlying software from the user.

Integration with external data storage

Data can be loaded into the COMBAT-TB Workbench via the web interface or via the command line IRIDA uploader tool. This allows data from external storage (for example the storage of a sequencing facility) to be loaded into the Workbench in bulk.

Analysis of data from MDR M. tuberculosis in Indonesia

Tania et al (22) collected M. tuberculosis samples from 30 patients with confirmed pulmonary tuberculosis treated in four hospitals in the Western region of the Indonesian island of Java. Phenotypic drug susceptibility testing (DST) was performed on the cultured samples and DNA was isolated, sequenced and deposited in the NCBI Sequence Read Archive (SRA) with BioProject access number PRJNA633244. Thirty samples were retrieved from EBI European Nucleotide Archive (ENA) and uploaded to the COMBAT-TB Workbench using the irida-uploader (23) command line tool. Per sample quality control was performed automatically (using FastQC(24)) after sample upload. Uploading and quality control took 9 and 3 minutes respectively.

The inferred ancestral M. tuberculosis reference genome produced by Comas et al (25) was downloaded from Zenodo (26) and uploaded to the COMBAT-TB Workbench as the default reference genome for the sequencing project. While the H37Rv reference genome can be used, we used the M. tuberculosis inferred ancestral reference as Goig et al (27) have previously shown that this genome is equidistant, in terms of sequence variants, from all known M. tuberculosis lineages and thus provides a superior reference to H37Rv (NC_000962.3) for variant calling, especially if that variant calling is going to be used for phylogeny construction.

Per sample reporting

The M. tuberculosis Sample Report pipeline (see Materials and Methods) was run on all 30 samples. Detailed output from this pipeline is made available in a per-sample analysis report (Supplementary Figure 1) that includes a full report on variants identified in the sample (annotated using information from the COMBAT-TB NeoDB (28)), drug resistance prediction (from TB-Profiler) and quality control information on read mapping.

The pipeline updated the metadata associated with each sample, with the result that after pipeline completion each sample was annotated with percentage of reads mapped to genome, drug resistance prediction, drug resistance associated variants and lineage assigned (Figure 3).

• Download figure
• Open in new tab

Figure 3:

Sample Line List

Examining the read mapping outputs shows that only 9.63% and 1.47% of reads from samples SRR12416824 and SRR12416842 respectively mapped against the M. tuberculosis genome. Further investigation with kraken2 (29) using the Standard database from 14 April 2020 showed that the majority (63.55%) of reads from SRR12416824 were classified as belonging to the Mycobacterium avium complex (MAC) and the majority of reads (75.39%) from SRR12416842 were classified as the non-tubercular mycobacterium Mycolicibacterium fortuitum. These samples were thus excluded and not used in subsequent analyses.

As shown in Table 1, the results were broadly concordant with those found by Tania et al. Some small differences are likely a result of different versions of the TB-Profiler software used. Without our own analysis, updating the TB-Profiler version from 2.8.4 to version 3.0.6 reduced the discordance between streptomycin resistance predicted by the COMBAT-TB Workbench and that predicted by MGIT by 2 samples. This illustrates the evolution of drug resistance prediction from WGS data. As noted previously, the COMBAT-TB Workbench pipeline uses snippy for mapping and variant calling and applies TB-Profiler only to the results of the mapping step. Our results are, however, wholly concordant with running the complete TB-Profiler pipeline.

The lower degree of concordance found between WGS and MGIT resistance testing for ethambutol can be explained by the presence of variants impacting amino acid 306 of the embB gene in 7 out of the 11 samples that were reported as resistant by WGS but susceptible by MGIT. These variants (M306V in 2 and M306I in 5) lead to a change in the MIC of ethambutol (30).

Phylogeny on all samples

The 28 samples that previously passed quality control were selected in the web interface and submitted to the M. tuberculosis phylogeny pipeline (https://github.com/COMBAT-TB/irida-plugin-tb-phylogeny). This pipeline computes a Maximum Likelihood Phylogeny using the single nucleotide variants identified in each sample (Figure 4).

• Download figure
• Open in new tab

Figure 4:

COMBAT TB Workbench phylogeny viewer

The IRIDA Advanced Visualization view (Figure 5) allows metadata from the sequence analysis project’s sequencing store to be associated with tips (i.e. samples) in the phylogeny view. In addition to the metadata computed by the M. tuberculosis Sample Report pipeline, an Excel spreadsheet was generated associating each sample ID with the sample IDs and hospital collection sites identified in the Tania et al. paper. This spreadsheet was used to load additional metadata into the IRIDA project.

• Download figure
• Open in new tab

Figure 5:

Advanced Phylogeny Visualisation with Metadata columns

The distance between hospital sites varied from 10km (Mampang Prapatan Hospital, Jakarta to Cempaka Putih Islamic Hospital, Jakarta) to 107km (Drajat Prawiranegara Hospital, Serang to M. Goenawan Partowidigdo Pulmonary Hospital, Bogor). When visualised this way it was apparently that there was no clear relationship between phylogenetic relationship (and thus lineage) and hospital site, illustrating that the outbreaks occurring in the region were circulating in areas broader than hospital catchment areas.

Analysis of runtimes

The runtimes of the steps in the analyses are listed in Table 2. Galaxy scheduled the execution of analysis steps in parallel when data dependencies allowed. The similarities of the runtimes between the sample processing and the phylogeny pipelines are to some extent due to the fact that IRIDA currently starts all analyses from sequence reads rather than from assembled genomes. This limitation is currently being addressed by the IRIDA development community. (31)

Analysis of data from M. tuberculosis in Spain

Xu et al (32) collected M. tuberculosis samples from 117 patients and analysed them to understand dynamics of TB transmission in the Valencia Region, in Spain. While the key element of their analysis is identifying individual transmission patterns, their work also provides a convenient larger dataset to examine the performance of the COMBAT-TB Workbench

Uploading of the 117 samples took 21 minutes. Per sample analysis of the 117 samples took a total of 6 hours. A phylogeny was computed using the samples. Phylogeny generation took just over 8 hours and visualisation with the advanced phylogeny viewer revealed clustering by M. tuberculosis lineage (Figure 7) as expected. Unfortunately Xu et al did not include the output of their phylogeny construction in their results, so a direct comparison of the computed phylogenetic trees is not possible.

• Download figure
• Open in new tab

Figure 7:

Advanced Phylogeny viewer with Xu et al data (metadata shows lineage)

Implementation of the M. tuberculosis Sample Report pipeline

The TB sample report workflow starts by running Trimmomatic (v. 0.38.1) with the Sliding Window trimmer, truncating reads when the average quality within a four base window drops below quality score 30, followed by the Minimum length trimmer that discards reads shorter than 20 bases. Only reads which remain in read pairs after quality trimming are used in subsequent analysis.

After quality trimming, sequence reads are mapped to a user supplied reference genome. The workflow requires a genome with the same coordinate scheme as the M. tuberculosis H37Rv reference (RefSeq NC000962). Mapping and variant calling is done using snippy (v. 4.4.5), a Perl-based pipeline that combines the bwa-mem (v. 0.7.17) (33) mapper and the freebayes (v. 1.3.2) (34) variant caller. Snippy has been shown (35) to produce good quality variant calls in M. tuberculosis when run with default parameters. Samtools (v. 1.9) (36) flagstat is run to provide an overview of statistics from the mapping process.

After variant calling, variants are annotated with SnpEff (v. 4.3) (37) using the H37Rv reference genome annotation. Variants are then filtered using tb_variant_filter (v. 0.1.3). While snippy performs quality based filtering of the variants it predicts, this tool offers a variety of filtering options commonly used in M. tuberculosis variant filtering. In our workflow it filters out variants in the PE/PPE gene regions (38), in the repetitive and insertion section regions identified by UVP (39) and those with lower than 30 supporting reads or within 5 base pairs of an indel. These filtering options like all tool options in the workflow can optionally be altered by the user.

In parallel to the SnpEff annotation and variant filtering steps, the mapped reads are provided to TBProfiler (v. 2.8.4) which performs its own variant calling and lineage and drug resistance prediction.

Finally the filtered, annotated variants and the TBProfiler results are fed to tb_vcf_report (v. 0.1.7) which produces a report further annotated with information from the COMBAT-TB eXplorer database (28) in both text and HTML formats.

The final reports provided to the user are the variant reports from tb_vcf_report, text and JSON format reports from TBProfiler, variants in VCF format from SnpEff and mapping statistics from samtools flagstats. These reports include both user readable and raw data suitable for further downstream analysis. The metadata stored in the sample line list is updated with mapping %, M. tuberculosis lineage and spoligotype information and drug resistance information.

Implementation of the M. tuberculosis Phylogeny pipeline

The workflow in the phylogeny module starts with quality filtering of samples using fastp (v. 0.19.5) (40) with default settings. The filtered reads are then aligned to the user provided reference using snippy and predicted variants are filtered with tb_variant_filter as described above. In addition, only single nucleotide variants (SNVs) are retained as the phylogeny software used in the workflow cannot extract meaningful information from indels.

For each sample, the identified sequence variants are inserted into the reference genome, yielding one sequence per sample, of the same length as the reference but with variants from the sample inserted. These are concatenated into a multiple sequence alignment (MSA), which is used as input to snp_dists (v. 0.6.3). Variant and constant sites are identified using snp_sites (v. 2.5.1) and the multiple sequence alignment is filtered to retain only variant sites. A phylogeny is then built using IQ-TREE (v. 1.5.5.1) (41–43).

The final report includes the SNP distance matrix and (Newick format) tree. The tree is also visualised in a phylogeny viewer and can also be displayed with associated sample metadata.

The SNP distance matrix, the phylogeny and the filtered variants (in VCF format) are presented to the user as output.