Systematic assessment of homology-based methods for fine-grained functional annotation using diverse protein families

Systematic assessment of homology-based methods for fine-grained functional annotation using diverse protein families | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Systematic assessment of homology-based methods for fine-grained functional annotation using diverse protein families View ORCID Profile Rakesh Busi , View ORCID Profile Pranav Machingal , View ORCID Profile Nandyala Hemachandra , View ORCID Profile Petety V. Balaji doi: https://doi.org/10.1101/2025.09.17.676781 Rakesh Busi 1 Department of Biosciences and Bioengineering, Indian Institute of Technology Bombay , Mumbai, Maharashtra, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Rakesh Busi For correspondence: rakesh.busi{at}iitb.ac.in Pranav Machingal 2 Department of Industrial Engineering and Operations Research, Indian Institute of Technology Bombay , Mumbai, Maharashtra, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Pranav Machingal Nandyala Hemachandra 2 Department of Industrial Engineering and Operations Research, Indian Institute of Technology Bombay , Mumbai, Maharashtra, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Nandyala Hemachandra Petety V. Balaji 1 Department of Biosciences and Bioengineering, Indian Institute of Technology Bombay , Mumbai, Maharashtra, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Petety V. Balaji Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract The size of the protein sequence database is increasing without a consequent increase in the number of proteins with known molecular function, especially at the fine-grained level. Alignment-based approaches such as BLAST and profile hidden Markov models (HMMs) are widely used to infer homology and transfer annotation to be subsequently confirmed by experiments. The ability of BLASTp to distinguish orthologs from paralogs varies across protein families; for profile HMMs, this depends on the sequences considered for generating multiple sequence alignments. In this study, we systematically evaluated the performance of BLASTp and HMM-based methods for fine-grained function annotation using carefully curated protein datasets that are diverse in sequence-structure-function relationships. Expectedly, BLASTp performed well in detecting close homologs but failed to detect remote homologs. BLASTp detected homology between 22.6% and 100% of sequence pairs within different homologous protein families. The extent of sequence identity between trypsin and chymotrypsin sequences is high despite differences in fine-grained molecular function. Transferring function annotation based on homology inferred from BLASTp leads to errors in trypsin-chymotrypsin-like situations. Profile HMMs improved sensitivity and captured subtle homology signals even when sequence identity was low, though some known family members scored below threshold due to functional divergence or mutations at catalytic sites. We further showed that relying solely on homology for annotation transfer can lead to misleading conclusions when proteins have evolved divergent functions despite structural similarity. Our findings highlight that a cautious approach involving BLASTp, profile HMMs, and expert domain knowledge provides the most reliable strategy for functional annotation. After all, not every family member may be doing what we think they are doing. Introduction Gene or genome duplication followed by functional divergence leads to genes encoding new functions. Several models have been proposed to explain functional divergence during evolution. These include (i) neofunc-tionalisation, (ii) subfunctionalisation by duplication, degeneration, and complementation, (iii) subfunction-alisation via specialisation or escape from adaptive conflict, and (iv) innovation-amplification-divergence [ 1 ]. With reference to the function of the ’parent gene’, the new function may be a variant or completely new. Trypsin and chymotrypsin are examples of the former [ 2 ], [ 3 ], and lysozyme c and α -lactalbumin are examples of the latter [ 4 ]. Orthologs arise from speciation events and typically retain the same function in different organisms, showing high sequence identity when the species are closely related. Paralogs, in contrast, originate through duplication within a genome and may diverge in function. Both orthologs and paralogs are homologs, sharing sequence and structural similarities [ 5 ]. However, such similarities often lead to uncertainties in distinguishing paralogs and orthologs from each other. The BLASTp algorithm can be used to search for homologs of a query sequence in a database [ 6 ]. Homology is inferred if (i) a high percentage query sequence is aligned and (ii) sequence identity ≥30% [ 7 ]. Several homologs show high sequence divergence, and sequence identity between such pairs will be <30% 1 . However, such proteins share the same fold. Homologs that show only fold similarity and low sequence identity are considered remote homologs. It is important to note that BLASTp misses remote homologs due to their low sequence identity. Therefore, two sequences that are not aligned by BLASTp are either remote homologs or not homologs. BLASTp is inherently ”blind” to biological functions; it does not assess functional relationships but rather detects evolutionary relationships based on sequence identity. Hence, it is not possible to know whether the hit sequence is an ortholog or a paralog of the query unless additional information is available. Sequence alignments are often not end-to-end because of high sequence divergence (remote homologs), truncations, indels, and/or gene fusions. These factors limit the use of BLASTp for function annotation based on sequence similarity. Profile methods, e.g., PSI-BLAST [ 10 ] and profile HMM [ 11 ], exploit the availability of sequences of several homologs to identify remote homologs with higher confidence relative to what one can achieve by BLASTp. These methods are critically dependent on two factors: (i) ensuring that only homologs are used for building profiles and (ii) the quality of the multiple sequence alignment (MSA). If one were to use the profiles for function annotation, then it becomes important to use only structural and functional homologs for building a profile. Functional homology is essential to avoid inclusion of paralogs. As for the quality of the MSA, currently, there is no widely accepted metric other than visual analysis of sequence conservation and gap positions. BLASTp is being used extensively to detect sequence similarity since its introduction about 3 decades ago. And sequence similarity is used to infer homology and transfer of molecular function annotation. Several studies have shown that homologs can be quite divergent (remote homologs; inference based on fold similarity) and, in many cases, even the function diverges (paralogs). The focus of this study is fine-grained function annotation, i.e., the 4th digit of the EC number. The objective of the present study is to systematically assess if fine-grained molecular function can be inferred based on homology detected by BLASTp and profile HMM. Towards this objective, we have considered 15 datasets, and each dataset consists of two or more protein families. These datasets are divided into 5 groups on the basis of the fold-function relationship. Similarities and differences in fold and function between the protein families show wide variations ( Table 1 ). In all, 5224 single domain sequences representing 26 protein families and protein sets like non-SNARE, non-GPCR (TM), and non-GPCR (no TM) were considered. View this table: View inline View popup Download powerpoint Table 1. Grouping of datasets based on fold and function relationships. In Table 1 , the datasets with the same fold are considered homologous [ 12 ] and datasets with paralogous protein families, i.e., same fold but dissimilar function, were grouped into Group 2 and Group 4, with Group 4 exclusively containing GH-labelled sequences. The datasets in Group 3 contain orthologous families that share the same fold and have similar functions. In contrast, the datasets in Group 1 are not homologous, as they differ in fold and are likely functional analogs arising from convergent evolution. The datasets in Group 5 contain completely unrelated protein sets. Methods Datasets The choice of protein pairs for comparison was based on (i) domain knowledge and (ii) the number of sequences available. This led to the datasets Xylanase, Chitinase, β -Glucosidase, Lysozyme, β -Galactoside, GH2, GH3, and GH5 ( Table 1 ). These constitute Groups 1, 3, and 4. Lysozyme CGCh was specifically considered because of the availability of extensive experimental (ground truth) data. Here, ”extensive” means that we can find at least one literature reference (PubMed ID) to every sequence. Notably, some of the sequences are present in more than one Group ( Table 2 ). View this table: View inline View popup Download powerpoint Table 2. Summary of sequences present in multiple Groups. Protein pairs within any given dataset belonging to Groups 1 to 4 share similarities in at least one of these two characteristics: fold and function. Group 5 datasets were chosen as ’control’ wherein the comparison is between two sets: one set consists of homologs that have the same coarse-grained annotation, and the second set was chosen such that sequences do not share fold or function similarity. This group consists of SNAREs and GPCRs. The SNARE family proteins were chosen based on coarse-grained function annotation [ 13 ]. Proteins were considered without taking into cognizance the category to which they belong, i.e., v-SNARE, t-SNARE, R-SNARE, or Q-SNARE. This suggests that the SNARE family of proteins show a high degree of sequence divergence, and this is in consonance with (i) choosing SNARE proteins irrespective of the category to which they belong and (ii) the large variations in the sequence lengths. The only commonality to all GPCRs sequences considered in this study is the 7-transmembrane (7-TM) domain architecture. Non-GPCRs were chosen in two ways: (i) they have at least one TM domain (non-GPCR™) and (ii) they should not have even one TM domain (non-GPCR(no TM)). Clearly, in Group 5, comparison is between diverse sets of proteins, viz., GPCRs or not, and have/do not have SNARE activity. Fetching and pre-processing of sequences The GenBank IDs for the protein families present in the datasets of Groups 1, 3, and 4 were obtained from the CAZy database, while the corresponding sequences were retrieved from Batch Entrez of NCBI. For the protein families in the Protease dataset, the UniProt IDs were taken from Brenda [ 14 ], and the corresponding sequences were downloaded from UniProt [ 15 ]. The sequences for α -Lactalbumin, SNARE, and non-SNARE protein families were obtained from relevant literature. For all the remaining protein families, both the IDs and sequences were directly fetched from UniProt ( Table 3 ). View this table: View inline View popup Download powerpoint Table 3. Details of datasets, protein families, search keywords, and sources. The fetched sequences were subjected to the following preprocessing steps: (i) sequences with headers containing either the term partial or fragment were removed because we wanted to consider only full length sequences, (ii) sequences with two or more EC numbers were omitted because we did not want to include fusion proteins with ≥2 domains and protein which have more that one activity associated with the same active site, (iii) sequences with duplicate accession numbers between protein families were removed, while within each protein family, only a single copy of entries with duplicate accession numbers or identical sequences was retained, and (iv) sequences containing the letters B, J, O, U, X, and Z were removed, as these letters do not represent any of the 20 standard amino acids. Additionally, outlier sequences, if any, were removed based on the sequence length distribution of the protein family, as determined by a box-whisker plot. Pairwise sequence alignment All-against-all pairwise sequence alignment was performed to determine the percentage sequence identity (PSI) and percentage query coverage (PQC) for all protein sequence pairs of our datasets. NCBI’s Protein-Protein BLAST 2.9.0+ (BLASTp) command-line application was executed on a Linux server using default values for all the parameters unless mentioned to the contrary. The command executed is blastp -query dataset.fasta -subject dataset.fasta -outfmt “6 qseqid sseqid pident qcovhsp” -out result.txt. Here, both the -query and -subject parameters were set to the same file dataset.fasta to enable all-against-all pairwise comparison within a dataset. The output result.txt was generated in tabular format outfmt 6 with four columns: qseqid (query sequence ID), sseqid (subject sequence ID), pident (percentage of identical matches), and qcovhsp (query coverage per high-scoring segment pair (HSP)). An HSP refers to a local, gap-free alignment with one of the highest scores in the search. For bidirectional hits or multiple alignments of the same pair, we retained the result with the highest qcovhsp (maximizing alignment length). Two sequences were inferred to be homologs if the following criteria are met: (i) percentage sequence identity threshold ( T %ID ) ≥30% [ 7 ]; also used in AlphaFold [ 17 ] and Seq2Symm [ 18 ] and (ii) percentage query coverage threshold ( T %QC ) ≥70% ( Fig 1 in Supporting information S1 File). Download figure Open in new tab Fig 1. Flowchart of possible outcomes in fine-grained molecular function annotation. The figure illustrates the decision flow for assigning fine-grained molecular function to a protein using BLASTp against a curated functional database. Three possible outcomes are shown: (i) No homologs found - molecular function remains unassigned, though remote homologs may exist (unsuccessful). (ii) Homologs found with consistent annotations - annotation transfer depends on the granularity of the function (fine-grained or coarse-grained). (iii) Homologs with conflicting annotations - domain knowledge is needed to resolve conflicts, depending on ground truth data availability (resolve). Multiple sequence alignment Algorithms for multiple sequence alignments (MSA) assume that all the input sequences are homologs, i.e., they have evolved from a common ancestor. A column in an MSA therefore shows the residues in the input sequences that map to the same amino acid position in the ancestral sequence [ 19 ]. However, it becomes increasingly difficult to identify (align) homologous positions in input sequences as the level of sequence divergence increases. This is particularly so if the number of indel events is high. Unfortunately, there is no mathematical or statistical parameter (something similar to the e-value that BLAST uses) that one can use to assess the quality of MSAs. A set of benchmarks, such as HOMSTRAD and/or BAliBASE databases, is typically used to compare the performance of various MSA algorithms. The accuracy and cost of nine popular MSA programs, namely CLUSTALW, CLUSTAL OMEGA, DIALIGN-TX, MAFFT, MUSCLE, POA, PROBALIGN, PROBCONS, and T-COFFEE, were compared against the benchmark alignment dataset BAliBASE [ 20 ]. The accuracy of alignment was calculated as the sum-of-pairs and total-column scores. The computational costs were determined by collecting peak memory usage and time taken for the execution of the program. This comparative study found that consistency-based programs, viz., PROBCONS, T-COFFEE, PROBALIGN, and MAFFT, outperformed the rest. It further found that MAFFT can deliver faster and reliable alignments. In another study [ 21 ], the authors highlighted that MAFFT can consistently produce correct alignments with high accuracy in comparison to the other five MSA programs, viz., CLUSTALW, DIALIGN-T, MUSCLE, PROBCONS, and T-COFFEE. Hence, MAFFT v7.490(2012/Oct/30), installed on a Linux server, was used to obtain MSAs. The mafft command was used to generate an alignment file (.aln) for the given protein family. This alignment file was visually analyzed using Jalview software [ 22 ] and further examined for gaps and conservation. Default values were used for all the parameters. Hidden Markov models We utilized profile hidden Markov models (HMMs) [ 23 ], [ 11 ], [ 24 ] to determine if a given protein sequence belongs to the protein family represented by the profile HMM. HMMER 3.3 software package was installed on a Linux server. The workflow has four steps: A multiple sequence alignment (MSA) of the set of homologous sequences belonging to a protein family was generated using the MAFFT alignment tool [ 25 ]. A separate profile HMM was built for each protein family from the MSA using the hmmbuild command. The bit score threshold for each profile HMM was set by first performing self-scoring, i.e., scoring every sequence that was used to build the profile HMM against the same profile HMM. This gave a range of best 1 domain bit scores ( Tables 1 and 2 in Supporting information S1 File), and these ranges were used as guides to plot ROC curves ( Figs 2 -6 in Supporting information S1 File). Bit score thresholds ( T ROC ) were set such that there are no false positives when the search space is the collection of sequences from all the datasets used in the present study (Supporting information S1 Table). We note that T ROC is dependent on the search space. Every query sequence was scored against the profile HMM using the hmmsearch command. If the bit score of a query sequence was ≥ T ROC of the protein family’s profile HMM, the sequence was inferred to belong to that family. In other words, the query sequence was considered homologous to the sequences in the protein family. We used the best 1 domain bit score for further analysis for the following reasons: (i) Unlike E-value, bit score does not depend on the size of the search space. Hence, we get the same value whether we score a sequence by itself or as a part of a collection of sequences. (ii) The best 1 domain score contains the bit score for the single best scoring domain, unlike the score which sums up scores for all (partial) alignments as well as sequences repeats, if any, present in a sequence ( https://eddylab.org/software/hmmer/Userguide.pdf ). Download figure Open in new tab Fig 2. Sequence-structure-function relationships in protein families across datasets. The protein families included in the datasets exhibit different combinations of sequence, structure, and function similarity, which can be either similar or dissimilar. These datasets are organized into five groups, covering all combinations except three. The missing combinations are: (i) sequence similarity without structural and functional similarity, typically observed in metamorphic proteins (inadequate data) [ 32 ], [ 33 ]; (ii) sequence and functional similarity without structural similarity, a relationship which we are not aware (inadequate data); and (iii) similarity in sequence, structure, and function, which defines proteins belonging to the same family in this study. Results We present the results as follows: (A) Detection of homology of every sequence with other sequences (i) within the family, (ii) of protein families within the dataset, and (iii) from other families considered in this study. (B) To assess whether annotations transferred based on inference are correct. We have chosen proteins with known functions so that we can find out if transferred annotations are correct or not. Based on such an assessment, we discuss the suitability of BLASTp and profile HMM when the query sequences are of unknown function (fine-grained). Using BLASTp to assess if sequences considered in a protein family are homologs Datasets were created in such a way that proteins within each protein family are homologous to each other, as inferred from their shared fold. However, the percentage of alignment pairs exceeding the BLASTp threshold for homology within a family (% H intra ) varies from 22.6% (Cellulase; Group 4) to 100.0% (Lysozyme G; Group 1) (see % H intra in Tables 4 and 5 ). This variation is largely due to the sequence divergence caused by the broad taxonomic distribution leading to remote homology. For example, the Cellulase family in the GH5 dataset (Group 4) includes sequences from Eukaryota, Bacteria, and Archaea; some sequences are Unclassified. GH10 family in Xylanase dataset (Group 1) displays moderate conservation (% H intra = 52.1%), includes Eukaryota and Bacteria. In contrast, the sequences from the Lysozyme G family of Lysozyme CGCh (Group 1) are all from Eukaryota. These observations are consistent with the well-documented limitation of BLASTp to detect remote homology. View this table: View inline View popup Download powerpoint Table 4. Structural, functional, and sequence information of protein families in Groups 1–3. View this table: View inline View popup Download powerpoint Table 5. Structural, functional, and sequence information of protein families in Groups 4 and 5. Using BLASTp to assess if protein families of a dataset are homologs The protein families in the datasets of Groups 2, 3, and 4 are homologous to each other due to the same fold ( Table 1 ). BLASTp is able to detect homology in Group 2 (% H inter ranges from 56.9% to 89.3%, Tables 4 and 5 ) but not in Group 3 and Group 4 (% H inter either 0.0% or close to 0.0%). This is because of the inability of BLASTp to detect remote homologs caused by sequence divergence. Sequence divergence can be due to broader taxonomic distribution of sequences as well as the differences in rate of evolution of different protein families (section 3.3c in [ 26 ]). For example, in Group 2, the protein sequences in all datasets were from Eukaryota, whereas those in Group 3 and Group 4 were from Archaea, Bacteria, Eukaryota, and Viruses (only in Group 3; Lysozyme dataset), along with some unclassified. Protein families in Group 1 datasets are non-homologous by design ( Table 1 ). BLASTp not detecting any homologs is as expected. Even in Group 5, the choice of proteins is such that there is no homology between SNARE and non-SNARE (as described in the Methods section); so is the case with GPCR and non-GPCR (TM)/non-GPCR (no TM). Hence, % H inter = 0.0% as expected (see % H inter in Table 5 ). Alignment of every sequence with every other sequence considered in this study using BLASTp So far, we have analyzed results in the specific context of sequence identity within and between the protein families of a dataset. We also analysed results when sequences were aligned with sequences belonging to other families considered in this study. BLASTp finds homology of several sequences from Groups 1-4 with sequences belonging to the non-SNARE and non-GPCR sequences of Group 5 (Supporting information S2 Table). These are not unexpected because non-SNARE and non-GPCR sequences were chosen with the only condition that they are not SNAREs/GPCRs (with or without TM domain) as described in the method section. Assessment of annotations based on homology inferred from BLASTp alignments We present the assessment using the workflow shown in Fig 1 . The assessment considers a search space which is the collection of sequences from all the datasets used in the present study. This assessment is applicable as is even when faced with a protein of unknown function, its sequence is typically used as a query in a BLASTp search against a reference database. This database can be a curated collection (e.g., SwissProt, CAZy, UniRef, or a user-defined set) or a large aggregated resource containing all known sequences (e.g., NCBI nr). Consider a protein for which we are seeking to assign fine-grained molecular function. The sequence of this protein is used as the query for BLASTp search against a database consisting of proteins with curated molecular function. We can envisage the possible outcomes ( Fig 1 ): (i) No homologs are found. One may choose other search spaces to find homologs. Assuming that the search spaces represent all known molecular functions, BLASTp is of no help in assigning molecular function. We may note that the presence of remote homologs cannot be ruled out (Outcome: unsuccessful). (ii) Homologs are found with consistent molecular function. In this case, transfer of annotation will be justified, and the granularity of function annotation is determined by annotations of homologs (Outcome: fine-grained or Outcome: coarse-grained). (iii) Annotations of the homologs are not identical. Here, one needs to use domain knowledge to resolve the conflicting annotation (Outcome: resolve) but its success depends on the nature of the available ground truth data. Herein, we analyse the outcome of the above workflow for each sequence when queried against the search space comprising all other sequences considered in this study (Supporting information S2 Table). We also provide the detailed annotations of subject sequences that are different from the query sequence. These details are taken from the UniProt. Analysis of query sequences whose hits do not have the same annotations (Outcome: resolve) is given below, with additional details for this analysis given in Supporting information S3 Table. Lysozyme C, α -lactalbumin, and sperm acrosome membrane-associated protein 3 (Spaca3): these are sequence and structural homologs that perform completely different functions. Resolution by a domain expert depends upon the level of sequence identity of the query with these three functionally different groups of proteins and whether active site residues are conserved. Lysozyme C and GH22 family proteins: All these proteins catalyse the same reaction (same at the fine-grained level) but are involved in different biological processes. At the coarse-grained level, they are the same with respect to the reaction catalysed (glycoside hydrolase) but different with respect to the biological process. GH19 family chitinase and endochitinase 1: Both are glycoside hydrolases (same at coarse-grained level) but differ in their reaction specificity. Trypsin, chymotrypsin, and serine protease/peptidase: All these proteins catalyse the same reaction, hydrolysis of the peptide bond (same in coarse-grained level), but differ from each other in terms of their cleavage site and biological processes. FDX1 and FDX2: These two isoforms are highly specific for their substrates and participate in distinct biochemical pathways. The distinct conserved sequence motifs should be used to annotate FDX1 or FDX2 while resolving the conflicting annotation [ 30 ], [ 28 ]. GH1 family β -glucosidase and myrosinase 2: The former is a β -glucosidase, whereas the latter is a thioglucosidase. Whether each of the two enzymes can also act on the substrates of the other enzyme is not known. Hence, the available data is inadequate to resolve. GH3 family β -glucosidases and 1,4- β -xylosidases: β -glucosidases are exo-enzymes and some are known to have broader substrate specificity, i.e., hydrolyse β -D-galactosides, α -L-arabinosides, β -D-xylosides, and β -D-fucosides also (same at coarse-grained level). Hence, an unambiguous resolution of conflicting annotations is not possible. Cellulase and endo- β -mannanases: Both are glycoside hydrolases (same at coarse-grained level) but differ in substrate specificity. Hence, an unambiguous resolution of conflicting annotations is not possible. Using profile HMMs to assess if sequences considered in a protein family are homologs Sequences belonging to each protein family were scored against the corresponding profile HMM ( Table 6 ). The number of sequences whose bit score against the profile HMM is lower than T ROC even though these sequences were used to build the profile HMM, is denoted as N low . All the sequences of the family meet the bit score threshold in nearly half of the cases, i.e., N low = 0. For other families N low ≠ 0 because our aim is to choose T ROC such as to make FPR = 0 i.e., avoid propagation of wrong annotation even if it means the bit score for some of the sequences of the family < T ROC . View this table: View inline View popup Download powerpoint Table 6. Profile HMM thresholds and performance on protein families of all Groups. A possible reason as to why N low ≠ 0 for β -Glucosidase-GH3 HMM and 1,4- β -Xylosidase-GH3 HMM is the similarity shared by their substrates at the monomeric level: the substituent at C-5 is -CH 2 OH in glucopyranose and -H in xylopyranose. Several β -xylosidases score higher than β -glucosidases against the β -Glucosidase-GH3 HMM and vice versa for the 1,4- β -Xylosidase-GH3 HMM. In the case of β -Glucosidase-GH1 HMM, N low = 56 because myrosinase 2 (UniProt ID: Q9C5C2), a thioglucosidase (EC 3.2.1.147) from Non-GPCR (no TM) set, scores 527.2 bits score against the β -Glucosidase-GH1 HMM. The N low ≠ 0 for Lysozyme-GH22 HMM because of two possible reasons: (i) sperm acrosome membrane-associated protein 3 (UniProt: Q9D9X8) scores ≥ T ROC . It is a sequence and structural homolog of lysozyme C which has lost enzymatic activity because of the mutation of catalytic site residues [ 31 ]. (ii) 13 of the 135 sequences of this family score much below the T ROC because of high sequence divergence relative to the other 122 sequences. Using profile HMMs to assess homology across protein families and datasets While using BLASTp, thresholds for T %ID and T %QC to infer homology between the query and subject sequences are independent of the sequences being aligned. In contrast, each profile HMM will have a customized bit score threshold. As mentioned in the Methods section, bit score thresholds were set such that there are no false positives when the search space is the collection of sequences from all the datasets used in the present study. In view of this, none of the sequences score higher than T ROC for profile HMMs of other families. For example, if we consider the Protease dataset, none of the chymotrypsins are hits for Trypsin HMM and vice versa. Because of these reasons, setting a bit score threshold carefully can eliminate the possibility of erroneous transfer of function annotation. The workflow for searching a database using a profile HMM for the purpose of transferring annotation (Fig 7 in Supporting information S1 File) is broadly similar to the one described for BLASTp ( Fig 1 ) except that the we use bit score threshold ( T ROC ) instead of percentage sequence identity threshold ( T %ID ) and percentage query coverage threshold ( T %QC ). Discussion and Conclusion High-throughput genome sequencing has led to the accumulation of a very large number of conceptually translated protein sequences. Knowledge of the function(s) performed by such proteins is essential to fully realise the value of high-throughput sequencing projects. It is neither practical to experimentally determine the function(s) of such a large number of proteins nor desirable. Nevertheless, the sequence data can be effectively used to assign molecular functions to proteins of interest, and this can be experimentally validated. BLASTp and profile-based methods (PSI-BLAST, profile HMM) are widely used to assign function based on homology. However, assignment of molecular function based on homology may be erroneous because these methods do not by default distinguish orthologs from paralogs. Hence, in this study, we have investigated the suitability or effectiveness of BLASTp and profile HMM to assign fine-grained molecular function based on inferred homology using carefully curated datasets that cover diverse sequence-structurefunction relationships ( Table 1 ) ( Fig 2 ). BLASTp thresholds ( T %ID and T %QC ) to infer homology between query and subject sequences are inde-pendent of the sequences being aligned. Because of this reason we do find homologs which do not share function, for example, in Group 2 datasets. Transfer of function annotation in these cases will obviously be erroneous. In Group 4 datasets, homology cannot be inferred because of high sequence divergence, and hence the question of transferring annotation does not arise. In contrast to BLASTp, profile HMMs have customized bit score thresholds for inferring homologs. A sequence that scores above this threshold is inferred as belonging to the protein family represented by the profile HMM. Thresholds may be stringent so as to make the profile ’specific’ (eliminate false positives) or relaxed so that the profile is sensitive (eliminate false negatives). In this study, we chose to prioritize specificity over sensitivity. It is possible to adopt a two-step workflow wherein a relaxed threshold is used in the first step so that there are no false negatives. In the second step, hits from the first step that do not meet the stringent threshold may be taken up for further analysis. One may choose to generate separate profiles for each paralogous family to eliminate the possibility of erroneous transfer of function as illustrated by the Protease dataset: none of the chymotrypsin sequences are hits for the Trypsin HMM and vice versa. We have presented a general workflow to annotate protein sequences using BLASTp or profile HMM against any search space. In our study, a collection of all sequences except the one being used as the query constituted the search space. The question of transferring annotation does not arise if the query does not have any homologs in the search space. If homologs are found, then the decision to transfer annotation is straightforward if all hits have the same annotation. This is not the case when hits with different annotations are obtained; differences which may be due to (i) mutation at the active site leading to loss or change of function (e.g., lysozyme C, α -lactalbumin, and Spaca3), (ii) inadequate experimental data leading to ambiguous annotation (e.g., GH1 family β -glucosidase and myrosinase), and (iii) function being same at coarse-grained level but different at fine-grained level (e.g., trypsin and chymotrypsin) due to differences in substrate utilized, reaction catalysed, and/or biological process. Such conflicting annotations have to be resolved by domain experts. However, this may or may not be possible as illustrated below. The case of Spaca3 meeting the threshold of Lysozyme-GH22 HMM is easier to resolve because of the mutation of active site residues in Spaca3. Checking for the conservation of active site residues will suffice to infer whether hits of this Lysozyme-GH22 HMM are Spaca3. β -glucosidase and 1,4- β -xylosidase satisfy the thresholds of both β -Glucosidase-GH3 HMM and 1,4- β -Xylosidase-GH3 HMM. Activities of all β -glucosidases against xylosides have not been tested; similarly, activities of all β -xylosidases have not been tested against glucosides. It may turn out that β -glucosidase and 1,4- β -xylosidase are orthologs with broad substrate specificity but these are considered as paralogs. The current knowledge of the sequencesubstrate specificity relationship is inadequate to resolve whether these enzymes are paralogs or orthologs without experimental data. Another such example is β -glucosidase of β -Glucosidase-GH1 HMM which acts on O-glucosides, and myrosinase 2 which acts on thioglucosides. In general, enzymes are known to display promiscuity with respect to substrate utilized, reaction catalysed, and/or product formed [ 34 ], [ 35 ]. Often, data on promiscuity is not available because of the additional resources required to test for promiscuity. Paucity of such data limits the use of profile HMMs for transferring annotation. Several proteins are known to have moonlighting function ( http://www.moonlightingproteins.org/ ) [ 36 ]. Moonlighting, unlike promiscuity, does not affect bit score thresholds. One can envisage two possibilities: within a family of homologous proteins, some but not all have moonlighting function. This leads to uncertainties about transferring annotations related to additional (moonlighting) functions. We also note that just because a protein can perform a function does not necessarily mean it does so in vivo, as shown for SipB [ 37 ]. This illustrates that although transfer of function annotation based on inferred homology seems to be justified, the protein may not carry out the annotated function in vivo. Addressing challenges such as these requires experimental data. In conclusion, sequence divergence, variations in protein lengths, and annotation granularity within a protein family mean that considering a larger number of proteins for generating profile HMMs makes the profile less reliable if the interest is in transferring function at the fine-grained level. Erroneous function transfer may occur, especially at the fine-grained level due to limited experimental data, subtle sequence-function relationships leading to substrate, reaction, and product promiscuities, and biological ingenuities such as moonlighting proteins. Nevertheless, the present study demonstrates the use of BLASTp and profile HMMs to annotate protein functions while highlighting a set of challenging cases that preclude automation. Conflicts of Interests None declared. Funding The authors received no specific grant from any funding agency. Supporting information S1 File S1 Table S2 Table S3 Table Acknowledgment We thank the Indian Institute of Technology Bombay for the infrastructure and other facilities. Footnotes ↵ [ 1 some researchers set the sequence identity threshold to 25% [ 8 ], [ 9 ] References [1]. ↵ Glasner ME , Truong DP , Morse BC . How enzyme promiscuity and horizontal gene transfer contribute to metabolic innovation . The FEBS journal . 2020 ; 287 ( 7 ): 1323 – 42 . OpenUrl CrossRef PubMed [2]. ↵ Hedstrom L , Szilagyi L , Rutter WJ . Converting trypsin to chymotrypsin: the role of surface loops . Science . 1992 ; 255 ( 5049 ): 1249 – 53 . OpenUrl Abstract / FREE Full Text [3]. ↵ Ma W , Tang C , Lai L. Specificity of trypsin and chymotrypsin: loop-motion-controlled dynamic correlation as a determinant . Biophysical journal . 2005 ; 89 ( 2 ): 1183 – 93 . OpenUrl CrossRef PubMed Web of Science [4]. ↵ Qasba PK , Kumar S , Brew K. Molecular divergence of lysozymes and α-lactalbumin . Critical reviews in biochemistry and molecular biology . 1997 ; 32 ( 4 ): 255 – 306 . OpenUrl CrossRef PubMed Web of Science [5]. ↵ Koonin EV . Orthologs, paralogs, and evolutionary genomics . Annu Rev Genet . 2005 ; 39 ( 1 ): 309 – 38 . OpenUrl CrossRef PubMed Web of Science [6]. ↵ Altschul SF , Gish W , Miller W , Myers EW , Lipman DJ . Basic local alignment search tool . Journal of Molecular Biology . 1990 ; 215 ( 3 ): 403 – 10 . Available from: https://www.sciencedirect.com/science/article/pii/S0022283605803602 . OpenUrl CrossRef PubMed Web of Science [7]. ↵ Rost B. Twilight zone of protein sequence alignments . Protein engineering . 1999 ; 12 ( 2 ): 85 – 94 . OpenUrl CrossRef PubMed Web of Science [8]. ↵ Sander C , Schneider R. Database of homology-derived protein structures and the structural meaning of sequence alignment . Proteins: Structure, Function, and Bioinformatics . 1991 ; 9 ( 1 ): 56 – 68 . OpenUrl [9]. ↵ Chung SY , Subbiah S. A structural explanation for the twilight zone of protein sequence homology . Structure . 1996 ; 4 ( 10 ): 1123 – 7 . OpenUrl CrossRef PubMed [10]. ↵ Altschul SF , Madden TL , Schäffer AA , Zhang J , Zhang Z , Miller W , et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs . Nucleic acids research . 1997 ; 25 ( 17 ): 3389 – 402 . OpenUrl CrossRef PubMed Web of Science [11]. ↵ Eddy SR . Profile hidden Markov models . Bioinformatics (Oxford, England) . 1998 ; 14 ( 9 ): 755 – 63 . OpenUrl CrossRef PubMed Web of Science [12]. ↵ Pearson WR , Sierk ML . The limits of protein sequence comparison? Current opinion in structural biology . 2005 ; 15 ( 3 ): 254 – 60 . OpenUrl CrossRef PubMed Web of Science [13]. ↵ Auriemma Citarella A , Di Biasi L , Risi M , Tortora G. SNARER: new molecular descriptors for SNARE proteins classification . BMC bioinformatics . 2022 ; 23 ( 1 ): 1 – 20 . OpenUrl CrossRef PubMed [14]. ↵ Chang A , Jeske L , Ulbrich S , Hofmann J , Koblitz J , Schomburg I , et al. BRENDA, the ELIXIR core data resource in 2021: new developments and updates . Nucleic acids research . 2021 ; 49 ( D1 ): D498 – 508 . OpenUrl CrossRef PubMed [15]. ↵ Consortium U. UniProt: a worldwide hub of protein knowledge . Nucleic acids research . 2019 ; 47 ( D1 ): D506 – 15 . OpenUrl CrossRef PubMed [16]. Srivastava J , Hembrom R , Kumawat A , Balaji PV . Searching sequence databases for functional homologs using profile HMMs: how to set bit score thresholds? bioRxiv . 2021 . [17]. ↵ Jumper J , Evans R , Pritzel A , Green T , Figurnov M , Ronneberger O , et al. Highly accurate protein structure prediction with AlphaFold . nature . 2021 ; 596 ( 7873 ): 583 – 9 . OpenUrl CrossRef PubMed [18]. ↵ Kshirsagar M , Meller A , Humphreys IR , Sledzieski S , Xu Y , Dodhia R , et al. Rapid and accurate prediction of protein homo-oligomer symmetry using Seq2Symm . Nature Communications . 2025 ; 16 ( 1 ): 2017 . OpenUrl PubMed [19]. ↵ Mount DW . Bioinformatics-sequence and genome analysis . 2nd ed. Cold Spring Harbor Laboratory Press ; 2004 . [20]. ↵ Pais FSM , Ruy PdC , Oliveira G , Coimbra RS . Assessing the efficiency of multiple sequence alignment programs . Algorithms for molecular biology . 2014 ; 9 : 1 – 8 . OpenUrl PubMed [21]. ↵ Golubchik T , Wise MJ , Easteal S , Jermiin LS . Mind the gaps: evidence of bias in estimates of multiple sequence alignments . Molecular biology and evolution . 2007 ; 24 ( 11 ): 2433 – 42 . OpenUrl CrossRef PubMed Web of Science [22]. ↵ Waterhouse AM , Procter JB , Martin DM , Clamp M , Barton GJ . Jalview Version 2—a multiple sequence alignment editor and analysis workbench . Bioinformatics . 2009 ; 25 ( 9 ): 1189 – 91 . OpenUrl CrossRef PubMed Web of Science [23]. ↵ Eddy SR . Hidden markov models . Current opinion in structural biology . 1996 ; 6 ( 3 ): 361 – 5 . OpenUrl CrossRef PubMed Web of Science [24]. ↵ Eddy SR . What is a hidden Markov model? Nature biotechnology . 2004 ; 22 ( 10 ): 1315 – 6 . OpenUrl CrossRef PubMed Web of Science [25]. ↵ Katoh K , Rozewicki J , Yamada KD . MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization . Briefings in bioinformatics . 2019 ; 20 ( 4 ): 1160 – 6 . OpenUrl CrossRef PubMed [26]. ↵ Creighton TE . Proteins: structures and molecular properties . Macmillan ; 1993 . [27]. Goettig P , Brandstetter H , Magdolen V. Surface loops of trypsin-like serine proteases as determinants of function . Biochimie . 2019 ; 166 : 52 – 76 . OpenUrl CrossRef PubMed [28]. ↵ Schulz V , Basu S , Freibert SA , Webert H , Boss L , Mühlenhoff U , et al. Functional spectrum and specificity of mitochondrial ferredoxins FDX1 and FDX2 . Nature Chemical Biology . 2023 ; 19 ( 2 ): 206 – 17 . OpenUrl PubMed [29]. Zhang M , Chen T , Lu X , Lan X , Chen Z , Lu S. G protein-coupled receptors (GPCRs): advances in structures, mechanisms, and drug discovery . Signal Transduction and Targeted Therapy . 2024 ; 9 ( 1 ): 88 . OpenUrl [30]. ↵ Sheftel AD , Stehling O , Pierik AJ , Elsässer HP , Mühlenhoff U , Webert H , et al. Humans possess two mitochondrial ferredoxins, Fdx1 and Fdx2, with distinct roles in steroidogenesis, heme, and Fe/S cluster biosynthesis . Proceedings of the National Academy of Sciences . 2010 ; 107 ( 26 ): 11775 – 80 . OpenUrl Abstract / FREE Full Text [31]. ↵ Mandal A , Klotz KL , Shetty J , Jayes FL , Wolkowicz MJ , Bolling LC , et al. SLLP1, a unique, intraacrosomal, non-bacteriolytic, c lysozyme-like protein of human spermatozoa . Biology of reproduction . 2003 ; 68 ( 5 ): 1525 – 37 . OpenUrl PubMed Web of Science [32]. ↵ Dishman AF , Tyler RC , Fox JC , Kleist AB , Prehoda KE , Babu MM , et al. Evolution of fold switching in a metamorphic protein . Science . 2021 ; 371 ( 6524 ): 86 – 90 . OpenUrl Abstract / FREE Full Text [33]. ↵ Dishman AF , Volkman BF . Metamorphic protein folding as evolutionary adaptation . Trends in biochemical sciences . 2023 ; 48 ( 8 ): 665 – 72 . OpenUrl CrossRef PubMed [34]. ↵ O’Brien PJ , Herschlag D. Catalytic promiscuity and the evolution of new enzymatic activities . Chemistry & biology . 1999 ; 6 ( 4 ): R91 – R105 . OpenUrl CrossRef PubMed Web of Science [35]. ↵ Tawfik OK , S D. Enzyme promiscuity: a mechanistic and evolutionary perspective . Annual review of biochemistry . 2010 ; 79 ( 1 ): 471 – 505 . OpenUrl CrossRef PubMed Web of Science [36]. ↵ Jeffery CJ . Moonlighting proteins . Trends in biochemical sciences . 1999 ; 24 ( 1 ): 8 – 11 . OpenUrl CrossRef PubMed Web of Science [37]. ↵ Strauss EJ , Falkow S. Microbial pathogenesis: genomics and beyond . Science . 1997 ; 276 ( 5313 ): 707 – 12 . OpenUrl Abstract / FREE Full Text View the discussion thread. Back to top Previous Next Posted September 20, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Systematic assessment of homology-based methods for fine-grained functional annotation using diverse protein families Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Systematic assessment of homology-based methods for fine-grained functional annotation using diverse protein families Rakesh Busi , Pranav Machingal , Nandyala Hemachandra , Petety V. Balaji bioRxiv 2025.09.17.676781; doi: https://doi.org/10.1101/2025.09.17.676781 Share This Article: Copy Citation Tools Systematic assessment of homology-based methods for fine-grained functional annotation using diverse protein families Rakesh Busi , Pranav Machingal , Nandyala Hemachandra , Petety V. Balaji bioRxiv 2025.09.17.676781; doi: https://doi.org/10.1101/2025.09.17.676781 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41911) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13371) Ecology (19887) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40364) Molecular Biology (17163) Neuroscience (88537) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15129) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9817) Zoology (2269)