Investigation of hub genes and their nonsynonymous single nucleotide polymorphism analysis in Plasmodium falciparum for designing therapeutic methodologies using next-generation sequencing approach.

BACKGROUND: Incidences of resistance to current drugs by Plasmodium is increasing, hence, it is necessary to investigate and explore new drug targets to combat malarial disease.
OBJECTIVE: Analysis of the transcriptome sequence information to characterize hub genes and their nonsynonymous single nucleotide polymorphisms (nsSNPs) to derive therapeutic objectives for Plasmodium falciparum.
MATERIALS AND METHODS: Differentially expressed genes between Ring and other stages of P. falciparum were identified using Cufflinks tool. Using DAVID and KAAS programs, the gene ontology and pathway analysis were performed. The networks of protein-protein interaction (PPI) were developed by Search Tool for the Retrieval of Interacting Genes/Proteins and Cytoscape, and the node degree in the network was calculated by using Network Analyzer, and MCODE plugins of Cytoscape. SIFT, PROVEAN, and PredictSNP programs were used to study the genetic variations, which affect protein functions.
RESULTS: A list of 4196 nonredundant genes was used for functional annotation cluster analysis, and 8 significant hub genes have been picked from the PPI network using MCODE plugins of Cytoscape. Various nsSNPs were identified in these 8 hub genes and were investigated both for its native and mutant stage for solvent accessibility and alteration in secondary structure protein residues.
CONCLUSION: Hub genes identified in this study serve as potential targets to develop therapy to suppress the pathogenic action of P. falciprum through experimental techniques. Copyright:

Introduction

Malaria, a deadly infectious disease, is caused by intracellular single-celled parasites belongs to the genus Plasmodium. According to World Malaria Report 2018 by World Health Organization (WHO), about 219 million malaria cases and 435,000 deaths were reported in 87 nations in 2017. In the human genome, it is recognized as one of the most powerful evolutionary selection. Plasmodium falciparum, Plasmodium vivax, Plasmodium knowlesi, Plasmodium malariae, and Plasmodium ovale are different species that cause disease in humans. The most severe malignant malaria particularly in children below 5 years is P. falciparum.[1] In India, almost 95% population lives in endemic areas of malaria, and 80% of recorded malaria in India is limited to 20% of the population areas who lives in rural, hilly, remote, or difficult-to-reach areas.[2] Traditional first-line therapies such as chloroquine and pyrimethamine/sulfadoxine have lost their efficacy in most countries, which resulted in the development of new and more effective anti-malarial medicines, like artemisinin-based combination therapy.[3] Although efforts were made to develop vaccines and medicines to fight malaria, there is still a problem with vaccine escape and drug resistance.[4] What remains in dearth is how plasmodium genetic variation can result in drug resistance or can provide new drug targets. It has been shown that genetic variation and recombination accelerate antigen heterogeneity, immune escape, the production of anti-malarial drug resistance, and comparison of entire genomes may aid in these efforts.[5] RNA-seq is a method that can use next-generation sequencing to analyze the quantity and sequences of RNA in a sample. Knowing the transcriptome is essential to linking the genome data to its functional protein expression. RNA-seq technology enables accurate gene isoform identification, translocation events, nucleotide mutations, and posttranscriptional modifications.

In P. falciparum, the genetic diversity exists in the form of single nucleotide polymorphism (SNPs), microsatellite repeats, insertions, deletions, and a variety of gene duplication. Several studies have been reported to investigate the SNPs of P. falciparum.[67] Nevertheless, to the best of our acquaintance, no research to date has used the impact of deleterious SNPs and solvent accessibility and secondary structure change of protein residues at the native and mutant level in P. falciparum.

López-Barragán et al. 2011 have sequenced seven bidirectional and four strand specific libraries for the analysis of gene expression and antisense transcripts. RNA-Seq data sets of seven bidirectional libraries were used in the study. Therefore, network analysis was performed to recognize the main hubs from the Ring (R), Schizont (Sc), gametocyte stage V (GV), gametocyte stage II (GII), early trophozoite (ET), Late trophozoite (LT), and Ookinete (Oo) stages of P. falciprum. In these prospective hubs, we have examined SNPs which can be proposed as key areas for vulnerability to affect the function of protein. The areas projected in this research can be further targeted to develop therapy to suppress the pathogenic action of P. falciprum through experimental therapeutic techniques such as gene knockout method and gene targeting.

Materials and Methods

Datasets

The RNA-Seq dataset was obtained from NCBI website, https://www.ncbi.nlm.nih.gov/sra/for all stages with accession number SRP009370[8] for analysis. Seven bidirectional reads from the 3D7 parasite have been taken for further investigation. The non-synonymous SNPs (nsSNPs) information of the P. falciparum genes were collected from PlasmoDB database (https://plasmodb.org/plasmo/). The protein sequence was obtained from the UniProtKB database (https://www.uniprot.org/) in the FASTA format.

