Installation

Recommended: bioconda

The simplest and recommend way to install protal is from bioconda. We recommend installing protal in a separate conda environment.

conda create -n protal protal -c bioconda

Run the following to test your installation

conda activate protal
protal --help

Compile from source

Protal can also be installed by compiling the sources from github. Instructions are provided in the github Readme.

Binary releases

Precompiled static binaries are available in the github releases for protal. These are useful if you are unable to use conda on your cluster. If your cluster contains separate nodes and the nodes with internet are different from the ones running protal, this could lead to library issues, which could be resolved by using the static binaries. Protal uses pigz for compression, so this should also be installed if using these static binaries.

In the directory you downloaded protal binary to, run the following to test your installation.

./protal --help

Download the database

Protal requires a copy of the database. Download and extract this using:

wget https://protal.earlham.ac.uk/data/dbs/protal-db-r226-0.5.1a.tar.gz && \
tar -xzvf protal-db-r226-0.5.1a.tar.gz

This extracts to a folder called protal_0.5.1_r226, which is the path you pass to --db in the examples below.

A reduced database, protal-mini-db-r226-0.5.1a.tar.gz, is also available on the downloads page. It extracts to protal_0.5.1_r226_mini and uses roughly five times less memory at near-identical precision and somewhat lower sensitivity. Use it in place of the full database by pointing --db at that folder instead.

Instead of passing --db on every run, you can point the environment variable $PROTAL_DB_PATH at the extracted database folder and omit the argument.

Usage

Now that you have installed protal and downloaded the database, we will illustrate usage, we will use data from .... (context)

First we will show a simple example of running protal on a single sample and looking at the taxonomic profile. Afterwards, we will run protal on four samples, and look at strain resolved output as well as the taxonomic profiles.

Download test data

Download paired-end short read sequencing from four samples using

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/002/SRR8797712/SRR8797712_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/002/SRR8797712/SRR8797712_2.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/003/SRR8797713/SRR8797713_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/003/SRR8797713/SRR8797713_2.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/000/SRR8797950/SRR8797950_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/000/SRR8797950/SRR8797950_2.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/005/SRR8797945/SRR8797945_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR879/005/SRR8797945/SRR8797945_2.fastq.gz

Run a single sample

We will run protal for one sample with the following command

protal \
    --db protal_0.5.1_r226/ \
    -1 SRR8797712_1.fastq.gz \
    -2 SRR8797712_2.fastq.gz \
    --prefix SRR8797712 \
    --outdir ./ \
    --threads 4

The arguments are:

--db - The path to the protal database
-1 - Path to the forward reads
-2 - Path to reverse reads
--prefix - Prefix for this sample's output files (optional). If omitted, protal derives a prefix from the longest common prefix of the two read filenames, which only works if both files are in the same folder.
--outdir - Output directory
--threads - Number of threads

Protal currently supports paired-end reads only, so -1 and -2 are both mandatory. Taxonomic profiling and strain analysis both run by default; use --no_profile and --no_strains to switch them off.

Or condensed in one line

protal --db protal_0.5.1_r226/ -1 SRR8797712_1.fastq.gz -2 SRR8797712_2.fastq.gz --prefix SRR8797712 --outdir ./ --threads 4

Run with multiple samples using a mapping file

For more complex runs, it can be tedious to provide all parameters on the command line. Additionally, running lots of sample individually is inefficient as the database will need to be loaded for each run.

Therefore, protal allows you to provide a mapping file, a tabular format which provides paths to input and output files. The mapping file equivalent for the previous example looks like this:

#OUTPUT_DIR	/path/to/output/
#SAM_OUTPUT_DIR	/path/to/output/alignments
#PROFILE_OUTPUT_DIR	/path/to/output/profiles
#STRAIN_OUTPUT_DIR	/path/to/output/strains
#MISC_OUTPUT_DIR	/path/to/output/misc
#INPUT_DIR	/path/to/input/
#SAMPLEID	FIRST	SECOND	SAM	PREFIX	PROFILE
sample1	SRR8797712_1.fastq.gz	SRR8797712_2.fastq.gz	sample1.sam	sample1	sample1.profile

