Functional Prediction of Hypothetical Proteins from Shigella flexneri and Validation of the Predicted Models by Using ROC Curve Analysis
Article information
Abstract
Shigella spp. constitutes some of the key pathogens responsible for the global burden of diarrhoeal disease. With over 164 million reported cases per annum, shigellosis accounts for 1.1 million deaths each year. Majority of these cases occur among the children of the developing nations and the emergence of multi-drug resistance Shigella strains in clinical isolates demands the development of better/new drugs against this pathogen. The genome of Shigella flexneri was extensively analyzed and found 4,362 proteins among which the functions of 674 proteins, termed as hypothetical proteins (HPs) had not been previously elucidated. Amino acid sequences of all these 674 HPs were studied and the functions of a total of 39 HPs have been assigned with high level of confidence. Here we have utilized a combination of the latest versions of databases to assign the precise function of HPs for which no experimental information is available. These HPs were found to belong to various classes of proteins such as enzymes, binding proteins, signal transducers, lipoprotein, transporters, virulence and other proteins. Evaluation of the performance of the various computational tools conducted using receiver operating characteristic curve analysis and a resoundingly high average accuracy of 93.6% were obtained. Our comprehensive analysis will help to gain greater understanding for the development of many novel potential therapeutic interventions to defeat Shigella infection.
Introduction
Shigella, refers to a genus of gram-negative facultative anaerobes that belongs to members of the family Enterobacteriaceae and is the causative agent of shigellosis, a severe enteric infection, one of the most common causes of morbidity and mortality among children in developing nations. The Global Burden of Disease (GBD) classified Shigella as the second leading cause of diarrheal deaths on a global scale in 2015 [1]. Shigellosis leads to the recurrent passing of small, bloody mucoidal stools with synchronous abdominal cramps and tenesmus caused by ulceration of the colonic epithelium [2]. In malnourished children, Shigella infection may lead to a vicious cycle of further impaired nutrition, frequent infection and growth retardation resulting from protein loss enteropathy [3].
The Shigella genus is divided into four species: Shigella flexneri, Shigella boydii, Shigella sonnei, and Shigella dysenteriae. These are further classified into serotypes based on biochemical differences and variations in their O-antigen [4]. A total of 19 different serotypes of S. flexneri have been reported so far by various research groups [5]. Among the four Shigella species, shigellosis is predominantly caused by S. flexneri in the developing world especially in Asia, and is responsible for approximately 10% of all diarrheal episodes among children of <5 years [6]. Recent multicenter study in Asia revealed that the incidence of this disease might even exceed previous estimations, due to Shigella DNA being detectable in up to one third of the total culture negative specimens [7]. Currently, no effective vaccine with the ability to confer adequate protection against the many different serotypes of Shigella has been developed and made available. Existing antimicrobial treatments are becoming compromised in terms of efficacy due to increased antibiotic resistance, soaring cost of treatment, and persistence of poor hygiene and unsanitary conditions in the developing world.
A particular study conducted on numerous isolates of Shigella collected over a time span of 10 years, multi-drug resistance (MDR) were found to be exhibited by 78.5% of the isolates. 2% of the isolates were found to harbor genetic information capable of conferring resistance to azithromycin, a final resort antimicrobial agent for shigellosis [8]. On the other hand, a recent whole genome analysis of a particular strain of S. flexneri revealed 82 distinct chromosomal antibiotic resistance genes while successive re-sequencing platforms elucidated several distinct single nucleotide polymorphisms that contributed to eventual MDR [9]. Therefore, the development of new drugs has risen to become a subject of immense magnitude to not only shorten the medication period but also to treat MDR shigellosis. The genome sequence of S. flexneri serotype 2a strain 2457T, available in the NCBI database consists of 4,599,354 bp in a single circular chromosome containing 4,906 genes encoding 4,362 proteins and has G + C content of 50.9% [10]. Among these, the functions of 674 proteins have not been experimentally determined till date and are termed as hypothetical proteins (HPs). A HP is one that has been predicted to be encoded by an identified open reading frame, but for which there is a lack of experimental evidence [11]. Nearly half of the proteins in most genomes belong to the class of HPs and this class of proteins presumably have their own importance to complete genomic and proteomic platform of an organism [12, 13]. Precise annotation of the HPs of particular genome leads to the discovery of new structures as well as new functions, and elucidating a list of additional protein pathways and cascades, thus completing our incomplete understanding on the mosaic of proteins [13]. HPs may possibly play crucial roles in disease progression and survival of pathogen [11, 14]. Furthermore, novel HPs may also serve as markers and pharmacological targets for development of new drugs and therapies [15]. Functions of HPs from several pathogenic organisms have been already reported using a plethora of sequence and structure based methods [14, 16, 17].
Functional annotation of HPs utilizing advanced bioinformatics tools is a well-established platform in current proteomics [18]. Cost and time efficiency of these methods also favoring their preference over contemporary in vitro techniques [19]. In this study, we have used several well optimized and up to date bioinformatics tools to assign functions of a number of HPs from the genome of S. flexneri with high precision [20]. Functional domains were considered as the basis to infer the biological functions of HPs in this case. The receiver operating characteristic (ROC) analysis [21] was used for evaluating the performance of bioinformatics tools executed in our study. We also measured the confidence level of the functional predictions on the basis of bioinformatics tools employed during the course of the investigation [22]. We believe that this analysis will expand our knowledge regarding the functional roles of HPs of Shigella and provide an opportunity to unveil a number of potential novel drug targets [17].
Methods
The computational algorithm used for this study has been illustrated in Fig. 1. The entire work scheme has been divided into three phases namely, phase I, II and III. Phase I involves the characterization and sequence retrieval of the HPs, following the analysis of the S. flexneri genome. Phase II comprises of the annotation of various functional parameters using well optimized series of tools. The probable functions of the characterized HPs were predicted by the integration of various functional predictions. In phase III, an approach was made for systematic performance evaluation of various bioinformatics tools used in this study. In this case, S. flexneri protein sequences with known function were used as control. Finally, expert knowledge was applied for annotation of HPs at a considerable degree of confidence.
Computational algorithm used for annotating function of 39 hypothetical proteins (HPs) from Shigella flexneri. The framework has been divided into three phases: PHASE I, sequence retrieval from online databases; PHASE II, the extensive analysis of sub-cellular localization, physicochemical parameters, virulence, function and domain present in HPs; PHASE III, assessment of the predicted functions using the protein with known function from S. flexneri and reliable prediction of possible functions of HPs.
Phase I
Accession of genome and sequence retrieval
Complete genome sequence of S. flexneri 2a str. 2457T was retrieved from NCBI database (http://www.ncbi.nlm.nih.gov/genome/) and was found to code for a total of 4,362 proteins (accessed July 5, 2017). Fasta sequences of the complete coding sequence of 682 proteins, characterized as HPs were retrieved from UniProt (http://www.uniprot.org/). Finally, a total of 674 proteins were retained for downstream analysis following exclusion of duplicates.
Analysis of the conserved domains
Domains are often identified as recurring (sequence or structure) units, and can be thought of as distinct functional and/or structural units of a protein. During molecular evolution, it is assumed that domains may have been utilized as building blocks and have encountered recombination to modulate protein function [23]. A domain or fold might also exhibit a higher degree of conservancy when compared with the entire sequence [24].
In our study, five bioinformatics tools namely: CDD-BLAST (Conserved Domain Database-Basic Local Alignment Search Tool) [25–27], PFAM [28], HmmScan [29], SMART (Simple Modular Architecture Research Tool) [30], and SCANPROSITE [31] were used. These tools are able to search for the defined conserved domains in the targeted protein sequences and further assist in the classification of putative proteins in a particular protein family. HPs analyzed by five aforementioned function prediction web tools revealed the variable results when searched for the conserved domains in hypothetical sequences. Therefore, different confidence levels were assigned on the basis of collective results of these web-tools. One hundred percentage confidence level was considered upon obtaining the same results from the five distinct tools. Finally, we obtained 39 such proteins from 674 primary collected proteins, which were taken for further analysis (Supplementary Table 1).
Phase II
Physicochemical characterization
Theoretical physiochemical parameters such as molecular weight, isoelectric point, aliphatic index, instability index and grand average of hydropathicity (GRAVY) of these HPs were analyzed using ProtParam server of the Expasy tools (http://web.expasy.org/protparam/). Results of this analysis have been listed in Supplementary Table 2.
Determination of sub-cellular localization
For the identification of a protein as a drug or vaccine candidate, determination of the sub-cellular localization of the protein becomes particularly important. Surface membrane protein can be served as a potential vaccine target while cytoplasmic proteins may act as promising drug targets [32]. We used CELLO [33], PSORTb [34], and PSLpred [35] for the denotation of sub-cellular localization of the query proteins. TMHMM, SOSUI, and HMMTOP were applied for the prediction of query proteins for being a membrane protein, based on Hidden Markov Model [36–38]. SingnalP 4.1 [39] was used to predict the signal peptide and SecretomeP 2.0 [40] were utilized for the identification of proteins involved in non-classical secretory pathway. Results of these predictions are summarized in Supplementary Table 3.
Functional prediction of the query proteins
Various tools were used for precise functional assignments of all 39 HPs from S. flexneri (described in Table 1) such as CDD-Blast, Pfam, HmmScan, SMART, Scanprosite, MOTIF [41], INTERPROSCAN [42], CATH [43], SUPERFAMILY [44], and Protonet [45]. Results of these analyses have been outlined in Supplementary Tables 4 and 5.
List of annotated Hconf proteins from Shigella flexneri
The computational prediction of the structure of a protein from its amino acid sequences greatly facilitates the subsequent prediction of its function [46]. An online server PS2-v2 (PS Square version 2) [47], a template based method were used to predict the structure of the HPs. The modeling of proteins using this online server further substantiated the function of HPs. Besides, PFP-FunDSeqE [48] has been used to elucidate the protein fold patterns based on a combination of functional domain information and evolutionary information (Table 2).
Different types of folds identified in Shigella flexneri
Virulence factors analysis
Virulence factors (VFs) are described as potent targets for developing drugs because it is essential for the severity of infection [49]. VICMpred [50] and Virulentpred [51] tools were employed to predict VFs from protein sequences with an accuracy of 70.75% and 81.8%, respectively.
Functional protein association networks
The function and activity of a protein are often modulated by other proteins with which it interacts. Therefore, understanding of protein-protein interactions serve as valuable leads for predicting the function of a protein. In this investigation, we had employed STRING (Search Tool for the Retrieval of Interacting Genes/Proteins, https://string-db.org/) [52] to predict protein interactions partners of HPs. To predict functional association, only highest confidence score partner proteins were chosen in this study.
Phase III
Performance assessment
The predicted functions of HPs from S. flexneri and the accuracy of associated tools were validated using the ROC curve analysis. In this analysis, the diagnostics efficacy is evaluated at six levels where 1 and 0 classified as true positive and true negative respectively as binary numerals. In addition, the integers (2, 3, 4, and 5) were used as confidence ratings for each case. The ROC curves were carried out using 25 S. flexneri proteins with known function as control and were compared with the results obtained for the 39 HPs (Supplementary Tables 6 and 7). The results were submitted to web-based calculator for the ROC curves [53] in “format 1” form and the program thereby calculated the ROC curves. The results were expressed in terms of accuracy (Ac), sensitivity (Se), specificity (Sp) and the area under the curve (AUC) [54]. The average accuracy of the employed pipeline was found 93.6% (Table 3, Fig. 2).
List of accuracy, sensitivity, specificity, and ROC area of various bioinformatics tools used for predicting function of Hconf proteins from Shigella flexneri obtained after ROC analysis
Receiver operating characteristic (ROC) curve (area under the curve [AUC] = 0.774) for average accuracy of prediction.
Results and Discussion
Sequence analysis
Sequences of all the 674 HPs were analyzed for identification of the functional domains using five bioinformatics tools namely CDD-BLAST, Pfam, HmmScan, SMART, and SCANPROSITE. If the given five tools indicated the same domains for a protein, we considered it as 100% confidence level. In our study, all the five tools mentioned above revealed 39 such proteins and hence were grouped together. Only these HPs having 100% confidence level were considered for further analyses and termed as highly confident (Hconf) proteins. From the rest of the 635 proteins, no specific conserved domains were found for a total of 257 proteins. For other HPs (n = 378), specific domains were identified using several of these tools. To know accurate function of these proteins further studies are required.
The function of each of these 39 Hconf were successfully assigned by using different online tools, listed in Table 1. All sequence analyses were compiled and categorized into various functional classes constituting 9 enzymes, 10 binding proteins, 4 transporters, 4 lipoproteins, 6 which are involved in various cellular processes, while 6 proteins were predicted to exhibit miscellaneous functions (Fig. 3). Various functional classes of these classified Hconf proteins are described below.
Hypothetical proteins classified into different groups based on their functions.
Enzymes
Enzymes are key players in many leading biochemical processes in the living system and may facilitate the survival of pathogens in the host and making it viable for the course of infection. A total of 9 proteins out of 39 (23%) of our annotated Hconfs were characterized as enzymes. Among these, two proteins were characterized as transferases, among which, WP_000301054.1 is a lipopolysaccharide kinase (Kdo/WaaP), involved in the formation of outer membrane (OM) of gram negative bacteria and is encoded by the Waap gene. The OM protects cells from toxic molecules and is important for survival during infection and is required for virulence of the pathogen [55]. According to reports made by Delucia [55], the depletion of WaaP gene was seen to halt the growth of the bacteria suggesting that WaaP is essential to produce the full-length lipopolysaccharide, recognized by the OM [49]. Therefore, WaaP may result in a potent target for the development of novel antimicrobial agents. The other transferase, protein WP_000778795.1 was found to consist of an acetyltransferase (GNAT) domain that uses acetyl coenzyme A (CoA) to transfer an acetyl group to a substrate, a reaction implicated in various functions for the development of antibiotic resistance of bacteria [56].
Three enzymes were predicted to be hydrolases, which plays key role in the invasion of the host tissue and evading the host defense mechanism and are thus associated with various VFs [57]. For instance, WP_005051685.1 marks the lysin-like motif/peptidase family M23, is found in proteins from viruses, bacteria, fungi, plants and mammals. It is present in bacterial extracellular proteins including hydrolases, adhesins and VFs such as protein A from Staphylococcus aureus. We report WP_001295493.1 protein as the endoribonuclease/YjgF family active on single-stranded mRNA that inhibits protein synthesis by cleaving mRNA [58]. YjgF family members are enamine/imine deaminases that hydrolyze reactive intermediates released by pyridoxal phosphate-dependent enzymes, including threonine dehydratase [59]. It has also been reported in the inhibition of transaminase B in Salmonella [60].
Among the other enzymes predicted, there has been two isomerase and one lyase enzyme. WP_001247854.1 constitutes the toprim (topoisomerase-primase),a catalytic domain involved in breakage and rejoining of DNA strand [61]. WP_001205243.1 marks the Xylose isomerase-like TIM barrel involved in the myo-inositol catabolism pathway [62]. Lyases also play a key role in bacterial pathogenesis due to their involvements in various biosynthesis processes. WP_000943980.1 was found to demonstrate synthase activity that causes hydrolysis of ATP with the formation of an amide bond between spermidine and the glycine carboxylate of glutathione. In the pathogenic trypanosomatids, this reaction is the penultimate step in the biosynthesis of the antioxidant metabolite, and is a resounding target for target mediated drug design [63]. The WP_000454701.1 protein was found to be a cystathionine b-lyase, an enzyme which forms the cystathionine intermediate in cysteine biosynthesis and may be considered as the target for pyridiamine anti-microbial agents [64].
Binding proteins
Ten proteins annotated as binding proteins among which 1 RNA binding, 3 protein binding, 3 lipid binding, 1 metal binding, 1 peptidoglycan binding, and 1 adhesion protein have been predicted. WP_000132640.1 protein was predicted to be SymE (SOS-induced yjiW gene with similarity to MazE). It has been reported to involve in inhibiting cell growth, decrease protein synthesis and increase RNA degradation and thus exhibit a vital role in the survival and propagation of pathogen in the host [65, 66]. Despite not manifesting any functional homology with other type I toxin proteins, SymE belongs to the type I toxin-antitoxin system. Its function resembles that of type II toxins such as MazF, which is able to perform the cleavage of mRNA in a ribosome independent manner. However, SymE shares homology to the AbrB-fold superfamily proteins such as MazE, which acts as transcriptional factors and antitoxins in various type II TA modules [67]. It seems probable that SymE has evolved into an RNA cleavage protein with toxin-like properties from a transcription factor or antitoxin [66]. In our study, we reported WP_000003197.1 as von Willebrand factor with a type A domain which has been reported responsible for various blood disorders [68–70]. The association of type A domain makes it liable to be involved in various significant activities such as cell adhesion and immune defense [71]. On the other hand, WP_000755956.1 has been predicted to belong to the band-7 protein family that comprises of a diverse set of membrane-bound proteins characterized by the presence of a conserved domain [72]. The exact function of this domain is not known, but concurrent reports from animal and bacterial stomatin-type proteins demonstrate the ability of binding to lipids and in the assembly of membrane-bound oligomers that form putative scaffolds [73]. We have also predicted WP_001269672.1 and WP_000749269.1 as the lipid binding domain called lipopolysaccharide (LPS)-assembly lipoprotein LptE and the YceI-like domain respectively. The LPS transport machinery is composed of LptA, LptB, LptC, LptD, and LptE. LptE forms a complex with LptD, which is involved in the assembly of LPS in the outer leaflet of the OM [74]. This OM is an effective permeability barrier that protects the cells from toxic compounds, such as antibiotics and detergents, thus conferring the bacteria with the capability to adapt and consequently inhabit several different and often hostile environments. Among the binding prtoteins, WP_000266171.1 was found to be a tetratricopeptide repeat containing protein which is involved in protein-protein interactions and thus plays an important role in virulence [75].
Cellular processes/regulatory proteins
A total of 6 HPs have been predicted to be involved in various cellular and regulatory mechanisms, which are vital cognates in the pathogenesis of S. flexneri and thus can be treated as possible drug targets [76]. For example, WP_000189314.1 predicted to be a member of the GIY-YIG family involved in many cellular processes including DNA repair and recombination, transfer of mobile genetic elements, and restriction of incoming foreign DNA [77, 78].
WP_001387238.1 and WP_001297375.1 have been found to be RadC-like domain belonging to the JAB superfamily of metalloproteins [79]. In most instances, this domain shows fusions to an N-terminal Helix-hairpin-Helix (HhH) domain and may also be function as a nuclease [79]. WP_000848528.1 has been predicted to be a leucine-zipper found in the enterobacterial OM lipoprotein LPP [80]. It is likely that this domain is involved in protein-protein interaction via subsequent oligomerization. WP_000597196.1 and WP_048814497.1 have been respectively found to be a Glycine zipper 2TM domain found in the Rickettsia genus and leucine-rich repeat involved in a variety of biological processes, including signal transduction, cell adhesion, DNA repair, recombination, transcription, RNA processing, disease resistance, apoptosis, and the immune response [81].
Lipoprotein
Bacterial lipoproteins are a set of membrane proteins with many different functions. Due to this broad-ranging functionality, these proteins have a considerable significance in many phenomena, from cellular physiology through cell division and virulence [82]. Lipoprotein of gram-negative bacteria is essential for growth and division [83]. In our analysis, we report a total of 4 lipoproteins from the group of HPs predicted in this study. It has been also revealed that lipoproteins may function as vaccines [82]. The knowledge of these facts may be utilized for the generation of novel countermeasures to bacterial diseases [82].
Transport
In our findings, we report the prediction of HP WP_000070107.1 to be a member of the ATP-binding cassette superfamily, largest of all protein families with a diversity of physiological functions [84]. It has recently been identified that these proteins may be involved in virulence and are essential for intracellular survival of pathogens [85]. We have found protein WP_001296791.1 to be an auto-transporter of the YhjY type involved in DNA repair [86]. Protein WP_001238362.1 has been found to exhibit the function of transport of nutrients, control of cell regulation, pheromone transport, cryptic coloration and in the enzymatic synthesis of prostaglandins. An example a protein with such function is the retinol-binding protein 4, which transfers retinol from liver to peripheral tissues [87].
Other proteins
Six HPs have been predicted to exhibit miscellaneous functions where most of them are protein with unknown function. Among them, WP_001382892.1 and WP_000691930.1 have been preicted to be domains of unknown function and are found in a number of bacterial proteins. WP_001125713.1 has been found to be YcgL domain with conserved class of small proteins widespread in gammaproteo bacteria. This group of proteins contain a 85-residue domain of unknown function and two alpha-helices and four beta-strands in the sequential arrangement [88]. We have also predicted WP_001237866.1 and WP_005049020.1 as YecR and YehR like family of lipoproteins found in bacteria and viruses and are functionally uncharacterized.
Virulent proteins
Gram-negative bacteria undergo frequent genomic alterations and consequent evolutions thus increasing their virulence inside the host environment [89]. We have found 2 HPs that showed positive virulence scores servers among the Hconf proteins. These have been listed in Supplementary Table 8. It had already been hypothesized that targeting VF provides a better therapeutic intervention strategy against bacterial pathogenesis [89]. Predicted HPs having virulent characteristics thus provide powerful target-based therapies for the mitigation of an existing infection and are further considered as an adjunct therapy to existing antibiotics, or potentiators of the host immune response [90].
ROC curve
The average accuracy of the employed pipeline was identified 93.6% in our analysis which indicated that outcomes of the functional annotation of HPs were predicted with a high degree of confidence. We have also found sensitivity of 100% and specificity 64.3% for the tools used in this study. Finally, area under the curve was found to be 0.774. AUC is an effective way to summarize the overall diagnostic accuracy of the test. It takes values from 0 to 1, where 0.7 to 0.8 is considered acceptable.
Conclusion
Using an innovative in silico approach, all 674 HPs from S. flexneri were primarily analyzed and then using the ROC analysis and confidence level measurements of the predicted results the functions of the 39 HPs were precisely predicted with a reasonably high degree of confidence and thereby were successfully characterized. Following this, the validation of the functions of these proteins were carried out by using different approaches including structure based PS2-v2 server, sub-cellular localization and physicochemical parameters. These are important for distinguishing the HPs from the rest of the protein. The protein-protein interaction also gave insights in elucidation of the involvement of such proteins in various metabolic pathways. Moreover, some virulence proteins had also been detected which are essential for the survival of this pathogen. This in silico approach for functional annotation of the HPs can be further utilized in drug discovery for characterizing putative drug targets for other clinically important pathogens. The outcomes of ROC analysis indicated high reliability of bioinformatics tools used in this study. Hence, the functional annotation of HPs is reliable and can be further utilized for other experimental research.
Supplementary Information
Acknowledgments
The authors are grateful to core donors which provide unrestricted support to icddr,b for its operations and research. Current donors providing unrestricted support include: Government of the People’s Republic of Bangladesh; Canadian International Development Agency (CIDA), Swedish International Development Cooperation Agency (Sida), and the Department for International Development, UK (DFID). We gratefully acknowledge these donors for their support and commitment to icddr,b’s research efforts.
Notes
Author’s contributions
Conceptualization: MAG
Data curation: MAG, SM, SMF, MGK
Formal analysis: MAG, SM
Funding acquisition: MAG, MRI, HR, SD
Methodology: SM, MGK, PP, MRI, HR, SD
Writing – original draft: MAG, MM, TA
Writing – review & editing: MAG, MM, TA
Conflict of interest
No potential conflicts of interest relevant to this article was reported.