The preprocess module can preprocess the given GEO data (GSEXXX_RAW.tar and GSEXXX_series_matrix.txt.gz) stored in the assigned working directory (e.g., PathwayKO_platform) by utilizing the oligo (21) and limma (22, 23) packages for RMA normalization, which will generate intermediate data (pData.csv and PREP.RData). The preprocess module can handle most types of GEO data produced by the major types of machines thus far. Meanwhile, the makeSPIAdata module adapted from the SPIA package (16) can parse KEGG pathways to generate intermediate data (mmuSPIA.RData) ready for pathway enrichment analysis by topology-based methods (SPIA, ROntoTools_PE, and ROntoTools_pDIS). Only topology-based methods can utilize the topology information of gene-gene interactions (5, 16). All resulting output files are stored in the new directories automatically created and named after the given data ( Supplemental Figure S2 ).

The main module pathwayko can automatically select DEGs through enabling the HES (high-edge-score) approach (17), in addition to classical approaches (5, 7, 8). The main module pathwayko drives the external changepoint package (18), by which the change-point analysis method (see Figure 1 ) makes a statistical decision on choosing the differentially expressed genes (DEGs) based on the distribution of edge scores ( Figure 2 ) that are constructed from data (17). The resulting output directory will be automatically created and named after the given data; the directory contains a set of output files including the list of DEGs, the list of knockout (KO) KEGG pathways, and the KO gene-associated subnetwork.

The principles behind the high-edge-score (HES) approach and the change-point analysis method for statistically selecting differentially expressed genes (DEGs) were modified from the literature (6, 17, 18, 23, 24), as briefly depicted below ( Figure 2 ). First, a global graph is constructed ( Figure 2A ) by the external KEGGgraph package (24) based on available KEGG pathways (6); this graph comprises all known interactions among the entire gene sets present in the KEGG pathways. Second, the edge scores for all edges in the global graph are calculated ( Figure 2B ). Given two genes (X and Y) connected by an edge in the global graph, suppose FC_X and FC_Y are the values of the expression fold-change (FC) of X and Y, respectively, while pX and pY are the probability of observing FC_X and FC_Y just by chance, then the edge score between X and Y is calculated by the following formula (1):

where the FC value for gene X (or Y), FC_X (or FC_Y ), is calculated by comparing the expression values of KO samples versus normal samples in each data. The p-value for gene X (or Y), pX (or pY), is calculated between the same groups (i.e., KO samples vs. normal samples) by using a moderated t-test (17). The FC values and p-values for all of these genes are calculated by using the functions eBayes and topTable from the external limma package (23). Third, a change-point is statistically determined ( Figure 2C ) by using the change-point analysis method from the external changepoint package (18), where a change-point is defined to be an inflection point after which the distribution curve becomes flat (18). Fourth, important edges are statistically determined ( Figure 2D ). They are connected by DEGs that are determined by a chosen HES threshold once initialized. Such DEGs connecting important edges are statistically selected, and ready for subsequent pathway enrichment analysis. Fifth, the distribution of edge scores is generated with three optional HES thresholds defined ( Figure 2E ). HES1 is the least HES (the lower boundary of the change point), the beginning of the flat area of a curve, recommended as default; HES2 is the extra 25% HES (25% safety margin of the change point), top 75% of the remaining scores to avoid selecting an overwhelming number of DEGs that may introduce false positives; and HES3 is the maximal HES (the upper boundary of the change point), the last edge connected by last two genes with the highest score. Finally, the differentially expressed genes (DEGs) are statistically selected by HES1, HES2 or HES3.

The main module pathwayko can conduct the desired computations in a pipeline fashion, i.e., across methods over one data (see Figure 1 ). These computations include (i) building an ROC curve, (ii) computing the key metrics ( Table 1 ), (iii) computing the AUC (with 95% CI, confidence interval) for the full area under the entire ROC curve, (iv) computing the partial AUCs (pAUC_SP and pAUC_SE) for the specific regions focusing on 90–100% specificity and sensitivity in both original and corrected formats (25), and (v) computing the 95% CIs for specificity and sensitivity. The entire set of key metrics ( Table 1 ) for one data are then formatted as a summary output file (SUM.RData). The resulting output files are stored in a new directory automatically generated and named after the data. Such batch-computations intensively consume computing resources (CPUs, memory and storage). This is a time-limiting step.

