Venomix: a simple bioinformatic pipeline for identifying and characterizing toxin gene candidates from transcriptomic data

Data acquisition and transcriptome assembly

Raw reads from four different analyses were downloaded from the short read archive (SRA) on GenBank (C. sponsalis: SRR260951 (Phuong, Mahardika & Alfaro, 2016); U. yaschenkoi: SRR1557168 (Luna-Ramírez et al., 2015); E. colaratus: ERR216311 –ERR216312 (Hargreaves et al., 2014); T. bicarinatum: SRR1106144 –SRR1106145 (Bouzid et al., 2014)). The previously published transcriptome level analysis for U. yaschenkoi and T. bicarinatum were restricted to only characterizing the venom gland transcriptome in their respective species (Bouzid et al., 2014; Luna-Ramírez et al., 2015). While the C. sponsalis and E. coloratus venom gland transcriptomes were investigated in a comparative framework alongside closely related taxa (Hargreaves et al., 2014; Phuong, Mahardika & Alfaro, 2016). In our study, all four transcriptomes were assembled in Trinity (Grabherr et al., 2011; Haas et al., 2013), using default parameters of its built-in Trimmomatic program to clean up the sequences (Bolger, Lohse & Usadel, 2014). For each transcriptome, expression values as transcripts per kilobase million (TPM) and fragments per kilobase million (FKPM) were calculated in the program RSEM (Li & Dewey, 2011) as part of the Trinity package.

Analysis pipeline and execution

The bioinformatic pipeline for Venomix is outlined in Fig. 1. The program requires three inputs provided by the user: (1) an assembled transcriptome, (2) gene expression information in the form of a tab-delimited output with transcript names in the first column, and (3) tab-delimited BLAST output using the ToxProt as query sequences against the trainscriptome. Although we used Trinity and RSEM in our study, most assembly and expression profile programs with these formats should work with Venomix. Following transcriptome assembly and expression level calculations, the final user provided file is created using tBLASTn from NCBI BLAST+ version 2.4.0 (Camacho et al., 2009), with the ToxProt dataset as the search query with the final BLAST alignment results shown in a tabular format (-outfmt 6). Query sequences from ToxProt are provided within the Venomix package; however, alternative curations of the ToxProt dataset may be used if the sequence identifiers are not changed. In our analysis, we implemented two BLAST search procedures; the first used a more stringent identification algorithm (E-value 1e−20) and a less stringent identification algorithm (E-value 1e−6).

