New alignment-based sequence extraction software (ALiBaSeq) and its utility for deep level phylogenetics

Algorithm and implementation

The workflow is shown in Fig. 1 and the terminology is outlined in Table 1. The input alignment table can be generated by a number of programs, with supported formats including BLAST (blastn, blastp, blastx, tblastn, tblastx), 15-, 18-, and 22-field HMMER (Eddy, 2011) formats, Phyluce-style LASTZ output format, as well as SAM/BAM alignment formats (refer to online manual for more information). The software parses the input table and groups all HSPs for each contig (hit). For each query-hit pair, HSP hit regions are merged when overlapping and joined in the case of no overlap, producing a pseudocontig. Contig regions derived from HSPs with overlap only in the bait sequence (different contig regions matching the same bait region) are separated into alternative “pseudocontigs”. This accounts for a case of a contiguous assembly where several homologous genes occupy the same chromosomal contig and could be located nearby. Thus, no preliminary contig splitting (for example, as done in Phyluce) is needed for genomic data. For the translated searches, the hits derived from the opposite strands are kept as separate “pseudocontigs”. By default, each particular region of every contig is allowed to match only one bait (to prevent the same sequence from being assigned to multiple baits).

An optional reciprocal best hit (RBH) check can ensure that each pseudocontig which matches a particular bait also matches the bait-matching contig in the reference assembly. This check removes out-paralogs and is also a part of Orthograph, FortyTwo, and Assexon, although is not implemented in most popular target-enrichment bioinformatics pipelines (e.g., aTRAM, HybPiper, and Phyluce). A strict RBH check however may represent a too conservative solution, especially when a DNA-based reciprocal search was performed on divergent taxa. The algorithm may result in too many false negatives, the truly homologous contigs that did not get any significant matches to a reference taxon. Thus, we added a relaxed option of retaining the contigs without matches to the reference.

After verifying that all pseudocontigs match best to the bait they were initially assigned to, they are ranked by alignment scores. Bitscore, e-value, and identity are used for ranking and the user has the option to specify in which order these metrics are used (in case of a tie, the next one will be used to resolve the order). We additionally provide the option (1) to use both bitscore and e-value and consider identity in case of ties; (2) to rescale the metrics by HSP length; and (3) to pick several suboptimal hits with a specifiable threshold (for more information please see the manual). Optional contig stitching is then performed by attempting to iteratively fill all regions of the bait by non-overlapping pseudocontigs. Only pseudocontigs derived from different contigs are allowed to be stitched. Stitched supercontigs receive combined scores of the included contigs to prevent longer incorrect matches ranking higher than the composite correct match. A similar operation is done by HMMER when combining the protein domain scores into an overall hit score. While potentially allowing for chimeric sequence creation, contig stitching greatly improves coverage of long, multi-exon loci mined out of discontiguous assemblies. Apart from ALiBaSeq, a contig stitcher is also implemented in Orthograph, while aTRAM and HybPiper are among the only read-based tools to perform stitching of shorter assembled contigs.

The resulting supercontigs are then ranked by scores again, and the best supercontig or several high scoring supercontigs are then selected. We also implemented an option to select several competing supercontigs with similar scores matching the same bait in case the user wants to perform a more thorough tree-based orthology prediction. The procedure outlined thus far runs without access to an assembly FASTA file and is solely based on the alignment tables. We currently support BLAST, HMMER, LASTZ searches, as well as read aligners with SAM/BAM outputs, with a helper script enabling a tblastn-like search for HMMER.

After the list of matching sequences is prepared, the target assembly file is read, and sequences are extracted. Several extraction modes are offered: whole contigs (no contig stitching allowed), only single best matching region, all matching regions, and all matching regions plus the sequence in between (typically to get introns and variable regions as well). Additionally, extraction of flanking regions is permitted, although no similarity check is implemented for those regions. By default, the output sequences are grouped by locus and appended to existing alignments, but they may also be grouped by the sample from which they originated, or can all be written into a single file. Output sequences can also be translated if a protein-based search is used and only matching regions are retrieved. Additionally, extensive logs and bait-sample correspondence tables with contig names, combined scores, and coordinates are saved.