The principles behind the main module pathwayko (see Figure 1 ) are modified from the literature (5, 19, 26), as briefly depicted below. (i) By our definition, a true positive KO (TPKO) signaling pathway is a pathway that contains the knockout (KO) gene and was correctly identified to be a significantly impacted by that KO gene (e.g., at the pathway-level p-value < 0.001). (ii) A false positive KO (FPKO) signaling pathway is a pathway that does not contain the KO gene, and was not significantly impacted by that KO gene, but was still identified to be significantly impacted. A true negative KO (TNKO) signaling pathway is a pathway that does not contain the KO gene and was not significantly impacted by that KO gene; thus, it was not reported to be significantly impacted. A false negative KO (FNKO) signaling pathway is a pathway that contains the KO gene and was significantly impacted by that KO gene, but was not reported to be significantly impacted. (iii) For such a true-false case, the response versus prediction with a probability allows us to employ the external pROC package (19) to compute a set of key metrics defined by the above terms ( Table 1 ). (iv) An ROC curve represents the tradeoff between specificity and sensitivity for every possible p-value (19). The Youden’s best p-value threshold (denoted as p-Threshold) is a p-value that defines an optimal point (specificity, sensitivity) on an ROC curve (26), where the sum of specificity and sensitivity is maximal (19, 25). Each point on an ROC curve represents a true KO (both true positive and false negative) signaling pathway in our cases. (v) Some key metrics with local properties (FDR, FPR, FNR, specificity, sensitivity, accuracy, precision and recall) collected at the p-threshold ( Table 1 ) can be used to conduct a local comparison; others with global properties (AUC, pAUC_SP and pAUC_SE) can be applied to perform a global comparison. And these key metrics are appropriate for evaluating the performance difference of pathway analysis methods and for assessing the quality of data, both in terms of the ROC curve-based statistics analysis.

The internal evidenceplot module adapted from the SPIA package (16) can construct the two-dimensional probability evidence plots. The Adj.p-value is used to control the false discovery rate (FDR) for multiple testing (27, 28). Pathways above the oblique blue (or red) line are significant at 5% after BH-FDR (or Bonferroni) correction of the global p-values (pG), i.e., an adjusted Fisher’s product of pPERT and pNDE (16). Such plots can differentiate each data under comparison after having been analyzed by the SPIA method in view of the resolution along the X- and Y-axes (16). The highest resolution along the Y-axis, as indicated by the fine distribution, suggests that the most likely target pathways are identified, thus coinciding with its superiority among top-ranked TPKO signaling pathways (e.g., a list of top-30 ranked).

The internal roctest module adapted from the pROC package (19) can conduct the ROC curve-based statistical hypothesis testing both on the two ROC curves themselves (via the venkatraman method) and on the values of AUC, pAUC_SP and pAUC_SE (via the bootstrap method). The output files are stored in a new directory (roctest) automatically generated and named after the module.

The internal combineresult module can combine individual results of each data across methods, and format them to be a summary output file (STATS.RData). Thereby, the internal violinplot module and the internal wilcoxtest module can conduct the ROC curve-based statistics analysis, comparing the individual key metrics ( Table 1 ) among the methods under comparison. All resulting output files are stored in the respective directories (combineresult, violinplot, and wilcoxtest), which are automatically created and named after the modules.

The internal filtertrue module can extract the signaling pathways, and the internal pathwayview module adapted from the Pathview package (20) can render them. The resulting output files are stored in the two directories automatically generated and named after the two modules (filtertrue and pathwayview). Key regulators, which are impacted by the same KO gene and are coordinately, significantly up- or down-regulated, can be differentially marked on the target signaling pathways. These rendered signaling pathways illuminate putative mechanisms underlying the KO phenotype at the systems-level.

The work-flow logic of utilizing individual modules of the PathwayKO platform is indicated below ( Figure 3 ). The modules can be executed in interactive mode or in pipelined batch-execution mode.

After a number of data having been preprocessed in a step-by-step manner and stored in the user’s working directory ( Supplemental Figure S2 ), the user may conduct a batch-execution by typing the pathwayko_batch () command within the continued R session ( Supplemental Figure S3 ). This command will start sequential operations: (i) to scan and identify data generated previously by the preprocess module; and (ii) to apply the main pathwayko module to conduct a batch-execution, performing the integrative KO pathway enrichment analysis in a pipeline fashion, both across methods and across data ( Supplemental Figure S3 ).