Sequence quality control and preprocessing

The quality of sequence data obtained from high throughput sequencing pipelines were checked using FASTQC tool (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) was followed to clean the reads with low base-call quality using a quality filter tool.

Sequence reads alignment and transcript assembly

The reference genome of P. falciparum 3D7 (PlasmoDB version 7.1) was used to map the sequence reads with TopHat tool version: 2.0.14 (https://ccb.jhu.edu/software/tophat/).[9] The minimum anchor length was seven base pairs for reads present at each side and a maximum size of intron 800 bp. The output of TopHat was filtered to keep only reads mapped from Ring to gametocyte stages with 0 mismatches and up to 1 mismatch in Oo stage to maximize the accuracy. The Cufflinks (http://cole-trapnell-lab.github.io/cufflinks/)[9] tools were used to further evaluate the matched reads using a multifasta files (plasmodium_falciparum.fa), which improves the reliability of abundance of transcript by detecting bias and an algorithm for correction. The relative affluence of transcripts has been depicted as fragments/kilobase of exon/million fragments mapped (FPKM) and the fold change of the FPKM value between Ring (R) and other stage was reported. Poorly expressed genes were removed from the dataset by eliminating genes with FPKM value <2 in all the stages. A list of nonredundant genes which differentially expressed was created from the RNA-Seq data by combining all identified genes of Ring and other stages and the duplicate gene were removed. Fold change ratio was calculated for R versus ET, R versus LT, R versus Sc, R versus GII, R versus GV and R vs Oo groups. Further, genes which are expressed differentially were manually checked in the PlasmoDB[10] and UniProt database.[11]

Functional analysis and pathway mapping

DAVID web server[12] was used to identify and select significantly enriched gene ontology (GO) terms and pathways. The functional annotation was determined with DAVID program (http://david.abcc.ncifcrf.gov/home.jsp). Those terms with count number of ≥5 genes, and P ≤ 0.05 was chosen for analysis. In DAVID, the terms GO cellular component (CC), biological process (BP), and molecular function (MF) were used to classify improved biological topics in lists of genes which differentially expressed. The KEGG Automatic Annotation Server (KAAS) has been used to map the pathway.[13] Using the best single-directional hit method for orthology assignment, the amino acid sequences of genes which were up regulated and under regulated were submitted as input to the KAAS server. KAAS offers functional gene annotation in the database KEGG GENES through a similarity search tool of BLAST for a manually curated orthology group sets. For datasets mapped to one of reference pathways of KEGG, KAAS assigned a KEGG Orthology (KO) number to the genes.

Investigation of protein–protein interactions

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) program was used to study the functional relationships between differentially expressed genes (DEGs).[14] The GO terms – MF, CC, BP, and pathways of KEGG were used to create the interaction network using reference organism P. falciparum. The physical and functional interactions, coexpression, colocation, pathways, predicted interactions, and protein domain similarity are included in the network relationship between genes. The network has been filtered by eliminating all interactions with weights below 0.1.

Identification of hub genes

Cytoscape 3.4.0,[15] a data integration and network visualization bioinformatics package, has been used to identify hub genes by measuring distribution of node degrees via Network Analyzer plugin. Clusters were created using Cytoscape molecular complex detection (MCODE) plug-in. It detects regions, which are highly interconnected in the network.[16] The statistical criteria for MCODE are as follows: K score = 2, Cutoff degree = 2, Cutoff node score = 0.2 as the default. In the current network, top 5 genes with the largest distribution have been regarded as hubs.

Prediction of deleterious nonsynonymous single nucleotide polymorphisms of hub genes

PlasmoDB database was used to extract the information of SNPs for hub genes identified for P. falciparum. Deleterious coding nsSNPs were predicted by using PROVEAN, SIFT, and PredictSNP. SIFT and PROVEAN are tools to determine how amino acid substitution affects protein function in case the score is lower than threshold value. The PROVEAN tool was used to perform BLAST hits clustering. For each supporting sequence, a delta alignment score was calculated. To calculate the final PROVEAN score, average scores are generated across clusters.[17] Deleterious is considered a score of −2.5 or higher, when anything below this cut-off rating, which has an effect neutral. SIFT tool based on the physical properties of amino acids and homology of sequences determines whether substitution of an amino acid can influence the function of protein.[18] The sequence tool SIFT makes SIFT predictions for a particular sequence of proteins in FASTA format. The sequence of protein queries and interest substitutions of nsSNPs and hub genes with default parameters have been submitted tohttp://sift.bii.a-star.edu.sg/www/SIFT_seq_submit2.html. SIFT program predicts substitutions with values <0.05 to be detrimental. PredictSNP[19] has been explicitly designed to combine the projected outcomes of several tools to form a prediction of consensus. For the predictions, prediction tools use a list of variants in the protein sequence as input.

Protein solvent accessibility and prediction of secondary structure

To predict protein accessible surface area (ASA), surface accessibility and secondary structure, NetSurfP program,http://www.cbs.dtu.dk/services/NetSurfP was used. The NetSurfP simultaneously predicts accuracy for each prediction by calculating the Z-score. It uses two types of neural networks, the first type networks are based on secondary structure predictions and sequence profiles, with two outputs with respect to buried or exposed and in combination with sequence profiles, the other networks use these outputs as inputs and are trained to assess the relative surface exposure of each amino acid residues.[20] The normal and predicted SNP sequences in FASTA format have been submitted for prediction to the NetSurfP. For prediction of secondary structure and ASA, protein encoding gene in normal and its predicted SNP substitutions have been uploaded individually to NetSurfP web server. Microsoft Excel®2016 has converted the most predicted secondary structure probabilities from NetSurfP to a single letter code representing helical (H), β-strand (E), and coil (C).

Results

The RNA-Seq dataset for all the stages of P. falciprum such as Ring SRR364849, ET SRR364848, LT SRR364847, Schizont SRR364843, GII SRR364840, GV SRR364838, and Oo SRR364834 has downloaded from NCBI SRA for analysis. SRA Toolkit modules, vdb-validate and fastq-dump were used to validate the integrity of downloaded SRA data and to convert SRA data into fastq format, respectively. FASTX Toolkit, a quality filter tool, was used to clean the reads with low base-call quality. Bowtie index was built from a set of DNA sequences of P. falciparum chromosomes. RNA-seq reads for every stage were mapped for P. falciparum genome using TopHat program. The output of TopHat was filtered to keep only reads, which were mapped with 0 mismatches from ring to gametocyte stages, and up to 1 mismatch in Oo stage to maximize the accuracy.

Prediction of differentially expressed genes

Transcripts were assembled and analyzed in FPKM for their relative expression levels by Cufflinks tool after sequencing reads mapped to the TopHat reference genome. Between Ring and other stages, fold change analysis was done to compare the gene expression. The data with FPKM values equivalent to zero were removed, and remaining values were subjected to further analysis. As a result, 7517 genes from ring (R), 6799 genes from ET, 7482 genes from LT, 5102 genes from schizont (Sc), 8731 genes from gametocyte stages (GII), 8831 genes from gametocyte stages (GV), and 5155 genes from Oo stages were identified. Then, six groups – RvET, RvLT, RvSc, RvGII, RvGV, and RvOo were created, and the common gene in each group has been identified [Table 1].

Table 1

Genes which are differentially expressed in different stages of Plasmodium falciparum

Stages	Common genes	DEGs (≥2 fold)
		Upregulated	Downregulated	Total
RvET	4402	731	1711	2442
RvLT	4354	1895	901	2796
RvSc	4022	1249	1686	2935
RvGII	3917	897	1910	2807
RvGV	2988	1116	1064	2180
RvOo	4009	2070	825	2895

DEGs=Differentially expressed genes

The common genes between RvET-4402; RvLT-4354; RvSc-4022; RvGII-3917; RvGV-2988; and RvOo– 4009 have been identified. The fold change in each group was estimated using these common genes, which was defined as the FPKM value ratio of RvET, RvLT, RvSc, RvGII, RvGV, and RvOo groups. To classify the DEGs, less expressed genes with FPKM value <2 have been removed from the dataset in all the stages. On the basis of above criteria, there are 2442 DEGs between RvET, 2796 between RvLT, 2935 between RvSc, 2807 between RvGII, 2180 between RvGV, and 2895 between RvOo groups were sorted out for the analysis [Table 1]. The complete workflow for the analysis is depicted in Figure 1.

Figure 1

The flow chart depicting the overall methodology adopted in this study

Functional annotations and pathway analyses

The functional annotation cluster analysis was conducted by DAVID tool on the DEGs. The GO terms – BP, MF, and CC were used for interpretation and only those terms with count number of ≥5 genes and P ≤ 0.05 were selected. Five top GO terms with significant P values for each group from functional analysis are presented in Table 2. Only those gene IDs which were mapped via DAVID tool is used for further study. From this data, it is clear that the GO terms are enriched in RvET2209, RvLT2546, and RvGV1990 genes, which represent the functions necessary for host cell plasma membrane, antigenic variation, and pathogenesis [Table 2]. The enhanced GO terms in RvSc2735 including functions related to pathogenesis, single organismal cell–cell adhesion, receptor activity, and cell adhesion molecule binding [Table 2], and RvGII2594 include functions related to single organismal cell-cell adhesion, pathogenesis, receptor activity, and cell adhesion molecule binding. Similarly, functions-related single organismal cell–cell adhesion, pathogenesis, and cell adhesion molecule binding are enriched in RvOo2726 samples [Table 2].

Table 2

Five top improved gene ontology terms detected by DAVID in genes expressed differentially between a) RvET2209, b) RvLT2546, c) RvSc2735, d) RvGII2594, e) RvGV1990, and f) RvOo2726