The main features of the software algorithm and implementation are: assembly-based (reads are not required as input, thus all types of assemblies can be used as input, including many assemblies already available on NCBI which were based on legacy sequencing methods); lack of need for sequence annotations (local alignment search is used to determine homology and find target regions, only matching regions, for example, exons, can be output, with non-matching regions, for example, introns, not being reported); RBH check (sample contigs found are checked for match to contigs in a reference assembly from which baits were derived or to which baits are known to be homologous); contig stitching (genes broken up into multiple small contigs due to low coverage or low complexity intronic regions can be pieced together); universality (ALiBaSeq can work with inputs produced by various preceding assembly and local alignment search steps, and the output of ALiBaSeq is suitable for various subsequent multiple sequence alignment, trimming, and phylogenetic reconstruction steps).

Samples

In order to evaluate the software’s performance, namely amount of data retrieved, proportion of false positives, and speed, at deep phylogenetic levels, we used the two species Rhodnius prolixus and Cimex lectularius. They belong to the two distantly related families, Reduviidae and Cimicidae, respectively, within the heteropteran infraorder Cimicomorpha (Hemiptera). The divergence between the two species was recently dated to the mid-Triassic, roughly 225 million years ago (Johnson et al., 2018). Since genomes and annotations are available for both (Mesquita et al., 2015; Rosenfeld et al., 2016), we tested the retrieval of Rhodnius sequences using divergent Cimex bait sequences for orthologs known to be shared by both species.

We downloaded the original assemblies of the Cimex and Rhodnius genomes (ClecH1 and RproC3) and transcripts/gene sets (ClecH1.3 and RproC3.3) from VectorBase (Giraldo-Calderón et al., 2015). To test and compare both read-based and assembly-based software, we simulated Illumina reads derived from the Rhodnius assembly using Art (Huang et al., 2011). Six simulations at various depths (1×, 3×, 5×, 10×, 20×, 40×) were performed, with uniform coverage across target and off-target regions. We used recently published results (Kieran et al., 2019) of empirical enrichment of a closely related species, Rhodnius robustus, as a target capture sample for the UCE dataset (see below). To review assembly-based software equally, we processed all reads in a similar way, that is, with Clumpify (Bushnell, 2014) and Trimmomatic (Bolger, Lohse & Usadel, 2014) and assembled them using SPAdes (Bankevich et al., 2012). Properties of all used assemblies were evaluated using QUAST (Gurevich et al., 2013) and are shown in Table 2. Whenever non-default settings were used, they were indicated in Table S2.

Baits

To ensure our software was applicable for two commonly-used data types, two sets of loci were used as baits, the Ultra Conserved Elements (UCE) Hemiptera 2.7K v1 set (Faircloth, 2017) and the OrthoDB single copy orthologs for Hemiptera (ODB SCO). UCE loci and ODB SCO orthogroups used are listed in Text S3.