Venomix was tested on the University of North Carolina at Charlotte COPPERHEAD Research Computing Cluster, while requesting computational resources that may mimic most personal laptops/desktops, specifically one processor and 4 GB of RAM. Using these settings the Venomix completed in less than twenty minutes for each of the focal transcriptome. The implementation of Venomix requires the scripting languages Python 2.7 (http://www.python.org/download/releases/27/) and R 3.1.1 (https://cran.r-project.org/bin/), in addition to other Biopython packages (Cock et al., 2009) and data from ToxProt and Genbank that are built into the Venomix pipeline (https://bitbucket.org/JasonMacrander/venomix). We included versions of MAFFT (Katoh & Standley, 2013), NCBI BLAST+ (Camacho et al., 2009), and Transdecoder (Haas et al., 2013) that can be run locally. Although there are two versions of MAFFT (64 bit and 32 bit), the default is the 64-bit, as this is more common for computers used in bioinformatic analyses. Modification to the version of MAFFT in the Venomix pipeline can be done in the support_files/Alignment.py file.

Once the user specifies the input files (transcriptome, expression file, and BLAST output), the Venomix pipeline automatically produces several potentially useful and informative files within each of the toxin group directories (Fig. 1). The outputs within each of the Toxin Group directories are intended to provide the user with curated information to focus future investigations and analyses. Specifically the outputs are as follows: (1) unaligned full and partial nucleotide sequences representing candidate toxins, (2) UniProt query sequences used to identify candidates, (3) taxonomic information associated with query sequences, (4) Genbank entries from query sequences, (5) translated protein sequences from full nucleotide sequences (no partial), (6) UniProt toxin and candidate protein sequences unaligned and (7) in alignment, (8) SignalP information, (9) mature peptides when signaling regions are present, (10) a neighbor joining tree based on the protein sequence alignment, and (11) expression information, if provided by the user. All of these files are contained in directories with representative toxin names retrieved from ToxProt proteins. As the pipeline does not definitively identify candidate toxins to toxin type, similar toxin gene families likely have similar candidates in different directories.

Depending on the next step of the analysis, some of the output files may be used for additional venom related downstream applications or simply a quick reference for the user. Venomix also creates two ancillary products that may be informative to some users: TPM.fasta (only transcripts with TPM values >1.0) and a large GenBank file with information from ToxProt BLAST hits in a format that may lend itself to quick searches or downstream annotation. The user may choose to rerun Venomix with TPM.fasta instead of their assembled transcriptome if they would like to characterize only transcripts with a TPM >1.0, but this is not recommended when looking for rare or extremely diverse toxin genes.

Venomix evaluation

For each assembled transcriptome, we identified candidate toxin genes using the Venomix pipeline using a stringent (E-value 1e−20) and less stringent (E-value 1e−6) search strategy in BLAST. Venomix outputs were compared for both search parameters in terms of the number of toxin groups, number of transcripts, and number of “candidate” transcripts identified by the pipeline. A transcript was considered a “candidate” if the transcript had significantly better E-value associated with a toxin than with a non-toxin protein in Uniprot. Candidate transcripts were translated into their protein sequences using Transdecoder (Haas et al., 2013) and further evaluated in ToxClassifier (Gacesa, Barlow & Long, 2016). If a predicted protein sequence received a score of >1 it was considered a toxin candidate. In addition to screening the overall transcriptome analyses, some toxin groups and candidate toxins identified in our analysis were subjected to additional screening beyond what is included in the Venomix pipeline representative of what some users may do in downstream analyses. Sequence alignments for candidate transcripts shown below were done using MAFFT (v.1.3.3) (Katoh & Standley, 2013) and visualized in Geneious (Kearse et al., 2012). Toxin gene tree reconstructions were done in Fasttree v2 (Price, Dehal & Arkin, 2010) using maximum likelihood tree reconstruction methods and bootstrap supports calculated over 1,000 replicates. For the Bouzid et al. (2014) dataset, Venomix was used to compare alternative assembly approaches (Oases/Velvet vs. Trinity), in addition to assessing both transcriptomes for overall completeness in BUSCO (Simão et al., 2015). Expression values for each transcriptome were calculated using RSEM (Li & Dewey, 2011) rather than raw read counts as originally proposed by Bouzid et al. (2014).

Transcriptome assemblies

Transcriptomes reassembled in Trinity (Grabherr et al., 2011) produced similar de novo assembly outputs when compared to the original studies (Table S1), with the only notable difference being in the number of transcripts for C. sponsalis, which may be due to repetitiveness and sequence complexity encountered during their initial assemblies (Phuong, Mahardika & Alfaro, 2016). The transcriptome for T. bicarinatum was originally done using Velvet/Oases (Li & Durbin, 2009); however, we compared this to our Trinity assembly because of its ease of use (Sanders et al., 2018) and frequency in the venom literature (Macrander, Broe & Daly, 2015), in addition to a lower redundancy and chimera rate (Yang & Smith, 2013).

Pipeline output

In the original published studies, species-specific transcriptomes of C. sponsalis, E. colaratus and U. yaschenkoi were not subjected to any BLAST searches using ToxProt, but instead were screened using taxonomic specific toxin datasets from venom proteins of closely related species. The Venomix pipeline recovered the majority of these lineage specific toxins, as well as additional transcripts that resemble toxin genes from other taxa. It is worth noting that if there were lineage specific toxins that shared high sequence similarity to other toxins from distantly related taxonomic groups, the toxin group name may be assigned an incorrect lineage classification, yet remain a toxin candidate. For example, analyses of the scorpion U. yaschenkoi resulted in four venom groups with “Snake venom” in the name; however, in these instances, lineage specific toxin names are often members of larger gene families that may not be lineage specific toxins. The number of identified toxin groups varied considerably across species and stringency parameters (Table 1), with the less stringent parameters dramatically increasing the number of toxin groups. Each species had multiple toxin groups that were separated based on sequence similarity and correspond to large gene families, including astacin-like metalloproteases, conotoxins, phospholipases A2s, CRISPs, Kunitz-type serine protease inhibitors, snaclecs, metalloproteinases, thrombin-like enzymes, and allergens. For each species, there was at least one toxin group that was not retained following the reciprocal BLASTp search after the toxin-like transcripts were translated into their open reading frame (Table 1).

Venomix outputs for C. sponsalis

Conotoxins represent some of the best-studied toxin genes across the genus Conus, comprising of multiple gene families with cysteine rich proteins (Buczek, Bulaj & Olivera, 2005; Kaas et al., 2012). Venomix identified 76 toxin groups comprising of 246 toxin gene candidates based on our preliminary low stringency BLAST survey; 20 of these groups cluster with various conotoxins (Fig. 2). In total, there were 179 of the 246 toxin gene candidates from the 20 conotoxin groups. There were three instances where Venomix recovered more candidate conotoxins than the original study (O1, conodipine, and conophysin superfamilies), however, a majority of the conotoxin superfamilies identified by Phuong, Mahardika & Alfaro (2016) were missing from our analysis (Fig. 2). This discrepancy is likely due to the different assembly approaches, as iterative assemblies used by Phuong, Mahardika & Alfaro (2016) were unable to recover known transcripts using Trinity alone. Beyond the conotoxins, there were 19 candidate toxin genes found within the Kunitz-type conkuitzin-S1 group, which included characterized toxin proteins from the venom Kunitz-type family of sea anemones, cone snails, and snakes (Supplemental Information 1).

Venomix outputs for E. coloratus

Snakes represent one of the best-studied animal lineages and we were able to validate Venomix’s efficiency with the E. coloratus transcriptome. With a low stringency search (E-value = 1E−6) Venomix identified 132 toxin groups (Table S2), Among the transcripts identified, 45 had a TPM value greater than 100, with 39 of these in the venom gland, four in the scent gland, and two in the skin. The majority of the highly expressed transcripts in the venom gland (TPM > 100) corresponded with toxin groups previously identified (Hargreaves et al., 2014), comprising mostly of C-type lectins, cysteine rich venom proteins, disintegrins, metalloproteinases, and several others (Fig. 3). In addition to these venom candidates, we found one highly expressed cystatin in the venom gland, although it was also highly expressed in other tissues and not thought to be a toxic component of the E. coloratus venom (Hargreaves et al., 2014). We also identified a peroxiredoxin, ficolin, and three latroinsectotoxins (Supplemental Information 2), all of which may have some role during envenomation now well characterized in snakes (Magazanik et al., 1992; Calvete et al., 2009; OmPraba et al., 2010), but require further investigation.

Venomix outputs for T. bicarinatum

Ants represent a poorly-studied lineage of venomous animals, with T. bicarinatum being the only transcriptome used in our study which identified toxins using the ToxProt dataset (Bouzid et al., 2014). The original transcriptome assembly was done in Velvet/Oases (Li & Durbin, 2009). Despite the alternative approaches, the BUSCO (Simão et al., 2015) scores were similar, with the Velvet/Oases assembly at 95.9% and Trinity at 92.2%. When considering TPM (rather than raw counts) the number of candidate genes in the venom gland following the approach by Bouzid et al. (2014) were similar to what was originally published (Table 1). Among these 527 candidates, there were only 44 predicted ORFs from Transdecoder, and only three of these were given a score of 1 or greater in Toxclassifier (Table 1). The BLAST screening, however, resulted in 62 candidate toxins identified when the E-value threshold was set to 1E−3, but 287 when the E-value threshold was set to 10. As E-values were not specified by Bouzid et al. (2014) both are reported here (Table 1). In the less stringent analysis, the largest number candidate toxin genes corresponded to 280 transcripts in the alpha-latroinsectotoxin-Lt1a group, but overall expression within the ant carcass and venom gland across these transcripts were approximately the same. Among those more highly expressed transcripts in the venom gland, six had TPM scores greater than 100 in the venom gland; four corresponding to Venom Allergen 3 and two to cysteine-rich venom protein Mr30. Subsequent BLAST searches against UniProt indicated that all six of these toxins are likely Venom Allergen 3 toxins, making up ∼92% of the cumulative toxin genes expressed within the transcriptome (Supplemental Information 3).

When we compared the Venomix outputs for our Trinity assembly to the Velvet/Oases assembly that was previously published by Bouzid et al. (2014) we recovered some unexpected results. Although we used RSEM instead of BWA, only 33 of the original 69 candidate toxin sequences recovered from their analysis had higher TPM expression in the venom gland than in the ant carcass. Of these 33 potential toxins, only 10 had a 1,000 fold higher expression based on expected count values (Table S3). Such a drastic difference in the reported fold change cross these toxins were unexpected when considering only the mapping program was different. Further examination revealed that thrombin-like enzymes from the Velvet/Oases assembly were likely only partial sequences (Supplemental Information 4), which may have contributed to the larger than 50% discrepancy between our analyses and is likely the reason why only 47 of the 8,688 transcripts identified in Venomix were translated into their open reading frame.

Venomix outputs for U. yaschenkoi

Scorpions represent one of the oldest and best-studied venomous lineages, with the original study by (Luna-Ramírez et al., 2015) using scorpion-specific toxins as query sequences. Ultimately this original study identified 210 transcripts representing 111 unique scorpion toxins, venom gland enzymes, and antimicrobial peptides (Luna-Ramírez et al., 2015). By expanding the query sequences with the ToxProt dataset, we recovered 117 toxin groups representing 689 unique transcripts (Table 1). Of the 117 identified toxin groups only 10 were from scorpion derived query sequences. Of these 10, the Toxin-like protein 14 had the same number of candidate toxins as the original study, with Venomix recovering the complete protein sequence, when Luna-Ramírez et al. (2015) did not (Fig. 4A). When using exclusively scorpion venom proteins from ToxProt as query sequences, the number of candidate toxins identified by Venomix was approximately the same as that identified by Luna-Ramírez et al. (2015).

The most abundant toxin-like transcripts within the less stringent search found 190 transcripts within the delta-latroinsectotoxin-Lt1a group, the 182 transcripts alpha latroinsectotoxin-Lt1a group, and 91 transcripts within the Neprilysin-1 group (Supplemental Information 5). Query toxins which are used to form these toxin groups were previously identified in spiders (Graudins et al., 2012; Garb & Hayashi, 2013; Undheim et al., 2013; Bhere et al., 2014), and not included in the original transcriptome analysis (Luna-Ramírez et al., 2015). Preliminary screening based on reciprocal BLAST hits indicated that the latroinsectotoxin groups identified in Venomix might not be toxins. Maximum likelihood gene tree reconstructions were used as post-processing steps to further screen the Neprilysin-1 group as potential toxin sequences revealing that candidate toxins from the Neprilysin-1 group formed a well-supported cluster with neprilysin toxins from other scorpions at high expression levels (Fig. 4B).