README.md

phap - Phage Host Analysis Pipeline

A snakemake workflow that wraps various phage-host prediction tools.

• Uses Singularity containers for execution of all tools. When possible (i.e. the image is not larger than a few Gs), tools and their dependencies are bundled in the same container. This means you do not need to get models or any other external databases.
• Calculates Last Common Ancestor of all tools per contig.

Current tools

Tool (source)	Publication/Preprint	Comments
HTP	Gałan W. et al., 2019	ok
RaFAh	Coutinho F. H. et al., 2020	ok
vHuLK	Amgarten D. et al., 2020	ok
VirHostMatcher-Net	Wang W. et al., 2020	ok
WIsH	Galiez G. et al., 2017	ok (unnecessary?)

Installation

Software dependencies

To run the workflow your will need

• snakemake > 5.x (developed with 5.30.1)
• singularity >= 3.6 (developed with 3.6.3)

The following python packages are also required to be installed and available in the execution environment

• biopython >= 1.78 (developed with 1.78)
• ete3 >= 3.1.2 (developed with 3.1.2)

Conda environment

It is recommended to use a conda environment. The file environment.txt can be used to recreate the complete environment used during development.

The provided environment.txt contains an explicit list of all packages, produced with

conda list -n phap --explicit > environment.txt

This ensures all packages are exactly the same versions/builds, so we minimize the risk of running into dependencies issues

To get a working environment

# Clone this repo and get in there
$ git clone https://git.science.uu.nl/papanikos/phap.git
$ cd phap

# Note the long notation --file flag; -f will not work.
$ conda create -n phap --file=environment.txt

# Activate it - use the name you gave above, if it is different
$ conda activate phap

# The (phap) prefix shows we have activated it
# Check the snakemake version
(phap) $ snakemake --version
5.30.1

Data dependencies

• RaFaH, vHULK, HTP

For these tools there is no need to pre-download and setup anything - all data and software dependencies required for running them are bundled within their respective singularity image.

• VirHostMatcher-Net, WIsH

Databases and models need to be downloaded from the VirHostMatcher data repo (see here). WIsH models for the 62,493 host genomes used in their paper are also provided and are used here for WIsH predictions.

NCBI Taxonomy

The ete3.NCBITaxa class is used to get taxonomy information and calculate the LCA of all predictions, when possible. This requires a taxa.sqlite to be available either in its default location ( $HOME/.ete3toolkit/taxa.sqlite ) or provided in the config. See more about that on http://etetoolkit.org/docs/latest/tutorial/tutorial_ncbitaxonomy.html

Singularity containers

Definition files, along with documentation of how to use them to build the containers are in resources/singularity.

The pre-built containers are all available through the standard singularity library.

These are pulled at runtime (or used from cache).

Alternatively, you can pull all .sif files from the cloud, store them locally and use these in an offline mode (see below).

Configuration

The config_template.yaml file provided with this repo has all available configurable options. Short explanations are provided as commented blocks for each option.

A configuration for the workflow must be available as a config.yaml within the config directory. A separate my_config.yaml overriding the options in the default config.yaml can be supplied at runtime, e.g.

$ snakemake --configfile=path/to/my_config.yaml \
<rest of snakemake options>

You can either copy the config_template.yaml and rename it to config.yaml or make your edits straight on the template and rename it to config.yaml.

All fields included in the template must be specified, unless otherwise stated in the comment.

Input data

The tools wrapped in this workflow expect phage sequences as input. You should try to make sure that the input sequences you want to analyze correspond to phage genomes/contigs (or at least viruses).

You can probably input any valid fasta file but the GIGO concept is probably applicable.

A separate workflow to identify phage/viral genomes/contigs is What the Phage.

The current workflow can handle multiple samples. For each sample, all viral contigs to be analyzed should be provided as a single multifasta (can be gzipped). A mapping between sample ids and their corresponding fasta file is provided as a samplesheet (see below).

Sample sheet

You must define a samplesheet with two tab (\t) separated columns. The header line must contain two fields, sample fasta. Values from the sample column must be unique and are used as sample identifiers. Their corresponding fasta values must be valid paths to (multi)fasta files with the phage sequences for that sample.

An example

$ cat samples.tsv
sample	fasta
s01	/path/to/s01.fna
s02	/path/to/another.fna.gz

Note

There is no need to follow any convention for the fasta file name to reflect the sample id. The values in the sample column are the ones to worry about, as these are the ones used as wildcards within the Snakefile.

You can

• Fill in the location of the samplesheet within the config/config.yaml.
• Drop the file in the workdir - Attention: It should be named samples.tsv
• Use snakemake's --config samplesheet=/path/to/my_samples.tsv when executing the wofkflow.

Usage

A dry-run (always a good idea before each execution)

(phap)$ snakemake -n --use-singluarity
--singularity-args "-B /path/to/databases/data:/data"