For UCEs, we extracted loci shared between Cimex and Rhodnius from the Hemiptera UCE 2.7K v1 set (Faircloth, 2017) and merged probes to obtain the overall bait regions (since typically two probes tile a 160 bp region). Since 99.9% of the kit was shown to target exons of protein coding genes in the taxon with the most well-annotated genome (Kieran et al., 2019), we prepared DNA and AA sequences for both nucleotide and protein modes of the programs we were assessing. A few particularities of the probe set caused us to perform additional filtering and processing: (1) some introns flanking UCE exons are part of the probes despite their relative high rate of mutation; (2) several UCE “loci” in Cimex and/or Rhodnius are overlapping for about a third of their length, with probes targeting opposite strands (likely due to intron drop out in either of the taxa compared to the other probe set taxa); (3) the bait region borders (probe start/stop coordinates) are not always homologous between Cimex and Rhodnius (because of issue 1); (4) some UCE loci are too close to each other on the scaffold to exclude the possibility of capturing both loci with probes derived from only one (we did not address this issue, since it did not impede our evaluation analyses). Thus, to obtain homologous protein coding regions of the probe set, we first overlapped Cimex probes of the same region and aligned the resulting sequences using BLAST to Cimex genomic scaffolds. We randomly removed one of the UCE loci which were overlapping since these duplicate UCEs could be flagged by some of the software we were testing. We compared the mapping coordinates with the assembly annotation and truncated the probe regions to stay within the CDS. We then compared the CDS with the Cimex proteome to find the start and stop of the ORF and get the protein sequence of the probes. Using protein BLAST, we aligned Cimex AA probes to Rhodnius DNA probes and truncated the latter to make sure that homologous Rhodnius regions are used for evaluations. The resulting UCE dataset had 2,037 loci, with locus lengths ranging from 60 to 160 bp, nucleotide pairwise sequence distances (n-psd) ranging from 1.59% to 50.69% (average 22.44%) and 0% to 60.42% (average 6.88%) on the protein level.

For ODB SCO, we retrieved orthologs that are single copy for Hemiptera from OrthoDB v9. We then extracted Rhodnius and Cimex protein sequences, aligned them with MAFFT (Katoh & Standley, 2013), and trimmed off lateral gapped regions. We subsetted the alignments to retain only those with >100 AA and less than 10% gaps between the two species, then used the resulting truncated protein sequences to obtain the corresponding CDS from the VectorBase sourced transcriptomes of the two species. We obtained 1,027 loci and initially used those for benchmarking. However, when investigating false positive results, we discovered that some results seemed to be correct based on the pairwise distance and were in fact correct according to the VectorBase orthology predictions. We thus compared orthology assignment of 1,027 selected loci between OrthoDB and VectorBase. We found cases where one (primarily Rhodnius-sourced) gene was assigned to an incorrect orthogroup on OrthoDB, as well as cases where a gene was not single-copy (had paralogs in Cimex and/or Rhodnius) according to VectorBase, while being listed as single-copy on OrthoDB. A total of 43 such loci were found and discarded from the evaluations. We additionally inspected the Rhodnius-Cimex alignments and discarded 10 other loci that had poor alignment between the two species, suspecting them of being not true orthologs. The resulting ODB SCO dataset had a total of 974 loci, with locus lengths ranging from 216 to 20,862 bp, and nucleotide pairwise sequence distances ranging from 15.22% to 52.89% (average 29.84%) and 0% to 61.62% (average 20.15%) on the protein level.

Software and settings

For read-based approaches, we tested four software packages, Assexon, aTRAM, HybPiper, and Kollector (Table 3). In order to make the comparison equal, we used processed reads for the SPAdes assemblies (see above) as an input for the read-based software. The default settings of the evaluated software were drastically different, and in some cases inappropriate for the deep level divergence we were looking at, causing suboptimal performance. We thus attempted to standardize them and set parameters of the programs to a very generous minimal level (Table S2). Due to specific requirements of Assexon on NGS read format, we could not run it on our simulated reads. aTRAM’s performance with the SPAdes assembler was poor due to coverage cut-off values and timeout problems; when the time parameter was increased, the analyses took days to complete; thus we ran the protein-based analyses with Velvet (Zerbino & Birney, 2008) for evaluation purposes. HybPiper’s protein-based search speed on the ODB dataset, particularly when using 20× and 40× reads, was extremely slow (did not complete in 8 days), thus we set the filtering settings to default for these trials to maximize locus recovery. We note that changing the HybPiper settings to default did not change the results of the low depth trials (1x–5x), while improving the results on the 10x dataset.