Gene Ontology (GO) term	Total number of genes	P value
a) RvET2209
GO: 0009405~pathogenesis	60	1.50E-11
GO: 0004872~receptor activity	55	2.49E-10
GO: 0020033~antigenic variation	112	2.61E-10
GO: 0020002~host cell plasma membrane	122	7.27E-09
GO: 0020013~modulation by symbiont of host erythrocyte aggregation	55	1.09E-07
b) RvLT2546
GO: 0016337~single organismal cell-cell adhesion	55	1.25E-08
GO: 0009405~pathogenesis	57	8.41E-06
GO: 0004872~receptor activity	53	1.69E-05
GO: 0020033~antigenic variation	111	2.22E-04
GO: 0020002~host cell plasma membrane	117	0.00279306
c) RvSc2735
GO: 0016337~single organismal cell-cell adhesion	51	1.41E-04
GO: 0020035~cytoadherence to microvasculature, mediated by symbiont protein	44	4.92E-04
GO: 0050839~cell adhesion molecule binding	46	0.001862072
GO: 0009405~pathogenesis	53	0.008046292
GO: 0004872~receptor activity	49	0.050202845
d) RvGII2594
GO: 0016337~single organismal cell-cell adhesion	50	7.85E-05
GO: 0020035~cytoadherence to microvasculature, mediated by symbiont protein	44	1.16E-04
GO: 0050839~cell adhesion molecule binding	45	6.05E-04
GO: 0004872~receptor activity	49	0.00977134
GO: 0009405~pathogenesis	50	0.015754517
e) RvGV1990
GO: 0009405~pathogenesis	56	3.05E-09
GO: 0004872~receptor activity	53	3.57E-09
GO: 0020002~host cell plasma membrane	114	3.83E-09
GO: 0020033~antigenic variation	106	2.00E-08
GO: 0020013~modulation by symbiont of host erythrocyte aggregation	53	5.99E-07
f) RvOo2726
GO: 0016337~single organismal cell-cell adhesion	47	9.00E-04
GO: 0020035~cytoadherence to microvasculature, mediated by symbiont protein	39	0.008244114
GO: 0020030~infected host cell surface knob	38	0.008860222
GO: 0050839~cell adhesion molecule binding	40	0.04538152
GO: 0009405~pathogenesis	47	0.056291162

