Pharokka: a fast scalable bacteriophage annotation tool

1As the number of bacteriophage (phage) sequences increases, there is a need for bioinformatic tools that enable fast, consistent and scalable genomic annotation. Existing tools such as RAST, PHASTER and CPT Galaxy are web-server based which may be laborious in particular when multiple phage genomes require annotation. On the other hand, currently available customizable bioinformatics pipelines such as multiPhATE2 require significant understanding of bioinformatics to implement and have lengthy run times. Furthermore, command-line programs designed for viral discovery in metagenomic datasets such as Cenote-Taker 2, Hecatomb and MetaPhage require significant computational resources and storage for database installation.

2As a result, one-line prokaryotic genome annotation tools such as Prokka are often used for phages, especially where tens or hundreds of phages need to be annotated simultaneously. However, such tools implement prokaryotic gene prediction tools that are based on models that are not designed for phage genomes. In addition, phage genes are often lacking in their default databases, resulting in incomplete or hypothetical functional annotation of phage genes.

3Inspired by Prokka, here we created Pharokka as a one-line tool tailored to phages that provides annotations in a fast, scalable and consistent fashion. Pharokka identifies predicted coding sequences (CDS), transfer RNAs (tRNAs), transfer-messenger RNAs (tmRNAs) and clustered regularly interspaced short palindromic repeats (CRISPRs), providing functional annotation for CDS using the PHROGs database. With accessible bioconda installation, Pharokka is easily integrated into complex bioinformatic pipelines such as those created with Snakemake or Nextflow.

4The Pharokka workflow is outlined in Figure 1.

5Pharokka requires assembled DNA sequences in FASTA format (Fig. 1A). For phage isolates, this usually consists of one complete contig, but Pharokka can also annotate incomplete assemblies or metavirome samples with multiple contigs in the multiFASTA format. Furthermore, metagenomically assembled phage genomes and genomic contigs, derived using programs such as Virstorter2, Hecatomb and Cenote-taker 2, are also suitable to be annotated using Pharokka using meta mode.

6Pharokka uses PHANOTATE by default to predict CDS, as it is the only existing tool that is designed to predict genes in phage genomes (; Fig. 1B). PHANOTATE considers unique features of phage genomes such as small gene size, high coding density and alternative start codons. Moreover, PHANOTATE predicts significantly more genes that are on average smaller than other gene prediction tools. Small phage genes are more prevalent than predicted based on existing non-phage specific annotation tools and may encode for anti-CRISPR and antimicrobial resistance proteins. Alternatively, Pharokka users can specify Prodigal, a gene predictor designed for prokaryotic genes, be used for CDS prediction (; Fig. 1B). Prodigal may be useful if users wish to annotate large metavirome datasets quickly using meta mode (Table 4). In addition, Pharokka employs tRNAscan-SE 2 to predict tRNAs, Aragorn to predict tmRNAs and MinCED to predict CRISPRs (; Fig. 1D). When meta mode is specified, Pharokka runs PHANOTATE on each contig separately, leading to considerable speed improvement when multiple threads are specified.

7Functional assignment of predicted CDS is conducted using the PHROGs database (; Fig. 1C). The PHROGs database contains 38 880 PHROGs (protein orthologous groups) containing 868 340 proteins from 17 473 complete genomes of viruses infecting bacteria or archaea that are grouped together using remote homology detection. Each PHROG is assigned to one of nine functional categories. Each gene within each PHROG has a specific product description if known. All predicted CDS are translated into a protein sequence and assigned to the closest matching protein and accompanying PHROG in the PHROGs database using the protein searching tool mmseqs2. An e-value threshold of 10−5 is used by default. Users can also specify their own e-value if desired. This is particularly useful if the phage is novel with few or no similar sequences in the PHROGs database, as a lower threshold can be used to detect less similar matches. If no PHROG match is found, then the CDS is annotated as ‘No PHROG’ and ‘hypothetical protein’. All CDS that are annotated as ‘terminase large subunit’ are extracted into separate output files, as the terminase large subunit is commonly used for phylogenetic analysis.

8While virulence factors are commonly found on prophages, lytic phages rarely encode antibiotic resistance genes. Nonetheless, screening is required for all phages that are intended to be used for phage therapy. Pharokka adds antimicrobial resistance and virulence gene detection using the Comprehensive Antibiotic Resistance Database (CARD) and the Virulence Factor Database (VFDB) (; Fig. 1C). All protein CDS are assigned to the closest matching protein in each database with mmseqs2 if the match passes strict thresholds of 80% identity over 40% coverage recommended by.

9Pharokka’s output files are outlined in Table 1. The primary output of Pharokka is a.gff file that is suitable for use in downstream pan-genome programs such as Roary and panaroo (; Fig. 1E). Other files include a.tbl file, which is a flat-file table suitable to be uploaded to the NCBI's Bankit, a cds_functions.tsv file, which includes counts of CDS, tRNAs, CRISPRs and tmRNAs and CDS within each functional category for each contig in the input FASTA file, a length_gc_cds_density.tsv file, which outputs the length, GC percentage and CDS coding density for each contig, a cds_final_merged_output.tsv, which gives the combined parsed output from mmseqs2 searches against the PHROGs, VFDB and CARD databases, and terL.faa and terL.ffn output files that contain the amino acid and nucleotide sequences of any predicted terminase large subunit genes. When run on metavirome input, Pharokka’s contig-level summary output allows users to identify specific contigs within the metavirome that possess unusual features such as virulence factors, antimicrobial resistance genes or potential stop codon reassignment as indicated by low CDS coding density.

10To test the performance of Pharokka, we compared the run-time and annotations of Enterobacteria phage lambda (Genbank accession J02459), Staphylococcus phage SAOMS1 (Genbank accession MW460250) and 673 crAss-like metagenome-assembled phage genomes from the human gut with default Pharokka v1.1.0 using PHANOTATE as a gene predictor, Pharokka v1.1.0 specifying Prodigal as gene predictor and Prokka v1.14.6 using a version of the PHROGs HMM database that has been reformatted for use with Prokka (https://millardlab.org/2021/11/21/phage-annotation-with-phrogs/). For the 673 crAss-like phage genomes, Pharokka’s meta mode was employed. Benchmarking was conducted on an Intel® Xeon® CPU E5-4610 v2 @ 2.30 GHz specifying 16 threads for Pharokka and 16 cpus for Prokka.

11Tables 2 and 3 show that for phages Lambda and SAOMS1, Pharokka is slower than Prokka but still fast, finishing within 5 min regardless of gene predictor used, with PHANOTATE predicting more CDS with higher coding density than Prodigal. Table 4 shows that for the crAss-like phage genomes, Pharokka is considerably faster than using Prokka regardless of gene predictor, due to Pharokka using mmseqs2 for database searching rather than HMMER3 employed by Prokka. As the size of the input increases, the extra-time cost incurred by Pharokka is due to the gene prediction and tRNA calling, rather than the database searching. For extremely large metavirome datasets, Pharokka in Prodigal meta mode is therefore recommended.

12In addition, for the crAss-like phage genomes, the low coding density output of some contigs identified by Pharokka indicates that stop codon reassignment may be occurring in these contigs.