For assembly-based approaches, we tested Assexon, FortyTwo, Phyluce, and our software ALiBaSeq (Table 3). We note that Orthograph is very close in capabilities to FortyTwo, with both programs being primarily designed to work with transcripts. However, Orthograph has a substantial list of dependencies, and requires a relational database, thus, we have not included it in our evaluation. Despite the program BUSCO’s capability to process genomic data and retrieve target sequences with high precision, this software is hard to customize. The lowest level relevant gene set available are orthologs at the level of Insecta. The presence of Rhodnius in the profiles would not allow us to fully test the performance with a divergent reference. We thus did not include BUSCO in our comparison. Settings used are listed in the Table S2. For ALiBaSeq we ran the software both with strict and a relaxed RBH check (see implementation). Figures represent the best performing combination of such analyses (strict for 1-3x assemblies, relaxed for others), for details see Tables S3 and S4. For Phyluce, we generally used guidelines in tutorial 3 (https://phyluce.readthedocs.io/en/latest/tutorial-three.html). We had to use the 0 bp flank option to maximize the locus recovery since many loci were being discarded due to co-location on the same contig. The capture sample was processed directly starting with the SPAdes assembly, bypassing the typical read processing and Trinity assembly for comparison purposes. The phyluce_assembly_match_contigs_to_probes command for the capture sample was run with the same settings as for genomic samples.

All programs were run on the UCR High Performance Compute Center cluster. We capped all read- and assembly-based programs to run maximally with 32 threads to make the performance comparison fair. To assess accuracy of the recovery, obtained DNA (or protein, in the case of FortyTwo) sequences were aligned back to the Rhodnius sequences that were homologous to the Cimex bait sequences. The alignment was done with MAFFT G-INS-i or E-INS-i (Katoh & Standley, 2013); and the—adjustdirectionaccurately option was used for outputs of Kollector, since Kollector does not maintain hit sequence direction with respect to the query sequence. Since only ALiBaSeq has the capability to report the source contig names and coordinates for the retrieved sequences, we had to rely on sequence similarity to determine whether an ortholog or a paralog was recovered. We observed that sequence similarity of less than 100% was too strict a criterion for assessing false positive results. Small alterations in simulated reads, assembly of such reads, incorrect flanking bases of the exon-intron junctions, and alignment errors lower the apparent sequence similarity when the retrieved sequence is compared to the original sequence. Thus, we used a threshold of 10% difference to track strongly false positive results. Since our empirical UCE capture sample belonged to a different species (R. robustus), we used a larger threshold of 20% to account for additional interspecific variation.

We also evaluated software performance for difficult loci, which either have multiple introns, or are represented by multiple contigs in the assembly (for assembly-based tools only). The evaluation of the performance on multi-intron and multi-contig loci was evaluated on the 40x ODB dataset. To assess the impact of introns on the software performance, we queried the annotation file for Rhodnius to obtain the number of introns per SCO, and then regressed the SCO coverage produced by each software on the number of introns using a linear model. To evaluate the impact of SCO fragmentation in the assembly, we performed a blastn search for Rhodnius bait homologs in the Rhodnius 40x read assembly, processed the results with ALiBaSeq to estimate the number of contigs representing each SCO locus, and then regressed the SCO coverage produced by each software on the number of contigs.

Additional evaluation methods

Methods of the software evaluation using the plant dataset are available in the Text S1.

The UCE dataset

Number of loci recovered (Fig. 2A)

On the original assembly, ALiBaSeq recovered up to 2,032 out of 2,037 UCE loci. Protein based approaches (tblastn, tblastx, and phmmer) performed best, with relaxed nucleotide-based searches (blastn with wc9 and discontinuous megablast) performing almost as well. Both Assexon and Phyluce recovered 1,436 and 1,577 loci respectively, while FortyTwo only found ~300 loci. On the assemblies of simulated reads, ALiBaSeq generally performed best, or as well as other assembly-based software. ALiBaSeq, Phyluce, and FortyTwo reached their maximum performance at assemblies of 10x and deeper and also performed well on the empirical target capture data. Read-based methods using nucleotide space searches performed very poorly, with the exception of Assexon (~1,500 out of 2,037 loci), which we could only test on the capture sample. Protein-based aTRAM performed the best of these methods in finding the maximum number of loci on the simulated genomic reads at 5x coverage and less. Protein-based HybPiper only performed well on 40x and target capture data, with only ~1,300 out of 2,037 loci on the 20x data, and much less on the shallower data. Interestingly, the genomic pipeline of Phyluce performed better (1,939 loci) on the target capture data of Kieran et al. (2019) than the standard Phyluce capture pipeline (1,471 loci), a result closer to that reported in the original study (Kieran et al., 2019) for this sample (1,508 loci).

Time (Fig. 3)

Depending on the search type, the ALiBaSeq pipeline on the original assembly completed a forward BLAST search in between 2 s (dc-megablast) and 4 min (tblastx) and a reciprocal search in between 8 and 90 min. HMMER searches ran much longer (up to 7 h). The main script ran for 15 s to 2 min. The total processing time for ALiBaSeq thus was between 8 min and 1.5 h with BLAST searches, and up to 8 h with HMMER searches. Assexon and Phyluce worked faster than ALiBaSeq searches (30 s–2 min total time), while FortyTwo was slower (several days). The same relationship was generally true for the searches on simulated read assemblies. On smaller datasets (1–5x), ALiBaSeq with nucleotide blast was as fast as Phyluce and about as fast as Assexon with protein blast. DNA-based read-based tools worked very quickly, but that can be explained partially by much lower locus recovery. Adding the assembly time to the analysis time, assembly-based tools’ speed (1–1.5 days) is higher on the larger datasets than that of protein-based aTRAM (20x and 40x with 4–10 days). The processing time of protein-based HybPiper on the UCE dataset (9 min–5.5 h) is smaller than the total time of assembly-based tools (2.25 h minimum). Interestingly, the SPAdes assembly time of the capture sample was much lower (~7 min) despite the comparable read number with 1x whole genome data, likely due to lower K-mer diversity. This gave the assembly-based tools a clear advantage in speed (7–17 min total time) over the read-based tools that found a considerable amount of loci (Assexon, protein-based HybPiper, and aTRAM with 19 min, 47 min, and several days respectively).

The ODB SCO dataset

Number of loci (Fig. 4A)

ALiBaSeq found the highest number of loci (up to 933 out of 974) on the original assembly by a very large margin, twice as many as the best alternative (404 out of 974 for Phyluce). On the assemblies of simulated reads, FortyTwo found most loci in the low coverage assemblies (1–3x). ALiBaSeq performed well on low coverage assemblies and the best of all assembly-based tools on assemblies of 10x and deeper. As in the UCE dataset, only protein-based versions of the read-based tools found any significant number of loci. The software aTRAM outperformed HybPiper on low coverage data (1–10x) but recovered fewer loci than HybPiper when analyzing 20x and 40x reads.

Time (Fig. 5)

Processing time on this dataset was longer compared to the UCE dataset, likely due to the larger size of the loci. However, relative performance of the tools was generally similar between the two datasets. The notable exception was HybPiper, with its exonerate (Slater & Birney, 2005) steps taking much longer on 20x and 40x data leading to a slower performance compared to aTRAM on the same data.

Number of introns vs coverage (Fig. 6)

ALiBaSeq in conjunction with HMMER and protein-based BLAST searches retrieved the largest regions of the multi-intron genes among all other assembly-based software. Phyluce and Assexon generally failed to recover genes with over 10 introns (N = 38), while FortyTwo retrieved less than 25% of the sequence of the genes with over 10 introns. Among the read-based software, protein-based aTRAM and HybPiper recovered the largest regions of the multi-intron genes, slightly outperforming ALiBaSeq.

Number of contigs vs coverage (Fig. 7)

Among the assembly-based software only ALiBaSeq performed well in finding and stitching contigs comprising multi-contig genes. Assexon barely reported any multi-contig genes, while FortyTwo and Phyluce output highly incomplete sequences for loci represented by three or more contigs in the assembly.

Additional results

Results of the software evaluation using the plant dataset are available in the Text S1.