Identifying disease-associated signaling pathways through a novel effector gene analysis

Background Signaling pathway analysis methods are commonly used to explain biological behaviors of disease cells. Effector genes typically decide functional attributes (associated with biological behaviors of disease cells) by abnormal signals they received. The signals that the effector genes receive can be quite different in normal vs. disease conditions. However, most of current signaling pathway analysis methods do not take these signal variations into consideration. Methods In this study, we developed a novel signaling pathway analysis method called signaling pathway functional attributes analysis (SPFA) method. This method analyzes the signal variations that effector genes received between two conditions (normal and disease) in different signaling pathways. Results We compared the SPFA method to seven other methods across 33 Gene Expression Omnibus datasets using three measurements: the median rank of target pathways, the median p-value of target pathways, and the percentages of significant pathways. The results confirmed that SPFA was the top-ranking method in terms of median rank of target pathways and the fourth best method in terms of median p-value of target pathways. SPFA’s percentage of significant pathways was modest, indicating a good false positive rate and false negative rate. Overall, SPFA was comparable to the other methods. Our results also suggested that the signal variations calculated by SPFA could help identify abnormal functional attributes and parts of pathways. The SPFA R code and functions can be accessed at https://github.com/ZhenshenBao/SPFA.


INTRODUCTION
Recently developed high-throughput functional genomics technologies have generated large amounts of experimental disease data and detected new biological information. Challenge for biologists is understanding the biological behaviors of disease cells using both newly generated disease data and existing biological knowledge. Signaling pathway analysis methods are used to better understand the biological behaviors of disease cells.
The understanding of biological behaviors of disease cells benefits to understand the pathological scenery and treatment. Over-representation analysis (ORA) based methods were initially presented as signaling pathway analysis methods to help biologists identify over-represented pathways from a list of relevant genes produced from experimental data. ORA-based methods merely count the number of differentially expressed genes in specific functional category gene sets such as the Gene Ontology (GO) (Blake et al., 2013), the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2016), BioCarta (Nishimura, 2001), and Reactome (Joshitope et al., 2005). Then they determine significance of the overlaps via statistical tests such as Fisher's exact test. Many tools are based on this method including Onto-Express (Draghici et al., 2003;Khatri et al., 2002), Fisher (Khatri, Sirota & Butte, 2012), and the Gene Ontology Enrichment Analysis Software Toolkit (GOEAST) (Zheng & Wang, 2008). However, ORA-based methods only take into account large changes in individual genes that significantly affect pathways and they do not account for smaller changes in sets of functionally-related genes (i.e., pathways) capable of significant effects. Functional class scoring (FCS) based methods have been used to avoid this limitation of ORA-based methods. FCS-based methods take into account the coordinated gene expression changes in pathways, such as gene set enrichment analysis (GSEA) (Subramanian et al., 2005), gene set analysis (GSA) (Efron & Tibshirani, 2006), and mean-rank gene set enrichment tests (MRGSE) (Liu et al., 2008). However, ORA-based and FCS-based methods are both limited because they do not take into account the complex interactions between genes or the complex topology of pathways. To overcome this limitation, pathway-topology-based methods were proposed. Pathway-topology-based methods integrate the complex interactions between genes using pathway topology information, specifically KEGG signaling pathway information.
Signaling pathway impact analysis (SPIA), one of the most widely-used pathwaytopology-based methods, considers both the number of differentially expressed genes (DEGs) in a given pathway and the topology information of that pathway (Tarca et al., 2009). Many improved methods based on SPIA have been proposed. Li et al. (2015) developed a method called sub-SPIA, which used a minimum spanning tree way to prune signaling pathways and used the SPIA method to identify significant signaling subpathways (Li et al., 2015). Bao et al. (2016) developed two SPIA-based methods called PSPIA and MSPIA. These two methods replaced +1 or −1 interaction strength in SPIA with the interaction strength of the Pearson correlation coefficients and mutual information, respectively (Bao et al., 2016). There are different pathway-topology methods that make use of the topological information of signaling pathways. For instance, Gene Graph Enrichment Analysis (GGEA) uses prior knowledge derived from directed gene regulatory networks (Geistlinger et al., 2011). Liu, Xu & Bao (2019 used a subgraph method to take advantage of pathway topological information (Liu, Xu & Bao, 2019). ROntoTools introduced a term of perturbation factor by considering the type of interactions to take the pathway topology into consideration (Tarca et al., 2009;Voichita, Donato & Draghici, 2012). Sebastian-Leon et al. (2014) developed a method using topology to detect liner subpathways in a signaling pathway (Sebastian-Leon et al., 2014).
These methods still have their disadvantages. Pathway-topology-based methods do not consider the importance of genes in pathways. Gene-weight-based methods have been proposed to overcome this limitation. Pathway analysis with down-weighting of overlapping genes (PADOG) uses the frequency of a present gene in the analyzed pathways to improve gene set analysis (Tarca et al., 2012). Functional link enrichment of gene ontology or gene sets (LEGO) measures gene weights in a gene set according to its relative association with genes inside and outside the gene set in a functional association network (Dong et al., 2016). Fang et al. (2017) proposed an improved SPIA method called SPIA-IS that measured and assigned the importance as the average output degree of the gene in the pathway.
A signaling pathway is a cascade of molecular reactions that bring out the functional attributes (e.g., cell proliferation, apoptosis) associated with the biological behaviors of disease cells using effector genes. Effector genes receive signals without outputting signals to other genes in an individual signaling pathway (Sebastian-Leon et al., 2014). Diseases are always related to the abnormal signal that the effector genes receive. Therefore, the signal that the effector genes receive can be very different under disease and normal conditions. The limitation of the previously mentioned methods, including gene-weight-based methods, is that they do not consider the signal variations between disease and normal conditions. Additionally, the functional attributes in the same signaling pathway may be very different from one another, and can sometimes be opposites. For example, there are two opposite functional attributes on the axon guidance pathway: axon repulsion and axon attraction (see the hsa04360 pathway in the KEGG dataset). We cannot determine which functional attributes contribute more to the disease using most current pathway analysis methods. Furthermore, some pathways consist of several parts, each with very different contributions. For example, the Wnt signaling pathway is significant across different diseases and can be divided into three parts. Most existing pathway analysis methods cannot determine which part of the Wnt signaling pathway most significantly contributes to a specific disease.
We propose a new method that considers the signal variations between normal and disease conditions that effector genes received in pathways: the signaling pathway functional attributes analysis (SPFA) method. SPFA calculates the gene expression changes in a given pathway using an ORA method and then combines the ORA method results with the signal variation results under two conditions (normal vs. disease). The signal variations can help identify functional attributes and abnormal pathways. We tested the capabilities of the proposed signaling pathway analysis method on a series of real datasets using three parameters. We also showed that the two types of probabilities considered in this method were indeed independent. Ultimately, we verified the usefulness of the signal variations the effector genes received under two different conditions using the SPFA method.

Data sources and preprocessing
Signaling pathway analysis methods require two types of input: a collection of pathways and a list of genes or gene products with accompanying expression values across different samples between the compared phenotypes. We used the KEGG signaling pathway as it is the most common manually-curated signaling pathway used for pathway analysis. We downloaded 213 signaling pathways from the KEGG PATHWAY dataset.
We acquired 33 disease gene expression datasets from the KEGGdzPathwaysGEO R-package and KEGGandMetacoreDzPathwaysGEO R-packages (Table 1) (Tarca, Bhatti & Romero, 2013;Tarca et al., 2012). Each disease gene expression dataset was matched with a corresponding disease KEGG pathway. For example, a colorectal cancer dataset was associated with the colorectal cancer pathway (Tarca et al., 2012). The corresponding disease KEGG pathways were called target pathways. Three rules were used to select the gene expression datasets: 1. The dataset's DEGs were available. If no DEGs were selected, other comparable methods would return null results.
2. The results of these datasets could be analyzed. Pathway analysis result p-values of 1 could not be analyzed.
3. The target pathways of these datasets were KEGG pathways since we used KEGG pathways as examples.
DEGs were selected if they contained more than 200 genes with FDR adjusted p-values < 0.05. Otherwise, we selected more than 200 genes with original p-values < 0.05 and absolute log (fold change) > 1.5. If DEGs still less than 200 genes, we selected the top 1% of genes ranked by p-values as DEGs.

SPFA algorithm design
To assess the signal variations between two conditions (normal vs. disease) that the effector genes received from upstream genes, we calculated the sum of signal variations from all upstream genes to effector genes. Given an effector gene g e and an upstream gene g s , the signal variation from the gene g s to the effector gene g e can be defined as: where cor disease ðg s g e Þ and cor normal ðg s g e Þ refer to the Pearson correlation coefficient between the gene expression data of gene g s and gene g e in the disease and normal states, respectively. d se is the network distance between gene g s and gene g e . The Pearson correlation coefficient is always used in gene co-expression networks to represent the strength of interactions between two genes. The Pearson correlation coefficient can also be used to represent the strength of an interaction between two gene expression values. Studies have shown that the genetic regulatory patterns in signaling pathways between genes are different under normal and disease conditions (Jung, 2018). If the genetic regulatory pattern between the two genes changes, the signal transmitted between the two genes will be very different. Thus, we used the Pearson correlation coefficient to calculate the signal variations that the effector genes received from their upstream genes. However, if an upstream gene does not directly transmit a signal, the signal may be attenuated. Therefore, we used the network distance d se between gene g s and gene g e as a penalty coefficient. For each effector gene g i in a given pathway, the accumulated signal variations between normal and disease conditions that the upstream genes received (total s genes in the upstream of the gene g i ) were calculated using the formula: The accumulated signal variation ASVðg i Þ of the effector gene g i in a pathway can help us distinguish among the functional disease attributes. Effector genes with high ASVðg i Þ demonstrate that these functional attributes significantly contribute to their corresponding diseases.
For a given signaling pathway, the total accumulated signal variation ASV can be defined as: where k is the total number of effector genes in the given pathway. Ultimately, the probability P sd used to measure the signal variations between two conditions (normal vs. disease) that those effector genes received from genes upstream in a given signaling pathway P x is based on the pathway's ASVðP x Þ. The same number of genes as the one observed on the given signaling pathway are randomly selected from all genes (random gene IDs) and have any possible expression data in all samples in the range of the experimenter. Therefore, the observed signal variations were obtained by permuting the gene IDs 2000 times. ASV per ðP x Þwas the total accumulated signal variation of the given pathway P x obtained in the per-th time. The probability P sd ðP x Þ of the given pathway was calculated as: where I is a function that returns 1 when the argument is true and 0 when the argument is false. The probability P sd does not measure the gene differential expression in a given pathway. Thus, we had to combine the probability P sd with the probability P de which can measure the total gene differential expression in a given signaling pathway. The probability P de of a given pathway P x can be calculated through the following hypergeometric test: where the whole genome contains a total of m genes, n genes are the number of DEGs in the m genes, and the given pathway contains t genes and r DEGs.
The probability P sd uses the Pearson correlation coefficient of the two genes' expression data, but the probability P de uses the number of DEGs in a pathway. Thus, the two probabilities are independent of each other. The significance of the given pathway was calculated following the SPIA method which combines the probabilities P sd and P de (Tarca et al., 2009). The formulas are: where c is a product of P de and P sd . P is the combined probability of the signaling pathway.
Significantly enriched pathway analysis using SPFA The SPFA procedure identifies significantly enriched pathways in two steps (Fig. 1). The first step measures the total gene differential expression in the signaling pathways. DEGs need to first be identified from the gene expression datasets. Then the DEGs are mapped onto the signaling pathways. Finally, the signaling pathway p-values are calculated using a hypergeometric test. The second step is to measure the signal variations between the two conditions (normal vs. disease) that effector genes received from upstream genes in the signaling pathways. This is completed by: 1. Finding all effector genes in each signaling pathway.
2. Ascertaining all paths from the upstream genes to the effector genes in each signaling pathway. If a path exists between the upstream genes and effector genes, an interaction must exist between them. The path's network distances are used to weight the corresponding interactions.
3. Using the Pearson correlation coefficient absolute difference values between the disease and normal samples to calculate the signal variations of the corresponding interactions.
4. Using the network distance of each interaction to decrease their signal variations.
5. Calculating the accumulation of the signal variations between the effector genes and upstream genes for each effector gene.
6. Calculating the sum of the accumulations of all effector genes in each signaling pathway.
7. Evaluating the statistical significance of each pathway based on their score.
Ultimately, the results of the two steps are combined into one p-value. We used the FDR adjust method on the combined p-value to determine the significance of each signaling pathway. The pathways with the adjusted combined p-values smaller than a threshold value were considered as significant pathways.

The distribution of effector genes in the signaling pathways
Studying the signal variations between two conditions (normal vs. disease) that the effector genes received leads to a deeper understanding of the biological behaviors of disease cells. Effector genes are widely scattered throughout the signaling pathways. If a gene has no signal inputs in an individual signaling pathway, the gene is not considered an effector gene. The distribution of effector genes in each signaling pathway can be seen in Fig. 2. One hundred and ninety-five of the 213 signaling pathways contained effector genes.
There is no universally accepted technique for the validation of the results of pathway analysis methods. Current pathway analysis methods use the results of a very small number of datasets based on searching corresponding published life literature. This approach has its limitations. First, the number of datasets used is small. Second, authors often search their own, leading to biased results. Third, complex biological phenomena always directly or indirectly correspond to multiple signaling pathways. Tarca et al. (2012) compiled an objective and reproducible approach based on multiple datasets (Tarca et al., 2012). This approach avoided a biased literature search and required testing on a large number of different datasets (at least 10). In this work, we followed The second measurement was the median rank of the 33 disease target pathways. The higher ranked methods were more accurate. To further validate the different pathway analysis method results, we used a third measurement: the ratio of significant pathways (using a significance threshold of 0.05 of the adjusted p-value) in the 33 datasets. This measured the method's ability to control false positive and false negative rates.

The independence between the two probabilities
The two probabilities P de and P sd are theoretically independent under the null hypothesis. We verified their independence by calculating the squared correlation coefficient between the two probabilities using the 33 gene expression datasets (Table 2). Our results showed that the average squared correlation coefficient of the 33 datasets was R 2 ¼ 0:029. Only four of the 33 squared correlation coefficients were slightly higher than R 2 ¼ 0:09. These results indicated essentially no correlation between the two probabilities.

SPFA method performance
We compared SPFA with the other seven methods using three measurements: the median p-value of the 33 target pathways, the median rank of the 33 target pathways, and the ratio of significant pathways. The signaling pathways with adjusted p-values ≤ 0.05 were significant. When comparing the median rank of the 33 target pathways, SPFA ranked first (Fig. 3). When comparing the median p-value of the 33 target pathways, SPFA ranked fourth  . 4). Notably, the methods with the highest ranking in one measurement did not necessarily rank the highest in the others. This is because different measurements analyze different abilities. For example, MRGSE was first in median p-value but was sixth in median rank. Fisher was second in median p-value but ranked fourth in median rank.
To better compare SPFA's performance against the other methods, we added the ranks of Table 2 The squared correlation coefficients between the two probabilities using the 33 gene expression datasets. The four squared correlation coefficients which are slightly more than 0.09 are shown in bold. the median p-value and median rank values from each method together. We found that the combined value of SPFA and PADOG was the smallest (Table 3).

