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-r214-0.1.0a.tar.gz && \
tar -xzvf protal-db-r214-0.1.0a.tar.gz

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-db-r214-0.1.0a/ \
    -1 SRR8797712_1.fastq.gz \
    -2 SRR8797712_2.fastq.gz \
    -3 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 * -3 - Sample name for output * --outdir - Output directory * --threads - Number of threads

Or condensed in one line

protal --db protal-db-r214-0.1.0a/ -1 SRR8797712_1.fastq.gz -2 SRR8797712_2.fastq.gz -3 SRR8797712 --outdir ./ --threads 4 --profile

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:

#STRAIN_OUTPUT_DIR  /path/to/output/strain
#OUTPUT_DIR /path/to/output/
#SAM_OUTPUT_DIR /path/to/output/sam
#PROFILE_OUTPUT_DIR /path/to/output/profiles
#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.

The final header line gives a list of column headers, and is followed by details for each individual sample. The columns are * SAMPLEID - A unique name for each sample. * FIRST - First read file, in fastq format. Protal will try to find this file in {INPUT_DIR}/{FIRST}. * SECOND - Second read file for paired end reads. As above. * SAM - The name for the SAM output. This will be stored in {SAM_OUTPUT_DIR}/{SAM}. * PREFIX - A prefix used before filenames relating to this sample. * PROFILE - The name for the profile output. This will be stored in {PROFILE_OUTPUT_DIR}/{PROFILE}.

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_r214/ \
    --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:

#STRAIN_OUTPUT_DIRr results/strain
#OUTPUT_DIR results
#SAM_OUTPUT_DIR results/sam
#PROFILE_OUTPUT_DIR profiles
#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_r214/ \
    --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.

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. The python script from the protal repo can do this for you:

wget https://...{url}...
python merge_profiles.py *.profile

This will create tab-separated tables abundance.tsv and taxononmy.tsv.

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 strain subdirectory in the results directory. For each species, several files are created; the primary ones of interest for analysis are the multiple sequence alignment files. These are {species}.msa.fna and {species}.processed.msa.fna. More detail is available in the output details section The processed MSA has some experimental additional filtering applied to try to improve results. You can use either for further analysis.

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.gene.log - Unknown
x.profile.genes.log - Unknown
x.profile.log - Unknown
x.profile.truth_annotated - Unknown

`alignments`

For each sample x

x.sam.gz - Alignment of reads to marker genes
x.sam.gz.err - Unknown

`strains`

For each species y

y.msa.fna - Multiple sequence alignment of genes between strains
y.msa.processed.fna - As above, but with some additional experimental processing applied to improve alignment
y.partition.txt - Relationship between marker genes and coordinates in y.msa.fna
y.meta.tsv - Variety of information about the alignments

y.meta.tsv columns

Sample name
GeneId
ObservedVerticalCoverage (sum of all reads / gene length)
Number of positions with at least vertical coverage 1
Number of positions with at least vertical coverage 2
Number of multiallelic positions (see misc)
Number Filtered positions (see misc)
MultiAllelic positions / Column 4 (rate of multiallelic positions)
Filtered positions / Column 4 (rate of filtered snps)
MultiAllelic positions / Column 5
Filtered positions / Column 5

`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. The filtering criteria are that a single allele needs to have at least 2 observations AND a phred quality sum across observations > 60 to be counted as valid allele.

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