previous episode

Viromics2026

next episode

Host Prediction II

Overview

Teaching: 0 min
Exercises: 300 min
Objectives

  • Run RaFAH and PhageTransformer to predict a host genus for each contig.

Host prediction

Today, we will use the tools RaFAH and PhageTransformer to predict a genus of a putative host for our contigs. RaFAH predicts proteins in each viral sequence and then assigns them to a set of orthologous groups. The annotated proteins are then used to predict a host genus with a pretrained random forest model.

Next, we will use PhageTransformer as an alternative way of linking viruses to their hosts. PhageTransformer employs a genome language model to predict one or more hosts for a given input sequence. It can also identify bacterial sequences and gives you an additional estimate of how reliable a prediction is.

Finally, we will compare the resulting predictions and discuss the results together.

RaFAH

Exercise - Use RaFAH to predict hosts for our contigs

RaFAH requires a single file for each contig to run. You first have to write a python script which separates the combined assembly into single files. You can use the package biopython installed in your virtual environment for this:

import os, sys
from Bio import SeqIO

# define file paths from the arguments
...

# loop through the records in the combined assembly
with open(assembly_path) as handle:
    for record in SeqIO.parse(handle, "fasta"):
        # set a filename per record
        out_fasta = os.path.join(out_dir, f"{record.id}.fasta")
        # write the record to the file
        with open(out_fasta, "w") as fout:
            SeqIO.write([record], fout, "fasta")

You can run the script in the same sbatch script as RaFAH. Remember to source and deactivate our python virtual environment accordingly. Here you can find a description of the parameters you can pass to RaFAH (the page is a bit hard to read). The tool is programmed in perl and you can find it here on draco:

# activate the conda environment with the dependencies RaFAH requires
source /vast/groups/VEO/tools/miniconda3_2024/etc/profile.d/conda.sh && conda activate perl_v5.32.1

# set a variable to call the RaFAH script
rafah='/home/groups/VEO/tools/rafah/RaFAH.pl'

# create an output folder or use the one you set for the slurm logs:
mkdir 10_results_hostprediction_rafah

# run RaFAH with the appropriate file paths and arguments
perl "$rafah" --predict --genomes_dir 10_results_hostprediction_rafah/split_contigs --extension .fasta

RaFAH is computationally expensive and can use multiple threads. We recommend you set the following sbatch parameters:

  • #SBATCH –cpus-per-task=20
  • #SBATCH –partition=standard
  • #SBATCH –mem=50G

RaFAH uses the random forest model to compute a probability for all host genuses included in its training. You can find these probabilities in the output file *Host_Predictions.tsv. The file *Seq_Info_Prediction.tsv contains per contig the genus with the highest probability.

  1. How many contigs have a high probability score and do you trust these predictions?

python script for splitting the assembly into separate files

import os, sys
from Bio import SeqIO

def main():

    # define file paths from the arguments
    assembly_path = os.path.abspath(sys.argv[1])
    assert assembly_path.endswith(".fasta")

    # set an output directory and create it if it does not exist
    out_dir = os.path.abspath(sys.argv[2])
    if not os.path.exists(out_dir): os.makedirs(out_dir)

    # loop through the records in the combined assembly
    with open(assembly_path) as handle:
        for record in SeqIO.parse(handle, "fasta"):
            # set a filename per record
            out_fasta = os.path.join(out_dir, f"{record.id}.fasta")
            # write the record to the file
            with open(out_fasta, "w") as fout:
                SeqIO.write([record], fout, "fasta")

if __name__ == "__main__":
    main()

sbatch script for host prediction with RaFAH

#!/bin/bash
#SBATCH --tasks=1
#SBATCH --cpus-per-task=32
#SBATCH --partition=standard
#SBATCH --mem=50G
#SBATCH --time=02:00:00
#SBATCH --job-name=rafah
#SBATCH --output=10_rafah/rafah.slurm.%j.out
#SBATCH --error=10_rafah/rafah.slurm.%j.err

assembly='../1.3_virus_identification/30_filter_contigs/assembly.fasta'
contigs='10_rafah/split_contigs'
mkdir -p "$contigs"

# RaFAH expects each genome in a separate file. Activate our virtual environment
# and run a python script to split our filtered assembly into single files
source ../py3env/bin/activate