Basic:

# From within this directory
# Make sure you have defined a samplesheet
(phap)$ snakemake -p --use-singularity -j16 \
--singularity-args "-B /path/to/databases/data:/data"

where /path/to/database/data is the directory containing tables, WIsH models and CRISPR blasts databases.

• The -j flag controls the number of jobs (cores) to be run in parallel. Change this according yo your setup
• The -p flag prints commands that are scheduled for executing. You can
• remove this
• Binding the data dir with the --signularity-args is required (at least in my tests). You must also provide it as a value in the config.yaml.

Output

All output is stored under a results directory within the main workdir. Results are stored per sample according to the sample ids you provided in the sample sheet. For each sample, results for each tool are stored in directories named after the tool. An example looks like this:

$ tree -L2 results/A
results/A
├── all_predictions.tsv
├── lca.tsv 
├── htp
│   ├── predictions.tsv
│   └── raw.txt
├── rafah
│   ├── A_CDS.faa
│   ├── A_CDS.fna
│   ├── A_CDS.gff
│   ├── A_CDSxMMSeqs_Clusters
│   ├── A_Genomes.fasta
│   ├── A_Genome_to_Domain_Score_Min_Score_50-Max_evalue_1e-05.tsv
│   ├── A_Ranger_Model_3_Predictions.tsv
│   ├── A_Seq_Info.tsv
│   └── predictions.tsv
├── tmp
│   ├── filtered.fa.gz 
│   ├── genomes
│   └── reflist.txt
├── vhmnet
│   ├── feature_values
│   ├── predictions
│   ├── predictions.tsv
│   └── tmp
├── vhulk
│   ├── predictions.tsv
│   └── results
└── wish
    ├── llikelihood.matrix
    ├── prediction.list
    └── predictions.tsv

Per sample

---

all_predictions.tsv

Contains the best prediction per contig (rows) for each tool along with its confidence/p-value/whatever-single-value each tool uses to evaluate its confidence in the prediction.

An example for three genomes:

contig  htp_proba       vhulk_pred      vhulk_score     rafah_pred      rafah_score     vhmnet_pred     vhmnet_score    wish_pred       wish_score
NC_005964.2     0.8464285626352002      None    4.068828        Mycoplasma      0.461   Mycoplasma fermentans   0.9953  Bacteria;Tenericutes;Mollicutes;Mycoplasmatales;Mycoplasmataceae;Mycoplasma;Mycoplasma fermentans;Mycoplasma fermentans MF-I2   -1.2085700000000001
NC_015271.1     0.995161392517451       Escherichia_coli        1.0301523       Salmonella      0.495   Muricauda pacifica      0.9968  Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;Raoultella;Raoultella sp. NCTC 9187;Raoultella sp. NCTC 9187       -1.3869200000000002
NC_023719.1     0.9999957241187084      Bacillus        0.0012575098    Bacillus        0.55    Clostridium sp. LS      1.0000  Bacteria;Firmicutes;Clostridia;Clostridiales;Clostridiaceae;Clostridium;Clostridium beijerinckii;Clostridium beijerinckii       -1.29454

lca.tsv

Last Common Ancestor of predictions, based on taxonomy

An example for the genomes above:

contig  name    rank    lca
NC_005964.2	Mycoplasma	genus	2093
NC_015271.1	Enterobacteriaceae	family	543
NC_023719.1	Firmicutes	phylum	1239

tmp (dir)

• Directory genomes: Contains one fasta file per input genome
• File reflist.txt: An intermediate file that holds paths to all produced genome fastas (used as intermediate file to ensure smooth execution)
• File filtered.fa.gz: Fasta files containing sequences > 5000 bp.

Per tool

htp

• File raw.txt: The raw output of htp per contig
• File predictions.tsv: Two-column separated tsv with contig id and probability of host being a phage.

rafah

• Files prefixed with <sample_id>_ are the rafah's raw output
• predictions.tsv: A selection of the 1st (Contig) , 6th (Predicted_Host) and 7th (Predicted_Host_Score) columns from file <sample_id>_Seq_Info.tsv

vhulk

• File results.csv: Copy of the results/sample/tmp/genomes/results/results.csv
• File predictions.tsv: A selection of the 1st (BIN/genome), 10th (final_prediction) 11th (entropy) columns from file results.csv.

vhmnet

• Directories feature_values and predictions are the raw output
• Directory tmp is a temporary dir written by VirHostMatcher-Net for doing its magic.
• File predictions.tsv contains contig, host taxonomy and scores.

wish

• Files llikelihood.matrix and prediction.list are the raw output
• File predictions.tsv has contig, host taxonomy and llikelihood scores.

Logs

Logs capturing stdout and stderr during execution of each rule can be found in workdir/logs/<sample_id>/*.log files.