STAR★Methods

Key resources table

REAGENT or RESOURCE | SOURCE | IDENTIFIER
Deposited data

Experimentally validated virulent proteins | PHI-Base | http://www.phi-base.org/
Victors | http://www.phidias.us/victors/
DFVF | http://sysbio.unl.edu/DFVF/
Non-virulent proteins | UniProt | https://www.uniprot.org/
Code availability | Present study (GitHub) | https://github.com/ekjotkaurm/FunVFPred
FunVFPred GUI | Present study | http://14.139.59.213:8501/funvfpred/
FunVFPred Docker | Present study | https://github.com/ekjotkaurm/FunVFPred-GUI

Software and algorithms

SPAAN | Sachdeva et al., 2005 | https://academic.oup.com/bioinformatics/article/21/4/483/202900
python3 | Python | https://www.python.org/
numPy | NumPy | https://numpy.org/
pandas | Pandas | https://pandas.pydata.org/
sklearn | sklearn | https://scikit-learn.org/stable/
Random Forest | RandomForest | https://scikitlearn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier/
Tensorflow | Tensorflow | https://www.tensorflow.org/
metaplotlib | Metaplotlib | https://matplotlib.org/
Torch | pytorch | https://pytorch.org/
biopython | biopython | https://biopython.org/
Joblib | joblib | https://joblib.readthedocs.io/en/latest/installing.html
tape-proteins | tape-proteins | https://pypi.org/project/tape-proteins/
Scipy | scipy | https://scipy.org/

Experimental model and study participant details

No experimental model system used for this study.

Method details

Data collection

For this investigation, three notable databases on fungal pathogens were used: the Database of Fungal Virulence Factors (DFVF),15 the Victors database,35 and the Pathogen–Host Interaction Database (PHI-Base).36 To ensure the quality and relevance of the collected data, stringent screening procedures were applied to the protein sequences extracted from these databases. Protein sequences experimentally validated as virulence factors (VFs) were collected from DFVF, Victors, and PHI-Base. These validations were based on functional evidence reported in the literature, including gene knockout experiments, infection models, or biochemical assays confirming pathogenic roles.

Initially, 813 experimentally validated pathogenic protein sequences or fungal virulence factors (VFs) were extracted from eight Candida spp., including Candida albicans (strain SC5314/ATCC MYA-2876), C. glycerinogenes, Pichia kudriavzevii (previously known as Candida krusei) (Taxon ID: 4909), C. dubliniensis (strain CD36/ATCC MYA-646/CBS 7987/NCPF 3949/NRRL Y-17841), C. tropicalis (Taxon ID: 5482), C. parapsilosis (Taxon ID: 5480), C. orthopsilosis (Taxon ID: 273371), C. glabrata (Taxon ID: 5478), and C. metapsilosis (Taxon ID: 273372). Redundancy removal was performed using CD-HIT,37 with a 100% sequence identity cutoff, resulting in 616 unique virulent proteins (Table 6).

The whole proteomes of these eight Candida spp. were downloaded from UniProt.38 Proteins overlapping with the positive dataset were excluded. The remaining proteins, which had not been experimentally characterized as virulent, were treated as putatively non-virulent. The final negative dataset comprised 57,841 proteins after applying a 100% sequence identity threshold and removing redundancy using the same CD-HIT approach.

This careful curation ensured that the positive dataset contained only experimentally validated virulent proteins, while the negative dataset was free from sequence overlap with the positive set. Minimizing redundancy within the two datasets resulted in high-quality input data for downstream feature extraction and predictive model development.

Data balancing

In real-world scenarios, there is a substantial imbalance between the number of virulent and non-virulent proteins, with non-virulent samples significantly outnumbering virulent proteins. Initially, our dataset contained 813 virulent and 57,841 non-virulent proteins. To investigate the impact of this imbalance on prediction performance, we created a reduced subset of the training data with a 1:10 ratio of virulent to non-virulent proteins. The same ratio was applied to the test and validation datasets.

This imbalanced subset was used to train the model and evaluate its performance. The results showed a strong bias toward predicting the majority (non-virulent) class, highlighting the importance of dataset balancing for accurate and fair classification.

To address this issue, we employed a random undersampling strategy to balance the classes at a 1:1 ratio. Specifically, non-virulent sequences were randomly sampled using a fixed seed to ensure reproducibility, and the final balanced dataset was prepared and saved using Biopython's SeqIO module. This approach ensured equal representation of both classes during model training, reducing prediction bias and improving classification balance.

Feature encoding