Note that the mapping file uses tab characters between columns and fields, not spaces. Help is available below, and via the --map_help argument.

Header lines are indicated with a #, and define locations for input and output files. For output locations, these headers define a directory to write all results of a type to. If a header is omitted, the corresponding directory defaults to a subfolder of #OUTPUT_DIR named alignments, profiles, strains or misc. Unknown # keys are silently ignored, so watch for typos -- a misspelled header does not raise an error, it just leaves the default in place.

The final header line gives a list of column headers, and is followed by details for each individual sample. The columns are

SAMPLEID - A row label for readability. Note that protal does not actually read this column; the sample name used throughout the output (MSA headers, meta.tsv, statistics files) is taken from PREFIX.
FIRST - Mandatory. First read file, in fastq format. Protal will try to find this file in {INPUT_DIR}/{FIRST}.
SECOND - Mandatory. Second read file for paired end reads. As above.
PREFIX - Mandatory. The sample name, and the prefix used before filenames relating to this sample.
SAM - Optional. The name for the SAM output, stored in {SAM_OUTPUT_DIR}/{SAM}. If the column is absent, the SAM defaults to {OUTPUT_DIR}/{PREFIX}.sam -- note that this is #OUTPUT_DIR, not #SAM_OUTPUT_DIR.
PROFILE - Optional. The name for the profile output, stored in {PROFILE_OUTPUT_DIR}/{PROFILE}. If the column is absent, the profile defaults to {OUTPUT_DIR}/{PREFIX}.profile, again relative to #OUTPUT_DIR.

Give every sample its own SAM and PROFILE name, otherwise the samples overwrite each other's output.

If you need to, you could leave the *_OUTPUT_DIR headers empty and specify the full path for each output file in the sample rows. This can get complicated quickly, and we would recommend not doing so.

All we need to do now is provide the mapping file to protal.

protal \
    --db protal_0.5.1_r226/ \
    --map test.map \
    --threads 4

Run four example samples using a mapping file

We will now use a mapping file to run all four samples we downloaded earlier. The mapping for for all these samples looks like:

#OUTPUT_DIR	results
#SAM_OUTPUT_DIR	results/alignments
#PROFILE_OUTPUT_DIR	results/profiles
#STRAIN_OUTPUT_DIR	results/strains
#MISC_OUTPUT_DIR	results/misc
#INPUT_DIR	data
#SAMPLEID	FIRST	SECOND	SAM	PREFIX	PROFILE
SRR8797712	SRR8797712_1.fastq.gz	SRR8797712_2.fastq.gz	SRR8797712.sam	SRR8797712	SRR8797712.profile
SRR8797713	SRR8797713_1.fastq.gz	SRR8797713_2.fastq.gz	SRR8797713.sam	SRR8797713	SRR8797713.profile
SRR8797945	SRR8797945_1.fastq.gz	SRR8797945_2.fastq.gz	SRR8797945.sam	SRR8797945	SRR8797945.profile
SRR8797950	SRR8797950_1.fastq.gz	SRR8797950_2.fastq.gz	SRR8797950.sam	SRR8797950	SRR8797950.profile

Either save this to a file called example.map, or you can also download a copy here.

We can profile these samples with protal using:

protal \
    --db protal_0.5.1_r226/ \
    --map example.map \
    --threads 4

The output files are described in more detail in the Output details section, but we will take a look at how to use two of the important outputs: species abundance profiles, and strain information.

Useful options

Protal has sensible defaults, so most runs need nothing beyond --db, the read files and --threads. The options below are the ones users reach for most often; protal --help lists them all, and protal --full_help additionally shows the developer options.

