Rephine.r: a pipeline for correcting gene calls and clusters to improve phage pangenomes and phylogenies

Overview of the pipeline

The Rephine.r pipeline (summarized in Fig. 1) assumes the researcher has already completed a workflow for predicting gene clusters in a pangenome, such as the combination of blastp (Altschul et al., 1990) and MCL (Enright, Van Dongen & Ouzounis, 2002) implemented by Anvi’o (Eren et al., 2021) and other programs (e.g., vConTACT Bolduc et al., 2017 and Roary Page et al., 2015). In what follows, we use Anvi’o as the basis for initial pangenomes, as Anvi’o is both a popular tool for bacterial pangenomes and includes several useful commands for facilitating our corrections. Future updates will expand Rephine.r’s compatibility with other tools.

Following initial gene clustering Rephine.r pipeline: (1) identifies and merges gene clusters containing distantly related homologs using HMMs, and (2) identifies fragmented gene calls that can be fused for the purpose of SCG inference and generating phylogenies. By default, Rephine.r will first run the cluster merging and defragmentation steps in tandem, produce a set of new clusters that combine the results of these corrections, and will then run a second round of defragmentation to identify any new cases that emerge due to the prior steps. Command line options are also offered for users that wish to run the HMM merging or fragment fusion steps individually. In addition to the main pipeline, we include two complementary scripts: getSCG.r returns the single-copy core genes and a concatenated alignment file for phylogenetics; fragclass.r categorizes the likely events that led to fragmented gene calls.

Merging gene families with HMMs

Gene clustering based on sequence similarity relies on threshold criteria for defining when two sequences are related and for clustering related sequences into groups. In Anvi’o, the default identity heuristic is defined by the “minbit” score, the ratio of the BLAST bit score between two sequences and the minimum bit score from blasting each sequence against itself. This metric generally performs well, and for bacteria, where homologs are typically over 50% identical, it is especially successful. For phages, however, this approach can miss more distant homologs. Even using a 35% amino acid identity threshold (Cresawn et al., 2011; Shapiro & Putonti, 2018), we may miss cases that only appear related when viewing alignments or comparing phage genes by structure or synteny. Unfortunately, it is not as simple as specifying a lower minbit threshold, since doing so will also increase the number of unrelated genes that are clustered together erroneously.

Given the initial gene clusters returned by Anvi’o, Rephine.r builds separate HMM profiles for each cluster using the hmmbuild function from HMMER (Eddy, 1998) and converts the concatenated HMM profiles into a database with hmmpress. The script then uses hmmscan to compare every original gene call against each HMM profile. This step is expected to be more sensitive for recognizing distant homologs than the initial blastp, as the HMM profiles make use of variation from multiple members of the same cluster. Significant hits are then defined as follows: for each original gene cluster, the “minimum self-bit” (or “selfbit”) score is recorded as the minimum of the bit scores for each of the gene calls that was initially assigned to that cluster by MCL. This selfbit score then serves as a profile-specific significance threshold. Any gene call that was originally assigned to another cluster but has a bit score greater than this value is then used to establish a putative connection between gene clusters. We also include the option of specifying an absolute minimum bit score as an additional criterion. These connections are recorded in the form of a network edgelist linking gene calls to gene clusters. Next, this edgelist is relabeled to define edges between the original gene clusters that share putative homologs. Finally, this edgelist is used to generate a network with the R (R Core Team, 2013) package igraph (Csardi et al. , 2006), and the connected components are returned with the function “components”. The result defines sets of the original gene clusters that are suitable for merging into a single, larger cluster.

Identifying fragmented gene calls

To find fragmented genes, Rephine.r first identifies every gene cluster that includes at least two sequences from the same genome. These sequences may represent true duplicates or paralogs, or they may be separate pieces of the same original sequence that have been split by one of several processes, including: a frameshift due to an indel, insertion of a selfish genetic element, or being artificially split across the ends of the genome when it was reported to GenBank. This third case may also arise as an artifact of the two other mechanisms. For any of these scenarios, the two pieces of the gene will be notable in two ways: (1) they will align with separate parts of the gene in a multiple sequence alignment, with one piece corresponding to an N-terminal fragment, and the other to the C-terminus; (2) they should have lower sequence similarity to each other than to the average comparison with other sequences in the multiple sequence alignment. Fig. 2A illustrates how a fragmented gene may appear in an alignment.

Given clusters with potential fragments, every gene call within an affected cluster is compared using blastp to every other gene call in the same cluster. For the two focal gene calls from a potential fragmented gene, the bit score from their blast alignment is compared to the mean bit score for other blast results within the gene cluster. We defined the ratio of this pairwise blast to the cluster average as the “relative bit” (or “relbit”). Mathematically, for potential fragments A and B within a gene cluster G, this is defined as: (1)${\displaystyle relbit\left(A,B\right)=\frac{bit\left(A,B\right)}{\overline{bit}\left(A,G\right)},relbit\left(B,A\right)=\frac{bit\left(B,A\right)}{\overline{bit}\left(B,G\right)}}$ where the overbar refers to the mean. The maximum of these relbit values is then used as a criterion for judging similarity between A and B. If this value is below a chosen threshold, the ORFs are considered to be sufficiently dissimilar.