To enhance the prediction of VFs, we used sequence-based features, including Amino Acid Composition (AAC), Dipeptide Deviation from Expected Mean (DDE), and UniRep embeddings. UniRep is a pretrained deep learning model trained on a large corpus of protein sequences and is capable of capturing structural, evolutionary, and biophysical information from protein sequences.

The rationale for combining these features is both biological and technical: AAC captures overall amino acid usage and associated physicochemical properties, DDE captures local sequence dependencies that may be relevant to functional motifs, and UniRep embeddings provide high-dimensional learned representations of protein sequence properties. Collectively, these complementary features allow models to better capture the biological diversity of virulence-associated proteins.

Amino acid composition (AAC)

The frequencies of the 20 standard natural amino acids (i.e., “ACDEFGHIKLMNPQRSTVWY”) within a protein or peptide sequence are expressed as the amino acid composition (AAC) feature.39 The formula below was used to calculate AAC:where N(a) indicates the number of a particular amino acid that exist are, N is the protein or peptide's sequence length, and f(a) is the final 20-dimensional feature vector that was generated.

f(a)=N(a)N,a∈{A,M,L,…..,P}

Protein sequences were read from a FASTA file using Biopython (version 1.81), and AAC values were calculated using NumPy for efficient computation. Each feature vector was then linked to its corresponding class label (virulent = 1, non-virulent = 0), as provided in a CSV file containing protein identifiers. The final dataset, combining AAC features and class labels, was saved in CSV format for use in ML models.

Dipeptide deviation from expected mean (DDE)

The DDE feature comprises three components: the theoretical mean (TM), dipeptide composition (DPC), and theoretical variance (TV).40 The TM characteristic was calculated as follows:where CN equals 61, which indicates the total number of potential codons excluding the three stop codons, and Ca and Cb, respectively, indicate the codon numbers encoding amino acids a and b. The DPC feature calculation makes use of the above explanation. This is how the TV feature is computed:

TM(a,b)=CaCN×CbCN

TV(a,b)=TM(a,b)(1−TM9a,b))N−1

Subsequently, the following formula is employed to calculate DDE(a, b):

DDE(a,b)=DPC(a,b)−TM(a,b)TV(a,b)

The DDE features were computed using an in-house Python script employing the Biopython package. Each protein sequence was parsed from a FASTA file, and its dipeptide composition was calculated. The generated features were merged with corresponding labels and saved in CSV format for downstream classification tasks.

Deep representation learning embeddings

UniRep32 is a deep representation learning method that uses a multidimensional Long Short-Term Memory (mLSTM) model with 1900 dimensions. It was trained on the UniRef50 protein database to learn complex features of amino acid sequences. High-dimensional embeddings that represent structural and functional characteristics of proteins are generated using this approach. UniRep provides a flexible representation of protein sequences and has demonstrated potential for applications in protein engineering and computational biology.

The UniRep feature calculation process can be summarized as follows:1.Sequence encoding: Each amino acid in the protein sequence was represented using a 20-dimensional one-hot encoding vector.2.Feature extraction using mLSTM: For each amino acid position t, a multiplicative long short-term memory (mLSTM) network processes the sequence to produce hidden states (ht).where, x__t is input vector (one-hot encoding of amino acid at position t), h__t-1 is previous hidden state and c__t-1 is previous cell state3.Sequence representation: For a protein sequence of length T, the final UniRep embedding can be either:4.The last state (h__T), or5.The average of all hidden states (mean pooling):

ht=mLSTM(xt,ht−1,ct−1)

UniRep=1T∑t=1Tht

The UniRep features were extracted using the TAPE (Tasks Assessing Protein Embeddings) library in Python. Protein sequences were parsed from a FASTA file and encoded using a UniRep-compatible tokenizer. The pre-trained UniRep model (babbler-1900) was utilized to generate 1900-dimensional embeddings by averaging the hidden states across each sequence. The resulting feature vectors were combined with protein identifiers and corresponding labels from the input dataset. The final feature matrix was saved in CSV format for model development.

The selection of AAC, DDE, and UniRep features to represent important biological features associated with virulence proteins in fungi was made because these features provided complementary information regarding composition (AAC), sequence order (DDE), and evolutionary (UniRep) properties of the proteins of interest.

Feature fusion

To improve the predictive performance of ML and DL models, various feature representations were integrated. The amino acid composition (AAC) encodes the overall frequency of individual residues in a protein sequence, while the dipeptide deviation from expected mean (DDE) characterizes statistical deviations in dipeptide frequencies. In contrast, UniRep embeddings offer high-dimensional representations learned from large-scale protein data, capturing intricate evolutionary and contextual patterns. These distinct features provide complementary insights—ranging from basic compositional profiles and local sequence dependencies to deep, abstract representations.