GO=Gene ontology

The KAAS has been used to map the pathways of DEGs in all six stages. DEGs amino acid sequences in FASTA format were submitted to the KAAS server. As a result, 277 pathways were predicted for RvET2209, 283 pathways for RvLT2546, 280 pathways for RvSc2735, 285 pathways for RvGII2594, 229 pathways for RvGV1990, and 272 pathways for RvOo2726 were recognized. The top 10 KEGG pathways for each six categories are shown in Table 3, and Supplementary Table 1 provides a complete list of pathways. It was observed from Table 3 that most DEGs have been linked with important biological processes, many of which are classified as metabolic pathways, secondary metabolite production pathways, ribosome, or being involved in biosynthesis of antibiotics.

Table 3

Ten top KEGG pathways in each of the six categories

	Pathway name	Number of mapped Genes
a) RvET
ko01100	Metabolic pathways	102
ko01110	Biosynthesis of secondary metabolites	44
ko03010	Ribosome	43
ko01130	Biosynthesis of antibiotics	32
ko01120	Microbial metabolism in diverse environments	20
ko00230	Purine metabolism	17
ko01200	Carbon metabolism	17
ko00240	Pyrimidine metabolism	14
ko05144	Malaria	13
ko04141	Protein processing in endoplasmic reticulum	12
b) RvLT
ko01100	Metabolic pathways	146
ko01110	Biosynthesis of secondary metabolites	63
ko01130	Biosynthesis of antibiotics	43
ko00230	Purine metabolism	33
ko01120	Microbial metabolism in diverse environments	27
ko05169	Epstein-Barr virus infection	26
ko00240	Pyrimidine metabolism	26
ko01200	Carbon metabolism	26
ko03030	DNA replication	24
ko00190	Oxidative phosphorylation	23
c) RvSc
ko01100	Metabolic pathways	129
ko03010	Ribosome	56
ko03040	Spliceosome	55
ko01110	Biosynthesis of secondary metabolites	53
ko01130	Biosynthesis of antibiotics	35
ko00230	Purine metabolism	32
ko00240	Pyrimidine metabolism	30
ko03008	Ribosome biogenesis in eukaryotes	29
ko01120	Microbial metabolism in diverse environments	29
ko03013	RNA transport	27
d) RvGII
ko01100	Metabolic pathways	137
ko01110	Biosynthesis of secondary metabolites	56
ko03040	Spliceosome	51
ko01130	Biosynthesis of antibiotics	39
ko03010	Ribosome	38
ko01120	Microbial metabolism in diverse environments	29
ko03008	Ribosome biogenesis in eukaryotes	27
ko05169	Epstein-Barr virus infection	27
ko03013	RNA transport	27
e) RvGV
ko01200	Carbon metabolism	25
ko01100	Metabolic pathways	96
ko01110	Biosynthesis of secondary metabolites	41
ko03010	Ribosome	35
ko03040	Spliceosome	34
ko01130	Biosynthesis of antibiotics	25
ko03008	Ribosome biogenesis in eukaryotes	20
ko00230	Purine metabolism	20
ko01120	Microbial metabolism in diverse environments	19
ko03013	RNA transport	17
ko05016	Huntington’s disease	17
f) RvOo
ko01100	Metabolic pathways	132
ko01110	Biosynthesis of secondary metabolites	60
ko03010	Ribosome	47
ko01130	Biosynthesis of antibiotics	44
ko03040	Spliceosome	36
ko01120	Microbial metabolism in diverse environments	31
ko00230	Purine metabolism	31
ko01200	Carbon metabolism	28
ko03013	RNA transport	27
ko00240	Pyrimidine metabolism	24