GEO ID
To further assess the performance of the eight methods, we collected the results from other general pathways typically associated with cancer using the 18 out of 33 datasets with a form of cancer in Table 4: Apoptosis and Pathways in cancer. When using the Apoptosis pathway and Pathway in cancer pathway instead of target pathways, SPFA's median ranks were both first, and the median p-values of MRGSE were also both ranked first. These results were in alignment with the target pathway results. However, when using the Apoptosis pathway and Pathway in cancer pathway instead of the target pathways, PADOG's median p-values were both ranked fifth. When using the Apoptosis pathway, SPFA's median p-value ranked third. When using the Pathway in cancer pathway, SPFA's median p-value ranked fourth. All these results suggest that SPFA had the best accuracy and a good sensitivity when compared with the other seven methods. Additionally, our results showed that SPFA's ratio of significant pathways was moderate, 0.16 (Fig. 5), compared to the others. MRGSE's ratio of significant pathways was almost 0.5, and it could be questioned whether a such number of pathways was realistic. GSA's ratio of significant pathways was lower than 0.05, and it reflected that the GSA method had a high false negative rate. The methods had a modest ratio of significant pathways indicated that the method had a modest false positive rate and a modest false negative rate. Thus, the discriminative ability of SPFA was good when compared with the other seven methods. In conclusion, our results strongly supported that SPFA was well-suited for signaling pathway analysis and confirmed previously reported results in Dong et al. (2016).