Rephine.r also compares the extent of overlap within the pairwise alignment space between each potential paralog. This step is needed, because dissimilar gene fragments may still have overlaps in the alignment due to alignment errors or if the original fragmentation event was caused by a short duplication. To quantify this overlap, the “percent overlap” is calculated as the size of the ORFs’ intersection within the alignment divided by the number of unique, aligned positions between the two sequences. In mathematical terms, for a gene with potential fragments A and B, we define: (2)${\displaystyle PercentOverlap=\frac{\left|A\cap B\right|}{\left|A\cup B\right|}}$ where the size terms are based solely on the aligned positions within the multiple sequence alignment.

Ultimately, sequence pairs with low relative bit scores (“relbit”) and low percent overlaps (“percoverlap”) are the likeliest to fit our expectations of a fragmented gene call. In practice, we implemented default parameters for these criteria of 0.25 for “relbit” and 0.25 for “percoverlap.” These choices are based on plotting values of each parameter (Fig. S1) from the test cases described below and identifying a set of points that weakly cluster together in the graph. When checked manually, each of these genes appeared to correspond to fragmented calls, whereas nearby points in the graph included potential errors. These parameters can be adjusted at the command line, and we would encourage others to visually inspect their alignments.

Once fragmented genes are identified, a new FASTA file is created in which the original pieces of the full-length gene are artificially spliced (or “fused”) into a single gene call. To preserve the original event that separated the sequences, the script inserts an “X” between the two pieces of the gene. New alignments are then made with MUSCLE (Edgar, 2004) for each affected gene cluster, with these X’s imposing a gap in the alignment (see Fig. 2B for an illustration of this step). If desired, the user can then use the additional script, getSCG.r, to return a list of the single-copy core gene clusters, along with a concatenated alignment file that is suitable for phylogenetics. The script, fragclass.r, can also be used to obtain a table summarizing predicted causes for each type of fragment based on the separation between the original gene calls.

Virus genomic data

Phages in the subfamily Studiervirinae (family Autographiviridae), the subfamily Tevenvirinae (family Myoviridae), and the genus Pbunavirus (family Myoviridae) were chosen as well-studied examples for testing Rephine.r. We downloaded all available RefSeq genomes from each of these taxa from the National Center for Biotechnology Information’s (NCBI) genome browser (as of February 2021). This data set included 145 Studierviruses, 127 Tevenviruses, and 38 Pbunaviruses (a full list of accessions is included in Table S1). The Studiervirinae (e.g., phages T3 and T7) and the Tevenvirinae (e.g., phage T4) are among the best-studied phage subfamilies and include characterized examples of introns and homing endonucleases (Chu et al., 1986; Belle, Landthaler & Shub, 2002; Bonocora & Shub, 2004; Petrov, Ratnayaka & Karam, 2010). These features made these two subfamilies ideal for testing methods for identifying distant homologs and fragmented gene calls. The Pbunaviruses were chosen due to the relatively large number of available genomes at the genus level, offering a less diverse contrast to the other phage groups.

Initial pangenome workflow with Anvi’o

We built an initial pangenome for each phage group using Anvi’o v6.2 (Eren et al., 2021) following the standard pangenomics workflow (https://merenlab.org/2016/11/08/pangenomics-v2/), which uses Prodigal (Hyatt et al., 2010) for gene calling. Alternative gene callers can also be used, and these gene calls can be imported into Anvi’o as part of the anvi-gen-contigs-database program. The “–use-ncbi-blast” flag was specified for the anvi-pan-genome command. Due to the large genetic diversity of phages, we set the minbit threshold to 0.35, based on prior work (Cresawn et al., 2011; Shapiro & Putonti, 2018).

Phylogenetics

Maximum likelihood phylogenies were estimated using IQTREE v2.0.3 (Nguyen et al., 2015) with ModelFinder (Kalyaanamoorthy et al., 2017) to automate choosing the optimal substitution model for each tree. For each of the three virus groups, trees were built based on concatenated alignments for the original SCGs and again following Rephine.r using the expanded SCGs. Tree summary statistics were computed in R using the ape package (Paradis, Claude & Strimmer, 2004) and drawn using ggtree (Yu et al., 2017).

Code Availability

All code for this work is provided on GitHub (https://github.com/coevoeco/Rephine.r).The code includes a walkthrough for running Rephine.r following a standard Anvi’o workflow, as well as utility scripts, getSCG.r and fragclass.r, that provide additional output of the SCG genes and predicted causes of fragmentation events.