STRING program has been used to investigate DEGs interaction. Only those genes showing significant interactions with weights higher than 0.4 were selected for network analysis. A network between DEGs was constructed for all the six groups, namely RvET2209, RvLT2546, RvSc2735, RvGII2594, RvGV1990, and RvOo2726. The six interaction networks resulted from STRING were then subjected to Cytoscape. The network consists of 1647 nodes and 23970 edges for RvET2209, 1969 nodes and 36152 edges for RvLT2209, 2236 nodes and 48514 edges for RvSc2735, 2010 nodes and 30113 edges for RvGII2594, 1473 nodes and 17986 edges for RvGV1990, and 2095 nodes and 37531 edges for RvOo2726. All genes in the network are represented in circles and their interactions represent edges. The interaction networks of all six groups are given in Supplementary Figure 1.

Further, all six networks have been analyzed using Network Analyzer and MCODE modules available in Cytoscape. Network Analyzer is calculating the node degrees in the network, whereas the MCODE module is creating the clusters in the network. The higher node degrees were regarded to be more significant genes and were referred as hub genes. The top 8 MCODE seed are PF3D7_1126700, PF3D7_0508100, PF3D7_1306000, PF3D7_1439500, PF3D7_0324900, PF3D7_1234300, PF3D7_1207100, and PF3D7_0705600. The interaction between these hubs and their first neighbors was presented as Figure 2.

Figure 2

(a) PF3D7_0324900: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster. (b) PF3D7_0508100: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster. (c) PF3D7_0705600: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster. (d) PF3D7_1126700: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster. (e) PF3D7_1207100: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster. (f) PF3D7_1234300: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster. (g) PF3D7_1306000: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster. (h) PF3D7_1439500: Interaction between hub gene and their first neighbors. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The hub is shown as purple diamond in each cluster

Single nucleotide polymorphism analysis

A SNP is a significant source of variance in a genome. SNPs may result in affecting protein function by decreasing protein solubility or by destabilizing structure of protein.[21] The PlasmoDB database has been used to retrieve the nsSNPs and SNPs for the identified hub genes. The PF3D7_0324900 have highest SNPs and nsSNPs information, i.e., 2416 SNPs and 1606 nsSNPs from the database among all hub genes; whereas PF3D7_1439500, PF3D7_0508100, PF3D7_1306000, PF3D7_0705600, PF3D7_1207100, PF3D7_1126700, and PF3D7_1234300 have 171, 128, 101, 62, 38, 28, and 11 nsSNPs, respectively. The total number of nsSNPs and SNPs recognized for the hub genes are illustrated in Figure 3.

Figure 3

Distribution of total number of single nucleotide polymorphisms and nonsynonymous single nucleotide polymorphisms identified for Plasmodium falciparum hub genes

Analysis of nsSNPs for protein function