Sources of improvement for SPFA
The main source of improvement in SPFA is that it uses signal variations that effector genes received under normal and disease conditions. SPFA is compared to the simpler ORA-based method used to calculate the probability P de without accounting for signal variations (Fig. 6). As shown in Fig. 6, the ORA-based method has a higher (worse) rank and p-value than SPFA for the target pathways. Table 3 The combined rank values of the ranks in terms of the median p-values and the median ranks of target pathways of eight methods.

Ranks of the median p-values
Ranks of the median ranks Sum  Table 4 The results of other general pathways: apoptosis and pathway in cancer typically associated with cancer using the 18 out of 33 datasets with a form of cancer. For each pathway, the values for the type of methods with the smallest median p-values and ranks (strongest association with the phenotype) are shown in bold. Validating the correlation between diseases and the signal variations that effector genes received under two different conditions The Alzheimer's disease dataset GSE16759 included four disease samples and four normal samples (Juan et al., 2010). The Wnt signaling pathway was altered in 90% of the colorectal cancer samples (Galamb et al., 2008). We assessed the signal variations that effector genes received in the Wnt signaling pathway using the GSE4183 dataset (Fig. 7). The results of (Galamb et al., 2008) coincided with our signal variation results (Galamb et al., 2008) reported that overexpression of TNS1 could induce the activation of JNK (ENTREZID: 5599, 5601, and 5602). The signal variation that the effector gene ENTREZID: 5602 received ranked first in our results. Galamb et al. (2008) detected that RBMS1 is another overexpressed gene and modulator of c-myc (ENTREZID: 4609). c-myc can regulate cell cycles and cause cells to transform pathways. The signal variation that the effector gene ENTREZID: 4609 received ranked second in our results. Galamb et al. (2008) also identified that TCF4 is an overexpressed gene that can participate in the transcriptional regulation of genes associated with colon carcinogenesis. These colon carcinogenesis associated genes include c-myc (ENTREZID: 4609), cy-clin D1 (ENTREZID: 595), PPARδ (ENTREZID: 5467), and MMP7 (ENTREZID: 4316). The signal variations that these effector genes received ranked second, fourth, fifth, and sixth, respectively. Many pathways can be studied in colorectal cancer datasets. For example, the PI3K-Akt signaling pathway plays a critical role in the growth and progression of colorectal cancer (Johnson et al., 2010). The effector genes ENTREZID:596, ENTREZID:842, and ENTREZID:1027 have the highest signal variations and are linked to cell cycle progression and cell survival (Fig. 8). The GSE4183 dataset results further confirmed the role of this pathway in colorectal cancer development.