Once a collection of benchmark data have been previously preprocessed and stored in a user’s working directory ( Supplemental Figure S2 ), the user should start a new R session to complete the desired batch-computations by following the offered tutorials, where some key parameters can be initialized in a step-by-step manner ( Supplemental Figure S3 ). All resulting output files should be finally stored in the new directories that are automatically created and named after the modules ( Supplemental Figure S4 ).

The original bench-experiments suggested that (i) wild-type (WT) mice had a high mortality rate after sterile 85% liver resection; (ii) single null Myd88 (Myd88 ^–/–) mice had greatly reduced survival rates; (iii) double null Myd88_ Ager (Myd88 ^–/– /Ager ^–/–) mice had even more greatly reduced survival rates; but (iv) single null Ager (Ager ^–/–) mice had drastically increased survival rates, which suggested opposing roles between MyD88 (encoded by Myd88) KO and RAGE (encoded by Ager) KO (3). The GSE22873 GEO data were created from tissues recovered from the remnants of mouse livers 2 hours after the sterile 85% liver resection in the indicated deficient and wild-type mice (3). These data were not subjected to bioinformatics analysis, similar to our perspectives, although some DEGs and pathways were suggested by primary analyses (3) with Pathway Express (16). To investigate the mechanisms underlying such phenomena raised by the original bench-experiments (3) from the perspective of systems immunology, we applied the PathwayKO platform to analyze the following six subtypes of data we assigned: GSE22873_M (KO Myd88 vs. WT), GSE22873_MA (KO Myd88_Ager vs. WT), GSE22873_A (KO Ager vs. WT), GSE22873_MAvM (KO Myd88_Ager vs. KO Myd88), GSE22873_MAvA (KO Myd88_Ager vs. KO Ager) and GSE22873_MvA (KO Myd88 vs. KO Ager). The results are presented as follows.

The distribution of edge scores varied drastically among data ( Figure 4 ). Each showed a distinction, as marked by an HES threshold (HES1, HES2 and HES3, respectively), and reflected the different effects of experimental treatments.

DEGs were statistically selected by using an HES threshold, as described in Methods (see Figures 2 , 4 ). By the HES1, HES2 and HES3 thresholds, respectively, GSE22873_M separately had 122 DEGs, 38 DEGs and 2 DEGs ( Supplemental Table S1 ) with the top 10-ranked genes (Myd88, Cxcl1, Fos, Nfatc3, Ccr1l1, Cxcr2, Ccr4, Rela, Jun and Xcr1); GSE22873_MA had 130 DEGs, 57 DEGs and 2 DEGs ( Supplemental Table S2 ) with the top 10-ranked genes (Myd88, Cxcl1, H2-Bl, Ccr1l1, Ccr10, Cxcr2, Rela, Xcr1, Ccr1 and Cxcr5); and GSE22873_A had 135 DEGs, 27 DEGs and 2 DEGs ( Supplemental Table S3 ) with the top 10-ranked genes (H2-Bl, H2-T24, H2-M3, H2-Q6, Cdkn1a, Tapbp, AP1g1, Ap1b1, Ap1m1 and H2-M10.1). Similarly, by the HES1 threshold, there were 55 DEGs (with top-ranked H2-Bl and H2-T24) in GSE22873_MAvM ( Supplemental Table S4 ); 67 DEGs (with top-ranked Myd88 and Cxcl1) in GSE22873_MAvA ( Supplemental Table S5 ); and 347 GEGs (with top-ranked H2-Bl, H2-M10.1, Myd88 and Cxcl1) in GSE22873_MvA ( Supplemental Table S6 ). Of note, rather than Ager itself, the MHC/antigen-related genes were ranked in the top 10 in GSE22873_A, which compared RAGE-null with wild-type mice. None of the above lists contained Ager, suggesting that (i) Ager itself was not yet significantly differentially expressed as early as 2 hours post-surgery when comparing RAGE-null against wild-type mice in GSE22873_A (KO Ager vs. WT), born by the transcriptome, in our current analysis; and (ii) RAGE (Ager)-mediated MHC/antigen-related genes play the most important roles, implying epigenetic effects.