Deleterious coding nsSNPs were predicted using SIFT, PROVEAN, and PredictSNP tools. SIFT is a sequence homology-based tool used to classify substitutions for amino acids. The SIFT tool predicts whether substitution of an amino acid can affect protein function for a given FASTA protein sequence or not. The SIFT predicts substitutions with values <0.05 to be deleterious. The SIFT sequence tool predicted 60, 52, 38, 33, 28, 19, 7, and 2 positions to be affect protein function for PF3D7_1439500, PF3D7_1306000, PF3D7_0508100, PF3D7_0705600, PF3D7_0324900, PF3D7_1207100, PF3D7_1126700, and PF3D7_1234300, respectively. To predict the final PROVEAN score, a delta alignment score is calculated for every supporting sequence and mean value across the clusters. By default, a score of −2.5 or above of it is taken as deleterious, while anything short of this cutoff rating has been considered as neutral. In the hub genes, the PROVEAN protein tool predicted 40, 4, and 2 positions to be deleterious for PF3D7_0324900, PF3D7_1306000, PF3D7_1207100 and 1 position for PF3D7_0705600 and PF3D7_1234300, respectively. However, PROVEAN tool not found any deleterious mutation in PF3D7_1439500, PF3D7_0508100, and PF3D7_112670. PredictSNP unambiguously designed to combine outcomes of several methods, mostly to annotate disease-variant relationships. According to PredictSNP tool, there are 18, 16, 11, 9, 6, 4, 2, and 1 mutations to be deleterious for PF3D7_1306000, PF3D7_0324900, PF3D7_1439500, PF3D7_0705600, PF3D7_1207100, PF3D7_0508100, PF3D7_1126700, and PF3D7_1234300, respectively.

It was confirmed and verified at least by two tools used above in the study that 25 positions for PF3D7_0324900, 20 for PF3D7_1306000, 11 for PF3D7_1439500, 10 for PF3D7_0705600, 6 for PF3D7_1207100, 4 for PF3D7_0508100, 2 for PF3D7_1126700, and 1 for PF3D7_1234300 were identified to influence the function of protein [Supplementary Table 2].

Solvent accessibility and secondary structure prediction

Further, the secondary structures and solvent accessibility of the hub genes was investigated using NetSurfP-1.1. The SNPs which were predicted by at least two tools to have a negative effect on the function of protein were used for the analysis of secondary structure and solvent accessibility. Moreover, the data were screened by selecting the residues, which showed change in ASA ≥10 Š2 from buried state to exposed state and also exposed state to buried state and its secondary structure change. The information related to ASA and secondary structure are shown individually in Supplementary Table 3.

In PF3D7_0324900, N615K mutation showed an ASA change to exposed state from buried but T1745P, S1747R, S2024R, P2026S and E2065K mutations showed the opposite change to buried state from exposed state. There was a change in most of the conformations from Coil to Alpha-Helix, due to these mutations [Table 4]. Due to mutations Y779D and Y862N, an ASA change to exposed state from buried was shown in PF3D7_1306000 whereas D1113Y showed an inverse change. Conformations in Y779D and D1113Y were changed from Alpha-Helix to Coil, while in Y862N, Beta-strand was changed to Coil. Mutations in the residue position N186H, S312N, N334S, Y338N, S348C, Y403S, S445L, and R570G in PF3D7_1439500 show a change to exposed state from the buried state and mutations S313N, T411I, and N746 K show a change from buried state in the majority of the cases. In most cases, it shows changes in secondary structure from coil to helix and coil to beta-strand. In PF3D7_0705600, mutations show almost the same change to exposed state from buried state and vice versa. It shows changes in secondary structure from helix to coil and helix to beta-strand in most cases and also from coil to helix in some cases. Due to mutation D377H, an ASA change primarily from buried state to exposed state was shown in PF3D7_1207100 whereas N235K, R277W, P528S, and N661Y show an inverse change. In R277W mutation, both F274 and E275 show change in C to H conformation. D377H shows change in conformation mainly from H to C. In P528S mutation, both L337 and Y340 show change from C to H conformation. In N661Y mutation, L337 show change from C to H and I410 show an opposite change in conformation. Mutations in the residues I800T, I1103K, S1411P and Y1474H, in PF3D7_0508100 show mostly from exposed to buried state. In most of the cases, it shows changes in secondary structure from helix to coil and vice versa. Some also change from helix to beta-strand conformations. In PF3D7_1234300, S434Y mutation, I137 show change in conformation from coil to helix. In PF3D7_1126700, no significant change was detected [Table 4].

Table 4

NetSurfP results showing change in accessible surface area from buried to exposed state and vice versa and also show change in their secondary structure

Mutation	Class change		Conformation change
	B-E	E-B	C-E	C-H	E-C	E-H	H-C	H-E