Pathway statistic
The Wnt signaling pathway is also closely related to the occurrence and development of Alzheimer's disease (Inestrosa et al., 2007). The signal variations that different effector genes received calculating based on the Alzheimer's disease dataset GSE16759 in the Wnt signaling pathway were shown in Fig. 9. The signal variations that the effector genes: ENTREZID: 595 and 896 received were considerably higher than the other effector genes in the Wnt signaling pathway. This result validated evidence of crosstalk between the Alzheimer's disease signaling pathway and the two effector genes' upstream genes in the Wnt signaling pathway.

DISCUSSION
Functional attributes (associated with biological behaviors of disease cells) are the responses that effector genes respond to the signal they received. Disease cells always have abnormal functional attributes. Thus, the signal that the effector genes received can be very different. However, no current pathway analysis method takes this factor into consideration. Most pathway analysis methods only include the activation and significance of pathways. Their results give us inadequate information on functional attributes that can help explain the biological behaviors of disease cells. Here, we proposed SPFA, a novel signaling pathway analysis method that takes into account signal variations that effector genes receive under disease and normal conditions. Our results showed that SPFA was comparable to seven other signaling pathway analysis methods. We also found that the signal variations that effector genes receive can reflect the contribution of different 0 10 20 5 9 5 8 9 6 5 6 0 1 4 6 0 9 5 6 0 2 5 5 9 9 5 5 7 8 8 0 6 1 5 5 7 9 4 7 7 3 8 1 8 8 9 4 8 1 6 8 1 7 4 3 1 6 4 7 7 5 3 7 2 5 5 4 6 7 4 7 7 6 5 5 8 2 8 1 5 4 7 7 2 9 4 7 5 6 4 7 7 Effector genes Signal variations received by the effector genes Figure 9 The signal variations received by effector genes from the upstream genes in the Wnt signaling pathway using Alzheimer's disease datasets (GSE16759).
Full-size  DOI: 10.7717/peerj.9695/ fig-9 functional attributes in the signaling pathway, deepening our understanding of disease cells' biological behaviors. Additionally, SPFA used the effector genes with high signal variations to find the abnormal part of the disease-related pathway. However, SPFA was weaker than MRGSE, Fisher, and PADOG when comparing the median p-values of target pathways. We assume this is due to the statistical models used. The probability P sd is evaluated by gene IDs permutation. Correlation differences are sometimes used to establish differential co-expression networks. This indicates that high correlation differences may exist in randomly-selected paired genes. The p-values may  Figure 10 The distribution of DEGs in Wnt signaling pathway using colorectal cancer datasets (GSE4183). The nodes with grey color mean that these nodes contain DEGs; the nodes with white color mean that these nodes do not contain DEGs. Full-size  DOI: 10.7717/peerj.9695/ fig-10 increase when paired genes with high correlation differences are randomly selected. Future studies should use a better statistical model to resolve this problem. Additionally, the 33 gene expression datasets used in this work were still limited. More experiments need to be conducted to further validate SPFA's performance. A large number of normal and disease samples are also needed to locate the effector genes with high signal variations in disease-related pathways. These genes could then serve as effective module biomarkers for accurately detecting or diagnosing complex diseases, or as drug discovery targets. SPFA depends on manually-curated signaling pathways which play a small role in complex cellular progression. More signaling pathways need to be discovered for SPFA's optimal performance.

CONCLUSIONS
In this study, we developed a new signaling pathway analysis method called SPFA. We compared this method's ability to identify altered signaling pathways against the other seven methods. SPFA showed better results than the seven other methods. Our results also showed that the SPFA method could help identify abnormal functional attributes under normal and disease conditions and the abnormal parts of a pathway during the disease biological process.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by grants from the National Natural Science Foundation of China (61871121, 61271055, and 61471112). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: National Natural Science Foundation of China: 61871121, 61271055, and 61471112.