The KO gene-associated subnetwork was comprised of the statistically selected DEGs ( Figure 5 ). The landscapes of these subnetworks explicitly revealed the global effects of experimental treatments at the systems-level. For instance, the top-ranked two genes (Myd88 and Cxcl1) connected the last edge, as highlighted by the HES3 threshold (see Figure 4 ); this result suggests that the KO Myd88-signaling stimulated the highest responsive expression of Cxcl1 in both GSE22873_M ( Supplemental Table S1 ) and GSE22873_MA ( Supplemental Table S2 ). Similarly, the top-ranked two genes (H2-Bl and H2-T24) connected the last edge (see Figure 4 ), suggesting that the KO Ager-signaling stimulated the highest responsive expression of H2-Bl and H2-T24 in GSE22873_A ( Supplemental Table S3 ). In addition, the remaining pairs displayed drastic variations, reflecting subtractive responses between the indicated pairs ( Figure 5 ).

The ROC curves with statistical features ( Figure 6 ) were generated in a pipeline fashion, both across methods and across data. These results indicate the diverse landscapes of experimental treatments (data), which can discriminate the performance difference of pathway analysis methods under comparison. (i) An AUC with 95% CI and Youden’s best p-value threshold (noted shortly as p-Threshold) were marked on each ROC curve, indicating an optimal point (specificity, sensitivity) for a method. (ii) The 95% CIs of both specificity and sensitivity were marked on each ROC curve, indicating the deviation for a method on one data point. For instance, the drastic variation in 95% CIs shown in GSE22873_A (KO Ager vs. WT) data suggested weak signals over noise, which might impact the statistical hypothesis testing (19). (iii) The partial AUCs (pAUC_SP and pAUC_SE) of 90–100% specificity and sensitivity were marked on each ROC curve, implying a bias toward a higher specificity and sensitivity, respectively. (iv) The multiple ROC curves with AUCs intuitively displayed the superiority of methods when benchmarked on the same data ( Figure 6 ).

The ROC curve-based statistical hypothesis testing was conducted to indicate the significance level of the difference between the two methods ( Figure 7 ). For instance, two methods (SPIA vs. GSEA) had a significant difference (e.g., p-value < 0.05) when benchmarked on GSE22873_M and GSE22873_MA, respectively, in view of both ROC curves themselves (via the venkatraman method) and the values of AUC, pAUC_SP and pAUC_SE (via the bootstrap method). SPIA was better than GSEA in these two cases ( Figures 7A, B ). However, the same two methods (SPIA vs. GSEA) showed no significant difference (e.g., p-value > 0.1) when benchmarked on GSE22873_A ( Figure 7C ). The reason is likely attributed to the drastic variations in the 95% CIs of specificity and sensitivity in GSE22873_A (see Figure 6C ), which impact the resampling-based bootstrap assessments used for the ROC curves-based statistical hypothesis testing (19).

The two-dimensional probability evidence plots were created for six data ( Figure 8 ), which discriminate these data in view of the resolution along the X- and Y-axes. The highest resolution along the Y-axis suggests that the most likely target pathways were identified by SPIA, as indicated by the fine distribution.

We identified diverse true KO signaling pathways, which contained target KO gene(s), including 21 from GSE22873_M (KO Myd88 vs. WT), 23 from GSE22873_MA (KO Myd88_ Ager vs. WT), 3 from GSE22873_A (KO Ager vs. WT) and 24 from GSE22873_MvA (KO Myd88 vs. KO Ager) ( Table 2 ; Supplementary Table S7 ). The majority of them were identified to be TPKO signaling pathways at a statistical significance level (e.g., pGFdr < 0.001 or 0.005) by our definitions described in Methods. (i) All 21 pathways from GSE22873_M were TPKO signaling pathways. (ii) Twenty-two out of 23 pathways were TPKO signaling pathways, whereas only 1 (mmu05010: Alzheimer disease) out of 23 pathways from GSE22873_MA was a false negative knockout (FNKO) signaling pathway; i.e., this pathway contained a target KO gene and should have been significantly impacted, but was not identified to be positive at a significance level (pGFdr > 0.005). (iii) Only 1 (mmu04933: AGE-RAGE signaling pathway in diabetic complications) out of 3 pathways from GSE22873_A was a TPKO signaling pathway, whereas the remaining two pathways were FNKO signaling pathways (pGFdr > 0.005). (iv) Twenty-three out of 24 pathways from GSE22873_MvA were TPKO signaling pathways, whereas only 1 (mmu05150: Staphylococcus aureus infection) out of 24 pathways was an FNKO signaling pathway (pGFdr > 0.005) ( Supplementary Table S7 ).