PF3D7_0324900
N615K	5	3	1	-	-	-	-	-
T1745P	4	4	-	1	-	-	-	-
S1747R	1	2	-	1	-	-	-	-
S2024R	1	2	-	1	-	-	-	-
P2026S	-	3	-	1	-	-	-	-
E2065K	1	2	1	-	1	-	-	-
PF3D7_1306000
Y779D	1	-	-	-	-	-	1	-
Y862N	1	-	-	-	1	-	-	-
D1113Y	2	1	-	-	-	-	1	-
PF3D7_1439500
N186H	83	25	5	8	2	-	2	-
S312N	61	54	4	4	3	-	2	-
S313N	64	69	7	6	2	-	1	1
N334S	85	25	2	6	2	-	3	-
Y338N	69	54	7	1	2	-	2	-
S348C	78	33	4	7	1	-	2	-
Y403S	70	68	2	12	2	-	2	-
T411I	59	66	7	7	3	-	2	-
S445L	69	59	4	9	5	-	2	-
R570G	67	58	5	8	5	-	2	-
N746K	69	73	9	8	2	-	3	-
PF3D7_0705600
N124Y	3	2	1	1	-	-	-	-
E623V	24	24	2	2	1	-	8	4
A626D	20	26	2	2	1	-	7	4
D645Y	22	27	2	5	1	-	6	5
T688R	20	22	2	6	1	-	8	5
T852S	14	13	1	1	2	-	1	1
PF3D7_1207100
N235K	6	15	1	3	-	-	1	-
R277W	3	9	-	2	-	-	-	-
D377H	16	8	-	-	-	1	2	-
P528S	4	3	-	2	-	-	-	-
N661Y	5	8	-	1	-	-	1	-
PF3D7_0508100
I800T	30	45	-	4	-	-	7	2
I1103K	43	51	3	5	2	1	10	2
S1411P	47	58	2	8	-	2	15	3
Y1474H	55	64	2	10	-	-	13	3
PF3D7_1234300
S434Y	3	4	-	1	-	-	-	-

B=Buried, E=Exposed, H=Alpha-Helix, E=Beta-strand, C=Coil

Further, a database has been developed to show the analysis done by PROVEAN, SIFT, PredictSNP, and NetSurfP software, which is available online at URL www.cdkd.org/pfsnp/.

Database was developed using PHP and JavaScript with user-friendly search environment using a range of options, such as simple searches and advanced searches. Users can search easily with any of the three fields such as gene name, UniProt ID, or SNP ID. Gene name with corresponding residues is available in advance search. Users can search for particular residues in any gene or genes available in this database (currently the hub genes). For displaying the SNP data, two-step approach is exploited by these two search methods. The first step is to display Gene name, UniProt ID, SNP ID, AA change, PredictSNP cuttoff, Confidence, PROVEAN Score, SIFT Score, and cutoff. The second step will reveal the details of the secondary structure, SNP ID, Residue N, Location, ASA N, Class N, SS N, Residue SNP, SNP, ASA SNP, Class SNP, and SS SNP. Users can also customize their search results by choosing any field like Class change, SS change, ASA change or Residue N or any combination for a particular gene. The relevant information is displayed dynamically. In the help section, all headers are defined and hyperlinked to this web portal from the search pages. In addition, every SNP ID is directly linked to its relevant entries in the PlasmoDB database.

Discussion

The most virulent pathogen of malaria and malaria mortality worldwide is P. falciparum. To control the disease in parasites of malaria, the study of genetic variation is of practical importance. A nonredundant list of 4196 genes was used for the functional annotation cluster analysis by the DAVID tool. The most important enriched GO terms identified in these genes by DAVID functional cluster analysis consist of functions required for the host cell plasma membrane, antigenic variation, and pathogenesis [Table 2]. The node degree was calculated for each gene in the network to explore the functional roles of genes involved in different processes and interaction networks were created using Network Analyzer and MCODE plugins of Cytoscape. PF3D7_0324900, PF3D7_1306000, PF3D7_1439500, PF3D7_0705600, PF3D7_1207100, PF3D7_0508100, PF3D7_1126700, and PF3D7_1234300 genes were considered as hubs genes in the network constructed from 4196 P. falciparum genes by Network Analyzer and MCODE plugins of Cytoscape. Hub genes are considered functionally important because these are highly interconnected with nodes in a system. Therefore, these can act as putative targets for drug designing.