Option	Description
`--no_profile`	Only produce alignments, skip taxonomic profiling.
`--no_strains`	Stay on species level. Do not output SNPs or MSAs.
`--knob 0.5`	Detection threshold. Taxa with a random forest probability >= `knob` are reported. Higher values favour precision, lower values favour sensitivity.
`--model`	Use an alternative PMML model file instead of `model.xml` in the database directory.
`--profile_dir`	Override the profile output directory. Takes precedence over the map file.
`--profile_only`	Comma separated list of existing sam files to profile without re-running the alignment.
`--map_range 1-10`	Run protal on a subset of the rows in a `--map` file. The first entry is 1 and the end is inclusive.
`--msa_species`	Restrict MSA output to a single species (`s__Genus_species`) or a comma separated list.
`--force`	Redo the alignment even if the sam file already exists.
`--strain_preset`	Aggressiveness of the qcmsa post-filter: `strict`, `default` or `sensitive`.
`--no_qcmsa`	Disable the qcmsa post-filter (it is on by default).
`--qcmsa_args`	Forward arguments verbatim to qcmsa, e.g. `--qcmsa_args "--gene-min-hcov 0.5"`. Run `qcmsa --help` for the full list.
`--verbose`	Verbose program output.

Species profiles

The estimated abundance of species in each sample can be found in results/profiles/{SAMPLEID}.profile. Look at the first few lines of the profile for sample SRR8797713

> head -3 results/profiles/SRR8797713.profile
RS_GCF_000585335.1      d__Bacteria|p__Synergistota|c__Synergistia|o__Synergistales|f__Synergistaceae|g__Cloacibacillus|s__Cloacibacillus evryensis 0.0157048
GB_GCA_022846095.1      d__Bacteria|p__Bacillota_A|c__Clostridia|o__Oscillospirales|f__Oscillospiraceae|g__UMGS1202|s__UMGS1202 sp900549565 0.00441174
RS_GCF_003201875.1      d__Bacteria|p__Bacillota_A|c__Clostridia|o__Lachnospirales|f__Lachnospiraceae|g__Hungatella|s__Hungatella effluvii  0.00372886

Each line gives information about a single species: the acccession, GTDB taxonomy, and relative abundance of the taxon.

Combining profiles into an abundance table

Each sample has a separate file giving their abundance profile. For further analysis it is usually convenient to merge these into a single table. Protal ships with protal_profile_utils, which is installed alongside protal, so it is already on your PATH in the conda environment:

