Sensommatic implements BLAST (Altschul et al., 1990) and AUGUSTUS (Stanke et al., 2004) to mine and predict sensory receptor genes from whole genome assemblies. It requires two input files: the genome assembly for a target species and a reference file of homologous receptors from closely related species (Fig. 1). We provide five reference files for each class of vertebrates, including mammals, amphibians, reptiles, birds, and fish. The pipeline outputs a single fasta file containing gene predictions for all sensory receptor families, including olfactory receptors (ORs), type I and type II taste receptors (TAS1Rs and TAS2Rs), type I and type II vomeronasal receptors (VN1Rs and VN2Rs), trace amine associated receptors (TAARS), and the visual opsins (OPNs and RHO).
2.1. Generation of a ‘pseudogenome’ to improve run times
To improve run times for poorly assembled genomes, Sensommatic initially performs a blastn search of all reference receptors against the genome assembly. Contigs with significant hits are written out to a ‘pseudogenome’ file. To further optimise speed, users can opt to implement the ‘esl-sfetch’ tool from the HMMER easel mini-apps suite (Eddy, 2011). If specified, Sensommatic utilises the esl-sfetch genome indexing tool to efficiently access and write out contigs.
2.2. Mining sensory receptor genes with BLAST
Using the reference file, blastn searches are performed on each contig. Hits are arranged in order of increasing e-value and decreasing identity. For each ordered hit, the identified region is fed into AUGUSTUS for gene prediction. The range of nucleotides added to the 3’ and 5’ ends of each hit is specified by the user. To account for the larger span that multi-exon genes occupy, Sensommatic uses a default range of 50 Kb for multi-exon genes. The operational default range for single-exon genes is 2500bp. The complementary 5’-3’ sequence is determined for hits on the reverse strand.
2.3. Gene prediction with AUGUSTUS
Hits are fed into AUGUSTUS for gene prediction. AUGUSTUS is trained on various vertebrate species. Users should specify the most closely related species available to improve prediction quality. Sensommatic applies the blat2hints feature from AUGUSTUS to generate hints from the reference file to aid gene prediction. Only the top five BLAST hits from the reference file are used to generate hints for the opsin and multi-exon taste receptors. Including more reference genes reduced the prediction accuracy for these specific gene families. The ‘intronless’ AUGUSTUS gene model is applied when predicting single exon sensory receptor genes.
2.4. Predicting single exon genes
If AUGUSTUS fails to predict a single exon gene, Sensommatic searches for the appropriate open reading frame using tblastn. The prediction is checked for start and stop codons at the 5’ and 3’ ends, respectively. If undetected, Sensommatic extends the sequence in-frame to search for the start and end positions (Supplementary Figure S1). If identified, a 'blastback’ verification step is carried out to compare the prediction to the reference gene. Predictions which share over 60% cover and 50% identity with the reference gene, post extension, are considered valid. If these thresholds are not met, or the start and end positions are not identified within the maximum extension limit, an internal search for the start and stop codons is performed (Supplementary Figure S2). An additional 'blastback’ verification step is then carried out. If the prediction fails this final checkpoint, the sequence is reverted back to the original BLAST coordinates and included for downstream analysis.
2.5. Validation of predictions with HMMER
Sensommatic implements HMMER to scan predictions for the presence of the characteristic seven transmembrane domains using profile hidden markov models (HMMs) for Class A and Class C GPCR receptors (Eddy, 2011). Predictions lacking transmembrane domains are filtered out to prevent the inclusion of off-target hits that may have been picked up during the initial BLAST searches.
2.6. Classification of sensory receptor genes
Predictions are classified as ‘functional’ or ‘pseudogene’ based on the presence of in-frame stop codons and prediction length. By default, predictions with lengths shorter than 860bp are considered non-functional. This threshold was chosen based on the length of the seven- transmembrane domain. Users can alter this threshold based on their specific requirements. A shorter length may be more appropriate for low-quality genome assemblies where genes are fragmented by contig boundaries. Predictions are assigned gene names based on their corresponding top BLAST hit in the reference file. Genomic locus coordinates and contig labels are included in the fasta header of each prediction.
2.7. Masking predicted regions
To prevent blastn hits repeatedly mapping to previously annotated regions of the contig, nucleotides in the predicted range are masked with an equivalent number of ‘N’ sites. By default, a range of 750bp upstream and downstream of each prediction is also masked to minimise the occurrence of single loci being annotated as multiple genes.
2.8. Benchmarking pipeline performance
To assess pipeline performance, we conducted tests on a total of fourteen vertebrate species including Homo sapiens (human), Pan troglodytes (chimpanzee), Pongo abelii (orangutan), Mus musculus (mouse), Rattus norvegicus (Brown rat), Canis lupus familiaris (dog), Bos taurus (cow), Monodelphis domestica (opossum), Ornithorhynchus anatinus (platypus), Xenopus laevis (African clawed frog), Anolis carolinensis (green anole), Anas platyrhynchos (mallard duck), and Danio rerio (zebrafish). Mammalian species were chosen based on the representatives available on ‘The Human Olfactory Data Explorer’ (HORDE) database (Olender, Nativ and Lancet, 2013). Coding sequences (CDS) of annotated sensory receptors were obtained from the NCBI RefSeq database for each species (O’Leary et al., 2015).
In each test, Sensommatic predictions were mapped and aligned to their corresponding NCBI annotations with blastn. Three key metrics were used to benchmark performance: recovery, prediction accuracy, and functional classification. Recovery was calculated as the percentage of NCBI receptors recovered by Sensommatic. Functional classification was measured as the proportion of recovered genes sharing functional status with their corresponding NCBI predictions. Prediction accuracy was determined by calculating the percentage of predictions which share 100% identity with their NCBI counterpart. Non-identical hits were further assessed to quantify the percentage of improved predictions. Predictions were considered ‘improved’ when additional transmembrane domains were incorporated, or where start or stop codons were included if lacking in the NCBI prediction. The mean score across these three metrics was used to quantify the final pipeline score. Pipeline specificity was evaluated by verifying that unmapped predictions contain transmembrane domains using deepTMHMM and the NCBI Conserved Domain Database (CDD) search tool (Marchler-Bauer and Bryant, 2004; Marchler-Bauer et al., 2010; Hallgren et al., 2022).
Three sets of tests were conducted to assess the pipeline. The initial test was used to evaluate performance when no adjustments were made to the reference file. In the second test, receptors from each query species were removed from the reference file, to assess how well Sensommatic recovers receptors in novel species. In the third test, NCBI predictions for each species were used as reference files to test pipeline performance in an ideal self-annotation scenario.
2.9. Evaluating pipeline run-times
To evaluate pipeline run-times, Sensommatic was run on the 14 test species and an additional 27 vertebrate species (Supplementary File S1). Genome assemblies were chosen based on varying assembly levels, with representatives from each class of vertebrates. Each test was executed with default parameters, using 20 threads with BLAST on a high-performance computing system equipped with dual Intel Xeon Gold 6152 processors (2.1GHz, 22 cores each) and 384 GB of available memory.