In this study, a number of nsSNPs have been identified for nearly all hubs, but few have had an impact on protein functions. PF3D7_0324900, PF3D7_1306000, PF3D7_1439500, PF3D7_0705600, PF3D7_1207100, PF3D7_0508100, PF3D7_1126700, and PF3D7_1234300 genes were annotated by at least two tools to affect protein function. The NetSurfP web server has been used to predict the protein secondary structure and surface accessibility in normal for each gene and its predicted SNP substitutions. The erythrocyte membrane protein 1 of P. falciparum, PfEMP1 (PF3D7_0324900) is considered as potential hub in this study, which mediates the attachment of infected erythrocytes to a range of host cells in the vascular lining during the blood stage of the infection with malaria. PfEMP1 is inserted into the RBC membrane and then laterally transferred to the preformed knob structures in the central region. Knob assembly may result in new ways of inhibiting PfEMP1 presentation on infected RBCs.[22] Strategies for overcoming PfEMP1 antigenic diversity would provide an exciting new opportunity for the development of malaria vaccines.[23] PF3D7_1306000 is a conserved protein of Plasmodium with unknown function and considered as hub gene. This gene was found as essential along with druggability index 0.5 in Tropical Disease Research (TDR) Targets database (https://tdrtargets.org/). Another hub CCAAT-binding transcription factor or oocyst rupture protein 2 (ORP2) (PF3D7_1439500) was involved in sporozoite egress. Sporozoite egress of the oocyst can be blocked by removing the N-terminal histone fold domain of ORP2.[24]

Hub gene RNA helicase, (PF3D7_0705600) play various roles, including the cell growth and development. Helicases are significant unwinding enzymes that are needed in the malaria parasite for nearly all the nucleic acid metabolism. RNA helicases could be used as reasonable targets to develop new antiparasite therapies and solve the problem of drug resistance.[25] Small subunit rRNA processing factor, (PF3D7_1207100) involved in the maturation of SSU-rRNA from tricistronic rRNA transcript was also found as a hub gene. It is a protease group of enzymes, which play key roles in the development and invasion of parasites. The ability to design particular protease inhibitors makes them promising objectives for drugs.[26] SET domain protein, (PF3D7_0508100), which enables histone-lysine N-methyltransferase activity was also observed as hub gene. Two types of protein methyltransferase enzymes (PMTs) are present in eukaryotic cells: lysine specific and arginine specific. They were both linked with a variety of diseases including neurodegenerative and inflammatory diseases, and cancer. PMT enzymes have emerged as a target class against human disease for drug discovery.[27]

The gene PF3D7_1126700 that was also observed as a hub, codes for autophagy-related protein 23. The ATG18 autophagy-related protein controls the biogenesis of apicoplast in apicomplex parasites and decline of ATG18 in P. falciparum showed in delayed death.[28] DNA polymerase epsilon subunit B (PF3D7_1234300) is involved in the DNA-dependent DNA replication and enables DNA-directed DNA polymerase activity. The main function of Pol epsilon is to extend the leading strand synthesis during replication.[29]

A total of 6 mutations (N615K, T1745P, S1747R, S2024R, P2026S, and E2065K) were identified for PF3D7_0324900, 3 mutations (Y779D, Y862N and D1113Y) for PF3D7_1306000, 11 mutations (N186H, S312N, S313N, N334S, Y338N, S348C, Y403S, T411I, S445 L, R570G and N746K) for PF3D7_1439500, 6 mutations (N124Y, E623V, A626D, D645Y, T688R and T852S) for PF3D7_0705600, 5 mutations (N235K, R277W, D377H, P528S and N661Y) for PF3D7_1207100, 4 mutations (I800T, I1103K, S1411P and Y1474H) for PF3D7_0508100, and only one mutation for PF3D7_1234300. These mutations brought a change in accessible surface area from buried to exposed state and vice versa and also change in their secondary structure as observed through NetSurfP tool.

A computational approach to systematically analyze nsSNPs was undertaken in this study to predict deleterious mutations. SNP within the proteins significantly affects stability of structure of a protein and its function. The role of nsSNPs in understanding the functional effects of mutations that may cause changes in amino acids in hub proteins was also examined. Further, determination of nsSNPs that interfere with the function of protein and cause a disease needs to be determined.[30]

Conclusion

In this study, a comprehensive approach has been used to derive potential therapeutic targets for P. falciparum using RNA-seq dataset. The differential expression of genes, functional, and pathway enrichment analysis of P. falciparum was apprized in detail. The present study results suggested that PF3D7_0705600, PF3D7_1207100, PF3D7_0508100, PF3D7_1126700, and PF3D7_1234300 hub genes serves as putative targets for drug designing. These hub genes are showing less mutation and no similarity with human proteins. These genes act as lysine methyltransferases, transcription, translation, RNA splicing, and other important cellular pathways in P. falciparum. Moreover, the study also provides clusters of hub genes and their network pathways analysis which can be used further studies to devise a therapeutic target that stabilizes their gene expression. In addition, nsSNPs and their functional impact of these hub genes were also calculated.

Financial support and sponsorship

Nil.

Conflicts of interest

There are no conflicts of interest.

Network of all six groups. A red node indicates query genes and the genes predicted by Search Tool for the Retrieval of Interacting Genes/Proteins are shown in cyan. The top hubs are shown as purple diamonds. (a) Interaction networks between R-ET differentially expressed genes and additional related genes. (b) Interaction network between R-LT differentially expressed genes and additional related genes. (c) Interaction network between R-Sc differentially expressed genes and additional related genes. (d) Interaction network between R-GII differentially expressed genes and additional related genes. (e) Interaction network between R-GV differentially expressed genes and additional related genes. (f) Interaction network between R-Oo differentially expressed genes and additional related genes

List of KEGG pathways of Ring verses other stages (R-ET, R-LT, R-Sc, R-GII, R-GV, and R-Oo) of P. falciparum

The nonsynonymous single nucleotide polymorphisms that were predicted to affect protein function in Plasmodium falciparum by using SIFT and PredictSNP tools

NetSurfp results showing change in accessible surface area from buried to exposed state and vice versa and also change in their secondary structure