# The script requires the assmbly path and a directory for outputting the contigs
python ../python_scripts/2.1_split_assembly.py $assembly $contigs

# deactivate the python environment, just to be sure not to cause problems with the conda environment RaFAH needs
deactivate

# activate the conda environment with the dependencies RaFAH requires
source /vast/groups/VEO/tools/miniconda3_2024/etc/profile.d/conda.sh && conda activate perl_v5.32.1

# set a variable to call the RaFAH script
rafah='/home/groups/VEO/tools/rafah/RaFAH.pl'

# rafah parameters (https://gensoft.pasteur.fr/docs/RaFAH/0.3/)
# --predict: run the pipeline for predicting hosts
# --genomes_dir: the directory with the separate files for each contig
# --extension: the extension of the contig files
# --file_prefix: can specify an output dir here and "run name" (rafah_1) here

perl "$rafah" --predict --genomes_dir $contigs --extension .fasta --file_prefix 10_rafah/rafah_1

# deactivate RaFAH's conda environment
conda deactivate

PhageTransformer

Exercise - Use PhageTransformer to predict hosts for our contigs

Unlike RaFAH, PhageTransformer takes the whole assembly as a single multi-fasta file, so no splitting step is needed. It reads the nucleotide sequences directly and does not depend on gene calling or protein annotation.

# activate the conda environment holding the PhageTransformer installation
source /vast/groups/VEO/tools/miniconda3_2024/etc/profile.d/conda.sh && conda activate phagetransformer_v0.2.1

# the pretrained model lives in its own directory on draco
model="/home/groups/VEO/tools/phagetransformer/v0.2.1"

# run the prediction; results are written to standard output, so redirect them into a file
phagetransformer predict --input <assembly> --model_dir $model > <output.tsv>

PhageTransformer is a language model and runs roughly 200x faster on a GPU than on CPUs. We recommend you set the following sbatch parameters:

  • #SBATCH –cpus-per-task=4
  • #SBATCH –partition=short,standard
  • #SBATCH –mem=10G

Since we only have a handful of contigs, the job also finishes in reasonable time on an ordinary CPU node. If the GPU partition is busy, drop the --gres line and use --partition=short,standard instead.

The output table contains one row per contig with the predicted host and a score expressing how confident the model is. PhageTransformer can also report that a sequence looks bacterial rather than viral, which is a useful sanity check on our filtered assembly.

  1. How many contigs get a host prediction, and how are the confidence scores distributed?
  2. Does PhageTransformer flag any of our contigs as bacterial? What would that mean for the filtering we did in the virus identification section?

sbatch script for host prediction with PhageTransformer

#!/bin/bash
#SBATCH --tasks=1
#SBATCH --cpus-per-task=4
#SBATCH --partition=standard
#SBATCH --mem=10G
#SBATCH --time=01:00:00
#SBATCH --job-name=pt
#SBATCH --output=./20_phagetransformer/slurm/phagetransformer.slurm.%j.out
#SBATCH --error=./20_phagetransformer/slurm/phagetransformer.slurm.%j.err

# source the conda env on draco, set model directory
source /vast/groups/VEO/tools/miniconda3_2024/etc/profile.d/conda.sh && conda activate phagetransformer_v0.2.1
model="/home/groups/VEO/tools/phagetransformer/v0.2.1"

# PhageTransformer takes the whole filtered assembly at once
assembly='../1.3_virus_identification/30_filter_contigs/assembly.fasta'
outtsv='20_phagetransformer/predictions.tsv'

# the tool writes its table to standard output, so we redirect it into a file
phagetransformer predict --input $assembly --model_dir $model > $outtsv

conda deactivate

Exercise - comparing RaFAH and PhageTransformer results

  1. For how many host predictions do RaFAH and PhageTransformer agree and upto which taxonomic level?
  2. How do you explain the differences between the predictions?
  3. Which predictions do we trust, how can we get more confident?

Key Points

  • RaFAH uses a random forest model to predict hosts to the genus level for phages

  • RaFAH returns a probability for each host genus it can predict.

  • PhageTransformer is a codon-aware genome language model that predicts a host from DNA sequence directly.

previous episode

next episode