By merging features from AAC, DDE, and pre-trained UniRep embeddings, feature fusion is accomplished as shown in Table 7. Concatenation serves to integrate these traits into a single, all-inclusive representation, improving their capacity to capture crucial biological data.

Feature integration was performed by aligning datasets on a shared identifier (protein_ids), ensuring accurate correspondence across records. Merging was carried out using Python's pandas library, followed by deduplication of label fields and standardization. Any missing entries introduced during the merge were imputed with zeros to maintain consistency across samples. This fusion of complementary representations enables models to leverage global composition, local sequence patterns, and deep learned structural or contextual information simultaneously, improving predictive performance. This approach enabled the construction of several comprehensive feature combinations as detailed in Table 7.

Data splitting

Following feature extraction and class balancing, the dataset was partitioned into three distinct subsets to facilitate model training, testing, and independent evaluation. Specifically, 70% of the data was allocated to the training set, which was used to develop ML and DL models. The test set, comprising 20% of the data, was used to evaluate model performance during development. An additional 10% was reserved as a validation set to assess the model's ability to generalize to unseen protein sequences, ensuring the robustness and reliability of the predictive framework.

To maintain equal representation of positive (virulent) and negative (non-virulent) classes, the dataset was stratified during the splitting process. Each class was separated and independently split using stratified sampling to preserve label distribution across subsets. After stratification, samples were randomly shuffled and recombined to form the final training, test, and validation sets. This careful partitioning strategy ensured that model evaluation remained unbiased while avoiding potential overfitting, thereby supporting reliable performance estimation.

Moreover, prior to dataset partitioning, duplication removal was performed to eliminate possibilities of bias in the similarities between additional sequences from the training, testing, and validation datasets.

Creation of an independent dataset for blind testing and generalizability assessment

In order to evaluate the prediction performance of our RF classifier, we prepared an independent dataset with sequences from a fungal species that was not included in the primary dataset (Candida spp.). PHI-Base, Victors, and DFVF were the same databases from which VFs from Neosartorya fumigata (Aspergillus fumigatus) were extracted as a positive, non-redundant dataset. This independent set, which is completely blind, was created exclusively to test the model's capacity to predict entirely new sequences. This enables us to evaluate the model's generalization ability to organisms that were not component of the training, test, or validation set. Random non-virulent protein sequences from the same fungal species were picked to create the negative dataset for independent set so as to match the 73 virulent proteins. The dataset was balanced to contain 73 virulent (positive) and 73 non-virulent (negative) sequences using a random undersampling technique. This dataset was subjected to the same feature extraction and preprocessing procedure and then split into three subsets: training (70%), testing (20%), and validation (10%). An impartial evaluation of the model's performance on the independent dataset was made possible by this comprehensive approach. Subsequently, five-fold cross-validation was implemented in conjunction with independent set to gauge the predictive ability of the RF model.

An independent dataset of data from other fungi, such as Neosartorya fumigata, not included in the main Candida dataset provided biological separation between the two datasets and decreased the chance of bias caused by a shared evolutionary history between the Candida and Neosartorya datasets.

Model building/classification

To evaluate the predictive ability of the extracted features, we trained classifiers using both deep learning (DL) and conventional machine learning (ML) algorithms. Specifically, we evaluated Deep Neural Networks (DNN), Multi-Layer Perceptrons (MLP), Random Forest (RF), and Artificial Neural Networks (ANN). The RF approach was selected because of its established ability to handle high-dimensional feature spaces while reducing the potential for overfitting through ensemble learning.41

The extracted features, including AAC, DDE, and UniRep, were used to train four different models (RF, ANN, MLP, and DNN), and their performances were compared to identify the most effective approach for distinguishing virulent from non-virulent fungal proteins. Given the limited availability of dedicated computational tools for predicting fungal virulence factors in human-pathogenic fungi, multiple ML and DL models were evaluated to assess their suitability for this task. Comparing DL architectures (MLP and DNN) with conventional ML approaches (RF and ANN) provided insight into the relative performance of these approaches for fungal VF classification.

The RF algorithm was implemented using 100 decision trees, bootstrap aggregation, and random feature selection at each split. Parallel computation was enabled across all available processor cores, and a random seed of 42 was used to ensure reproducibility.