We categorized target TPKO signaling pathways into the following two groups (1): the basal signaling pathway that is solely responsible for the KO-gene phenotype; and the composite signaling pathway that embeds the named basal signaling pathway. Here, the basal signaling pathway for KO Myd88 was the Toll-like receptor signaling pathway (mmu04620), while the composite signaling pathways were the remaining Myd88-containing signaling pathways ( Table 2 ; Supplementary Table S7 ). Similarly, the basal signaling pathway for KO Ager was the AGE-RAGE signaling pathway in diabetic complications (mmu04933), while the composite signaling pathways were the remaining Ager-containing pathways ( Table 2 ; Supplementary Table S7 ). Diverse target TPKO signaling pathways were mainly linked to innate inflammatory responses.

The target TPKO signaling pathways identified by the SPIA method (see Table 2 ; Supplementary Table S7 ) were rendered with key regulatory proteins (encoded by genes). Key regulators included ligands, receptors, adapters, transducers, transcription factors and cytokines, which were coordinately, significantly and differentially expressed at the systems-level. Such key regulators may constitute a signaling-axis, as marked on a target TPKO signaling pathway. Some top-ranked target TPKO signaling pathways are highlighted below as examples ( Figures 9 – 11 ).

Based on the above comprehensive characterizations of individual data in this study and in our previous publication (1), we conclude that GSE22873_M, GSE22873_MA, GSE24327_A, GSE24327_B and GSE24327_C have strong signal to noise ratios when generated from the original bench-experiments (2, 3). These datasets are qualified to be included in a collection of qualified data that serves as a proper benchmark used to exemplify the different tasks in the present study ( Supplemental Figure S2 ).

The operational scripts depict batch-executions, exemplifying how target KO signaling pathways are automatically deciphered at the systems-level in a pipeline fashion, both across methods and across data ( Supplemental Figure S3 ). It takes approximately 25 minutes to complete the batch-computations across seven methods over one data on our high-performance workstation with the assigned configuration (see Figure 1 ; Supplemental Figure S1 ). This is a time-limiting step.

The ROC curve-based statistical analyses (e.g., violin plot and Wilcox test) are suitable for assessing the performance difference of the pathway analysis methods under comparison ( Supplemental Figure S4 ). The violin plot graph suggested that the topology-based methods (SPIA, ROntoTools_PE and ROntoTools_pDIS) generally performed better than the non-topology-based methods (GSEA, GSA, SAFF and PADOG); and that SPIA was overall better than ROntoTools_PE and ROntoTools_pDIS ( Figure 12 ). These results coincide well with the findings of each individual data (see Figures 6 , 7 , plus Figure 2 therein Ref (1). However, the results of the Wilcox test (data not shown) on the global metrics (AUC, pAUC_SP and pAUC_SE) suggested that no significant difference was observed among methods when benchmarked on this collection of data.

Collected at the optimal p-Threshold, the key metrics (FDR, FPR, FNR, specificity, sensitivity, accuracy, precision and recall) reflect the local properties of the ROC curves and are suitable for a local comparison among methods (see Table 1 ). With the chosen benchmark, the topology-based methods (SPIA, ROntoTools_PE and ROntoTools_pDIS) performed better than the non-topology-based methods (GSEA, GSA, SAFF and PADOG) in terms of p-Threshold and FDR; and GSEA displayed the poorest performance among the seven tested methods in terms of specificity, sensitivity, accuracy and precision ( Figure 12 ).

The key metrics including AUC, pAUC_SP and pAUC_SE reflect the overall properties of the ROC curves and are appropriate for a global comparison among methods (see Table 1 ). The partial AUCs (pAUC_SP and pAUC_SE) of the 90–100% specificity and sensitivity may display drastic variations on a case-by-case basis, thus measuring possible bias toward a higher specificity (true negative rate) or higher sensitivity (true positive rate). With the chosen benchmark, ROntoTools_PE performed better than ROntoTools_pDIS in terms of AUC, pAUC_SP and pAUC_SE; and SPIA had a higher specificity (pAUC_SP) but a lower sensitivity (pAUC_SE) than GSEA ( Figure 12 ).