protal_profile_utils merge --input results/profiles/*.profile > abundance.tsv

This writes a single tab-separated table to standard output, with one row per taxon and one column per sample. Add --resolve-samples if your profiles come from different directories and the filenames alone are not enough to tell the samples apart.

Basic analysis and visualisation

This notebook on Google colab provides code for combining profiles, and visualising protal output.

Strain resolution

In addition to species profiles, protal creates multiple sequence alignments which can be used to create within-species phylogenies. These are written by default to the strains subdirectory in the results directory.

For each species, protal first writes its own unfiltered alignment as {species}.raw.msa.fna. The qcmsa post-filter then runs on that alignment -- this happens by default -- and writes the final {species}.msa.fna alongside a matching {species}.partition.txt. Use {species}.msa.fna for downstream analysis. It removes multi-allelicity outlier samples and genes, applies gene coverage and depth thresholds, and cleans up poorly covered sites, all of which improves the resulting phylogenies. Keep {species}.raw.msa.fna if you want to apply your own filtering instead.

The post-filter can be tuned without re-running the alignment:

--strain_preset strict|default|sensitive shifts how aggressively outliers are removed.
--qcmsa_args "..." forwards any qcmsa option verbatim, for example --qcmsa_args "--gene-min-hcov 0.5 --gene-min-samples 5". Run qcmsa --help for the full list.
--no_qcmsa disables the post-filter entirely, in which case only {species}.raw.msa.fna is produced.

Building a tree from multiple sequence alignments

This notebook on Google colab provides code examples for building and plotting a phylogenetic tree from a multiple sequence alignment.

Output details

Protal output is divided into four main types, typically written to different directories in the results. The contents of each of these is described here, with primary outputs marked in bold.

`profiles`

Relative abundance profiles of species in each sample. For each sample x, the following files will be created

x.profile - Relative abundance of species in the sample. Tab separated without column names. Columns are species_accession, gtdb_taxonomy, relative_abundance.
x.profile.log - Every candidate taxon, including those that did not pass the --knob threshold. Each row carries the prediction flag (1/0), the random forest probability, the accession, the lineage, the abundance, and the coverage and variance features the model was scored on.
x.profile.gene.log - Per-gene coverage detail behind each taxon in x.profile.log.
x.profile.genes.log - Per-gene table with a header row: Truth, Predicted, Probability, TaxID, Lineage, TaxVCOV, TaxaxAbundance, GeneID, GeneRefLength, TotalReads, TotalMappedLength, MAPQ, UniqueMers, UniqueTwoMers, UniqueTwoMersReads, UniqueTwoMerReads, ANI, VCov, VCovExp, HCovExp, HCovObs, HCovObsRel, Consistency.
x.profile.truth_annotated - Only written when the developer option --profile_truth is given. Annotates the profile taxa with TP/FP against a supplied truth set.

`alignments`

For each sample x

x.sam.gz - Alignment of reads to marker genes. Protal writes an uncompressed sam and then compresses it with pigz, using --threads.
x.sam.gz.err - Diagnostic output written by the profiler while reading the sam file.

`strains`

For each species y

y.msa.fna - Multiple sequence alignment of genes between strains, after the qcmsa post-filter. This is the file to use for phylogenies.
y.partition.txt - Relationship between marker genes and coordinates in y.msa.fna, in RAxML style. Coordinates are recomputed to match the filtered alignment.
y.raw.msa.fna - Protal's own unfiltered MSA, before qcmsa.
y.raw.partition.txt - Partitions for y.raw.msa.fna.
y.meta.tsv - Per-sample and per-gene alignment metrics, with a header row.
y.snp_stats.tsv - Per-sample SNP retention and filter breakdown: total variant positions, passing and filtered SNPs, retained SNPs, insertions and deletions, and how many variants each filter removed.
y.qcmsa_summary.tsv - Machine readable decision log from qcmsa: the parameters used, and how many samples, genes and sites went in and came out.

y.meta.tsv columns

sample
gene_id
vertical_coverage (sum of all reads / gene length)
counts_vcov1 - number of positions with at least vertical coverage 1
counts_vcov2 - number of positions with at least vertical coverage 2
multi_allelic - number of multiallelic positions (see misc)
filtered - number of filtered positions (see misc)
multi_rate_vcov1 - multi_allelic / counts_vcov1
filtered_rate_vcov1 - filtered / counts_vcov1
multi_rate_vcov2 - multi_allelic / counts_vcov2
filtered_rate_vcov2 - filtered / counts_vcov2
median_vcov - median vertical coverage across the gene
hcov - horizontal coverage (breadth) of the gene
gene_length - reference length of the gene
mean_vcov_nonzero - mean vertical coverage over covered positions only
median_vcov_nonzero - median vertical coverage over covered positions only

`misc`

The output in misc uses -1 and -2 to indicate specific situations

-2 is reserved for when a gene is not available for a species
-1 is reserved for when gene exists for a species, but has not been hit by a read

For each species y

y.hcov.tsv - Horizontal coverage (breadth) of marker genes in each sample. First column is the sample, second the total coverage across all genes, the remaining 120 columns are coverage per gene.
y.statistics.tsv - Summary of reads aligned to this reference in each sample, has column headings, tab separated format

Several files are related to the number of SNPs and how these are filtered. A variant is called as a SNP when it is supported by at least --snp_min_cov reads (default 2) and passes the quality check, which is an OR: the cumulative phred score across supporting reads must be at least --snp_min_phred_sum (default 90), or the mean base quality must be at least --snp_min_mean_qual (default 15). Two further filters apply: --snp_min_af sets a minimum allele frequency (default 0, i.e. off), and a strand-bias filter requires at least one supporting read from each strand -- this is on by default and can be switched off with --snp_no_strand.

y.snps_total.tsv - Count of number of SNPs passing the filtering criteria for species y. The first column is the sample, second the number of SNPs, the remaining 120 are count of SNPs per gene.
y.snps_filtered.tsv - Count of SNPs which did not pass the filtering criteria, structured as y.snps_total.tsv
y.snps_multiallelic.tsv - The number of biallelic, triallelic and tetraallelic positions after applying filtering criteria, structure as y.snps_total.tsv