The ANN model consisted of an input layer, two fully connected (FC) layers, and an output layer for VF prediction. To mitigate overfitting, the first FC layer included 64 neurons with ReLU activation and a 30% dropout rate. The second FC layer comprised 32 neurons with ReLU activation and a 20% dropout rate. The output layer contained a single neuron with a sigmoid activation function to predict probabilities for binary classification.

The MLP model for VF prediction consisted of three fully connected (FC) layers followed by an output layer. The first FC layer contained 64 neurons with ReLU activation, the second contained 32 neurons with ReLU activation, and the third contained 16 neurons with ReLU activation. The output layer contained a single neuron with a sigmoid activation function to predict probabilities for binary classification.

The DNN architecture used for VF prediction consisted of four fully connected (FC) layers followed by an output layer. To reduce overfitting and stabilize activations, the first FC layer contained 64 neurons with ReLU activation, batch normalization, and a 40% dropout rate. The second FC layer contained 32 neurons with ReLU activation, batch normalization, and a 20% dropout rate. Similarly, the third and fourth FC layers each contained 16 neurons with ReLU activation, batch normalization, and dropout rates of 10% and 20%, respectively. The output layer contained a single neuron with a sigmoid activation function to generate binary classification probabilities.

All models were trained using binary cross-entropy loss and optimized using the Adam optimizer. Hyperparameters were fine-tuned using the validation set, and the performance of each model was evaluated using the test set.

The architectures of the RF, ANN, MLP, and DNN models were designed to account for the high dimensionality of the feature space while maintaining a balance between model complexity and predictive performance. Hyperparameters, including the number of neurons, learning rate, batch size, and number of training epochs, were optimized using the validation set. To reduce overfitting, independent datasets for training, validation, and testing were used for model development and evaluation. In addition, comparative experiments were performed using individual and fused feature sets to assess the contribution of individual features to overall model performance and robustness.

Evaluation metrics

In addition to accuracy and AUC, additional metrics, including precision, recall (sensitivity), F1-score, and Matthews correlation coefficient (MCC), were used to provide complementary insights into model performance, particularly in the context of class imbalance. These metrics measure the correctness of positive predictions and the ability of the model to identify the minority class, represented by virulent proteins. False-positive predictions were also analyzed to assess potential downstream biological implications and to determine whether model predictions were consistent with biologically relevant patterns. MCC was included as a robust metric for evaluating binary classification performance, particularly for imbalanced datasets.

To thoroughly evaluate predictive performance, we used an independent test set and several standard evaluation metrics. These included recall or sensitivity (SN), specificity (SP), accuracy (ACC), F1-score, and Matthews correlation coefficient (MCC). The mathematical definitions of these metrics are given below:

SN=TPTP+FN

SP=TNTN+FP

ACC=TP+TNTP+FP+TN+FN

F1score=2×TP2TP+FP+FN

MCC=(TP×TN)−(FN×FP)(TP+FN)×(TN+FP)×(TP+FP)×(TN+FN)

Here, TP (true positives) represents the number of virulent proteins correctly predicted by the model, while TN (true negatives) refers to correctly predicted non-virulent proteins. FP (false positives) and FN (false negatives) denote incorrect predictions of non-virulent and virulent proteins, respectively.

Sensitivity refers to the model’s ability to correctly identify virulent proteins, and specificity measures its ability to correctly identify non-virulent proteins. Together with accuracy, F1 score, and MCC, these metrics provide a comprehensive and balanced understanding of the classifier's performance.

Quantification and statistical analysis

The performance of all machine learning models (RF, ANN, MLP, and DNN) was evaluated using five-fold cross-validation. For each model, the mean ± standard deviation of the performance metrics across the five folds was calculated. Model performance was assessed using accuracy, precision, recall, F1-score, and Matthews correlation coefficient (MCC). Independent test datasets were further used to evaluate the generalizability of the best-performing model.

No hypothesis-testing statistical analyses (such as t-tests or ANOVA) were performed because the objective of this study was to assess predictive model performance rather than compare experimental groups. Therefore, statistical significance indicators (e.g., asterisks) were not applicable.

Although VirulentPred, EffectorP, and FungalRV are commonly cited prediction tools, they were developed for bacterial virulence proteins, fungal effectors, or adhesin prediction and are not specifically intended for predicting the broad spectrum of virulence factors in human fungal pathogens. As a result, direct benchmarking against these methods using our dataset was not considered biologically appropriate. Instead, the performance of FunVFPred was evaluated through cross-validation and independent validation using experimentally characterized fungal virulence proteins.