Modeling causal signal propagation in multi-omic factor space with COSMOS

Article activity feed

1. 

   Review Commons

   Aug 17, 2026

   Note: This response was posted by the corresponding author to Review Commons. The content has not been altered except for formatting.

   Learn more at Review Commons

   ---

   Reply to the reviewers

   Reviewer #1: Evidence, Reproducibility And Clarity

   General comment

   Reviewer comment: Dugourd et al. present COSMOS+, a computational framework that integrates multi-omics data with prior knowledge networks to generate mechanistic hypotheses connecting signaling, transcriptional regulation, and metabolism. The study addresses a relevant challenge in the field: the difficulty of moving beyond purely data-driven factor analysis toward biologically interpretable and causally grounded insights, introducing MOON, a scalable iterative network scoring algorithm as a practical alternative to computationally expensive optimization-based approaches such as CARNIVAL.

   Note: This response was posted by the corresponding author to Review Commons. The content has not been altered except for formatting.

   Learn more at Review Commons

   ---

   Reply to the reviewers

   Reviewer #1: Evidence, Reproducibility And Clarity

   General comment

   Reviewer comment: Dugourd et al. present COSMOS+, a computational framework that integrates multi-omics data with prior knowledge networks to generate mechanistic hypotheses connecting signaling, transcriptional regulation, and metabolism. The study addresses a relevant challenge in the field: the difficulty of moving beyond purely data-driven factor analysis toward biologically interpretable and causally grounded insights, introducing MOON, a scalable iterative network scoring algorithm as a practical alternative to computationally expensive optimization-based approaches such as CARNIVAL.

   Importantly, the method has the potential to highlight errors in prior knowledge networks. This will not palliate the incompleteness of the existing prior knowledge but, at least, it can identify inconsistencies and help 'correct' the databases. It is not entirely clear whether inconsistencies stem from specificities of cell lines/samples that are not reflected in the general databases, or simply mistakes in the databases.

   In any case, reconstructing mechanistic hypotheses based on (multi)omics datasets remains an important endeavour, necessary to understand biological processes and also with clear applications in medicine.

   The code is available and clearly documented ensuring reproducibility. The literature survey is comprehensive and useful to understand the need for developments in this field.

   While the framework is flexible and the applications span a diverse range of contexts from cell line collections to clinical cohorts, some aspects of the work require stronger justification and additional evaluation tests might increase the usefulness of the methods for the community.

   Response:

   We thank the reviewer for the careful and constructive assessment of COSMOS+, including the observation that the framework may expose limitations in current prior-knowledge networks. We agree that the original manuscript needed clearer justification of several methodological choices and a more explicit account of the scope and limitations of its conclusions.

   In response, we added targeted analyses and clarifications on score interpretation, parameter robustness, multiple testing, benchmark coverage and reasons for failure, and prior-knowledge-network context dependence. These revisions frame COSMOS+/MOON outputs as standardized network-consistency scores and testable hypotheses, rather than layer-invariant hypothesis-test statistics, definitive causal mechanisms, or validated biomarkers. The detailed analyses and manuscript changes are described under R1.1–R1.10 below.

   R1.1: MOON score distributions across layers

   Reviewer comment: MOON scores are ULM t-values at layer 1, but at deeper layers they are t-values computed from previous t-values. It is unclear whether the propagated scores follow the same distribution across layers, which is a requirement to apply the same threshold of |score| > 1.5 uniformly. The authors should either provide a justification for why the threshold remains meaningful at deeper layers, or characterize through simulation how score distributions change across layers and propose a layer-specific thresholding strategy.

   Response:

   We thank the reviewer for raising this point. We agree that propagated MOON scores should not be described as formal t-distributed p-values with identical null calibration at every propagation level. The initial activity layer is inferred from molecular signatures, whereas upstream MOON scores are obtained by applying signed linear models to already-inferred downstream scores. We therefore interpret MOON scores as standardized network-consistency scores used for mechanistic prioritization, rather than as exact layer-invariant hypothesis-test statistics.

   To assess whether propagation level nevertheless introduced systematic score inflation in the CytoSig benchmark, we added an empirical score-distribution check using the estimated CytoSig MOON score. For each perturbation where the applied ligand was scorable, we extracted the matching ligand's MOON score and propagation level. These 549 matching ligand-experiment scores covered 63 ligands and levels 1–4. The applied-ligand scores were not significantly associated with propagation level (Spearman rho = 0.052, p = 0.223; linear-regression slope = −0.134 score units per level, p = 0.353, R2 = 0.0016). Thus, the benchmark did not show evidence that matching ligand scores were artificially larger simply because the ligand was farther upstream in the network.

   We also summarized all MOON scores across the 1,359 CytoSig experiments by propagation level. Score medians remained close to zero across levels, and deeper propagated layers did not show systematically broader score distributions than the direct-input or shallow propagated layers. For example, the median absolute score across experiments was 0.956 at level 0, 0.722 at level 1, 0.676 at level 2, and ranged from 0.567 to 0.686 across levels 3–10. This supports the use of a fixed absolute-score threshold as a pragmatic prioritization boundary in this benchmark, while avoiding the stronger claim that the same numerical threshold represents an identical calibrated p-value at every propagation level.

   Because the CytoSig analysis represents the effect of propagating signals across a prior knowledge network across a wide variety of distinct conditions, with almost every node of the network being scored, we expect the MOON score behavior to be generalised to other contexts as well. We have placed the detailed methods, results, and interpretation described in this response in Supplementary Text S2: “Empirical assessment of MOON score distributions across propagation levels.”

   We addressed the separate questions of maximum reachability steps and benchmark coverage under R1.7 and R1.8, respectively. The distinct use of score thresholds as candidate-input boundaries is addressed under R1.5.

   We have modified the following sections of the manuscript accordingly:

   CytoSig results, Section 2.2:

   Original manuscript excerpt:* “...COSMOS prior knowledge within the specified number of steps, therefore their score couldn’t be estimated).”*

   Updated manuscript excerpt:* “...COSMOS prior knowledge within the specified number of steps, therefore their score couldn’t be estimated). The MOON scores did not systematically increase with propagation level (see Supplementary Text S2), indicating that the same absolute score threshold can be used across levels in this benchmark.”*

   R1.2: Asymmetric consistency pruning

   Reviewer comment: The transcriptional consistency check removes TF -> target edges where TF score is incoherent with target expression, but no analogous check is applied to kinase -> TF or receptor -> kinase edges. The upstream signaling layer is therefore only constrained by the upstream anchor comparison while the transcriptional layer (TF -> target gene) undergoes aggressive pruning. The authors should justify this asymmetry explicitly or extend the consistency check to upstream layers where phosphoproteomic data is available.

   We are aware that prior knowledge on signs of activation might also be lacking, often it will be hard to know what is the real ground truth without doing specific experiments in the right cell types/conditions. Some more discussion about this would clarify the difficulty of this problem.

   Response:

   We thank the reviewer for raising this point. The TF-target coherence step is an edge-level consistency check that is used only when the directly downstream molecular readout is available. For a signed TF-target edge, expression of the target gene is measured; for example, a positive interaction between an active TF and a down-regulated target gene is incoherent with the data. This does not establish that the prior-knowledge edge is universally incorrect, but it supports removing the edge from the context-specific network for that analysis.

   By contrast, the effect of a kinase-TF interaction can depend on a particular phosphosite that is not measured, while a measured phosphosite may not be the site through which the kinase regulates the TF. An inferred kinase activity or a discordant phosphosite measurement therefore does not generally provide an equivalent direct test of the individual edge. Such discordance can reflect incomplete site coverage or context-specific regulation rather than an incorrect prior-knowledge interaction or sign, so we do not apply the TF-target pruning rule broadly to upstream signaling layers.

   Following this comment, we added two clarifications: one sentence in the Results section that directs readers to the Methods rationale, and a short Methods explanation with the kinase-TF example above. No new analysis or supplementary text was required for this methodological clarification.

   We have modified the following sections of the manuscript accordingly:

   Location: Results, Section 2.1

   Original manuscript excerpt:* “...(incoherence between sign of the TF activity score and the sign of the downstream measurement/factor weight input).”*

   Updated manuscript excerpt:* “...(incoherence between sign of the TF activity score and the sign of the downstream measurement/factor weight input). This check is specific to TF-target interactions because RNA data directly measure the expression of the target gene, whereas the available upstream measurements do not generally provide an equivalent direct readout of individual interactions (see Methods).”*

   Location: Methods, Section 4.4.1

   Original manuscript excerpt:* “If RNA data points (here, MOFA weights) are provided, a check can be performed to remove any interaction from the network that connects a TF and a downstream gene that has an incoherent expression sign with the TF MOON score. Thus, the moon function and the TF-target coherence check can be run in a loop until the output of the moon doesn’t contain any incoherence between TF scores and downstream targets. The algorithm to remove incoherent TF-target interactions is as follow:”*

   Updated manuscript excerpt:* “If RNA data points (here, MOFA weights) are provided, a check can be performed to remove any interaction from the network that connects a TF and a downstream gene that has an incoherent expression sign with the TF MOON score. Such a check can be performed for TF-target interactions because the expression of the direct target gene is measured: for example, a positive interaction between an active TF and a down-regulated target gene is incoherent with the data. In contrast, the effect of a kinase-TF interaction can depend on phosphorylation at a specific site that is not measured, and an inferred kinase activity or a measurement at a different phosphosite may therefore not provide an equivalent direct readout of that interaction. Thus, the moon function and the TF-target coherence check can be run in a loop until the output of the moon doesn’t contain any incoherence between TF scores and downstream targets. The algorithm to remove incoherent TF-target interactions is as follow:”*

   R1.3: MOFA scale_views parameter

   Reviewer comment: For the MOFA section, the authors state they used the parameter viewsscale_views = False. This means the three omics layers were not variance-normalized before input to MOFA. Since transcriptomics, proteomics, and metabolomics have different features ranges, the omic layer with higher absolute variance will dominate the factor structure. The authors should either justify this choice explicitly or show that the factor structure is not dominated by a single omic layer.

   Response:

   We thank the reviewer for highlighting this issue. To assess whether the NCI60 MOFA factor structure depended on scale_views = False, we trained a sensitivity model using the same configuration as the selected analysis (a maximum of 10 factors, yielding 9 active factors), changing only scale_views to True.

   The view-scaled model also retained 9 active factors. The total variance explained by RNA, metabolomics, and proteomics was essentially unchanged (59.5999%, 18.7037%, and 23.2947% in the original model versus 59.5989%, 18.6754%, and 23.2887%, respectively, in the view-scaled model). Factor scores showed a one-to-one correspondence between the two models, with absolute Pearson correlations of at least 0.9993 across all nine factors; opposite correlation signs reflect the arbitrary orientation of latent factors.

   These results show that the selected factor structure is not materially affected by view-scale normalization. We therefore retain the original analysis and interpretation reported in the manuscript.

   R1.4: Multiple-testing correction

   Reviewer comment: There is no multiple testing correction in the clinical association analysis. The ULM-based clinical metadata association tests each factor against each clinical category independently. With 9 factors and the number of clinical categories available in NCI60 (tissue of origin alone has ~10 categories, plus age, pathology, and other variables), the number of simultaneous tests is large. The same issue applies to the Cox survival analysis in the Paloma3 section, where a separate Cox model is fitted for every node in the MOON network and results are reported without any correction for multiple comparisons. The authors should apply FDR correction at both stages to confirm that reported associations are not false positives arising from the large number of simultaneous tests performed.

   Response:

   We thank the reviewer for raising this important point. We have now applied Benjamini-Hochberg correction to the NCI60 metadata-factor and PALOMA3 Cox test families.

   For the NCI60 analysis, we corrected 297 metadata category-factor tests (33 categories x 9 factors). The principal associations used to interpret the factors remained significant: Factor 2 was associated with melanoma origin (nominal p-value = 3.67 x 10^-15, FDR = 5.45 x 10^-13), Factor 4 was negatively associated with leukemia origin (nominal p-value = 2.19 x 10^-9, FDR = 2.17 x 10^-7), and Factor 2 was negatively associated with epithelial origin (nominal p-value = 1.69 x 10^-6, FDR = 1.00 x 10^-4).

   For PALOMA3, we distinguished the two-arm treatment-by-biomarker interaction analyses from the separately evaluated treatment-arm node-wise Cox scans. The latter did not yield an individual MOON score with FDR

   In response to this comment, we now report FDR values beside the principal NCI60 nominal p-values, revise the PALOMA3 wording to report the RB1 interaction FDR and separate it from the node-wise scan, and describe the correction families and verified ULM/MOON settings in the Methods.

   We have modified the following sections of the manuscript accordingly:

   Location: Results, Section 2.3

   Original manuscript excerpt:* “Factor 4 was showing a significant negative association with samples of Leukemic origin. We also saw that factor 2 was significantly associated with a Melanoma origin, and negatively associated with an Epithelial origin (Figure 3C).”*

   Updated manuscript excerpt:* “After Benjamini-Hochberg correction across the 297 metadata category-factor tests (33 metadata categories x 9 factors), Factor 4 was negatively associated with samples of leukemic origin (nominal p-value = 2.19 x 10^-9, FDR = 2.17 x 10^-7). Factor 2 was associated with melanoma origin (nominal p-value = 3.67 x 10^-15, FDR = 5.45 x 10^-13), and negatively associated with epithelial origin (nominal p-value = 1.69 x 10^-6, FDR = 1.00 x 10^-4; Figure 3C).”*

   Location: Results, Section 2.6

   Original manuscript excerpt:* “CDK2 and RB1 MOON scores were both found to be significantly associated with worse response in the treatment arm, but not at their expression level. Furthermore, RB1 coefficient relative direction is reversed between its expression and MOON score. A lower MOON score of RB1 in patients is significantly associated with worse patient response in the treatment arm, while a high expression was marginally associated with worse patient response (MOON interaction p-value = 0.008, RNA interaction p-value = 0.12). Since RB1 is a known inhibited target of CDK4 and 6 as well as a tumor suppressor (Knudsen et al, 2019), its MOON score direction is more consistent with expectation than its expression.”*

   Updated manuscript excerpt:* “CDK2 and RB1 MOON scores showed nominal treatment-by-biomarker interaction associations with worse response in the treatment arm. Furthermore, the relative direction of the RB1 coefficient was reversed between its expression and MOON score. A lower MOON score of RB1 in patients was nominally associated with worse patient response in the treatment arm (MOON interaction p-value = 0.008, FDR = 0.743), whereas high expression was not (RNA interaction p-value = 0.12). In the separately evaluated treatment-arm node-wise Cox analyses, no individual MOON score remained significant after Benjamini-Hochberg correction (FDR *

   Location: Methods, Section 4.11

   Original manuscript excerpt:* “Paloma3 cohort analysis Patient level RNA counts were z-transformed across patients. The ULM method of decoupleR was used with CollecTRI to estimate patient specific transcription factor activity signatures (XXX min target per TF). Patient specific COSMOS networks were generated using MOON with XXX steps up-stream from TFs. Using patient specific progression free survival values, COX survival models were computed for each gene of the one hand and each MOON score on the other hand, by splitting first patients by median gene expression or median MOON score, and then computing the COX hazard ratio between control and treatment group first, and second between high and low expression/moon scores.”*

   Updated manuscript excerpt:* “Paloma3 cohort analysis Patient level RNA counts were z-transformed across patients. The ULM method of decoupleR was used with CollecTRI to estimate patient specific transcription factor activity signatures, retaining TFs with at least 5 measured targets. Patient specific COSMOS networks were generated using MOON with a maximum of 10 upstream layers from the TF activities; for COX analyses, score matrices were restricted to nodes at levels 0–5. Using patient specific progression free survival values, COX survival models were computed for each gene of the one hand and each MOON score on the other hand, by splitting first patients by median gene expression or median MOON score, and then computing the COX hazard ratio between control and treatment group first, and second between high and low expression/moon scores. Nominal p-values from the interaction and treatment-arm Cox analyses were adjusted separately using the Benjamini-Hochberg method within the MOON-score and RNA-expression test families.”*

   R1.5: Uniform t-value thresholds

   Reviewer comment: The threshold of |t| > 2 is applied uniformly across TF activity scoring, kinase scoring, LR scoring, and clinical associations without justification. A t-value of 2 corresponds approximately to p

   Response:

   We agree that the role and statistical interpretation of the score thresholds required clarification. The numerical threshold is not treated as a common, calibrated p-value cutoff across the analyses. Rather, its role is analysis-specific. In network construction, score thresholds restrict the set of candidate inputs considered for mechanistic hypothesis generation. In the NCI60 clinical association analysis, |t| > 2 screen was used for descriptive heatmap prioritization, whereas statistical inference is now based on Benjamini-Hochberg correction across the 297 metadata category-factor tests (R1.4).

   To test the candidate-input boundary directly, we have now performed a focused sensitivity analysis of the NCI60 factor 4 TF-to-ligand MOON branch. We repeated this branch using absolute upstream TF score thresholds of 1.5, 2 (the submitted setting), and 2.5, while retaining the same prior-knowledge network, downstream ligand scores, and MOON settings. The re-estimated TF-to-ligand result was combined in each case with the same stored receptor-to-TF/metabolite branch.

   The threshold changed the number of TF candidates retained after prior-knowledge-network filtering from 104 at 1.5 to 84 at 2 and 64 at 2.5. Relative to the threshold-2 result, the runs with these alternative thresholds had Spearman score correlations of 0.9990 and 0.9996 and score-sign agreement of 99.84% and 99.89%, respectively; furthermore for both alternative thresholds the top 50 absolute-score nodes were preserved.

   The pathway-control analysis results were similarly stable: all 20 highest-ranked pathways were retained at the permissive threshold and 19 of 20 at the stringent threshold; 19 of the top 20 node-pathway pairs were retained at both alternatives. The principal focal-adhesion, neurotrophin, MAPK, ERBB, and cancer-related pathway-control signals were recovered in each run.

   We added a statement to the NCI60 Results and placed the detailed methods, numerical results, interpretation, and scope limitation in Supplementary Text S3, “Sensitivity of the NCI60 factor 4 MOON analysis to the upstream TF candidate threshold.” This branch-specific, reassembled sensitivity analysis supports robustness of the NCI60 factor 4 TF-to-ligand findings to this candidate-input boundary, but does not claim that a score of 2 has an identical p-value interpretation across datasets.

   We have modified the following sections of the manuscript accordingly:

   Location: Results, Section 2.3

   Original manuscript excerpt:* “To find which biological processes are captured in the moon network, we perform a pathway over-representation analysis with sets of nodes down-stream of top deregulated MOON score nodes (absolute MOON score > 1.5) Then, we can represent pathways that are significantly over-represented downstream of given top scoring nodes of the MOON thresholded network in a heatmap (Figure 3E). Reassuringly, this analysis found expected control mechanisms, such as JAK1 or IL6ST very significantly controlling the JAK-STAT signaling pathway or ITGB (integrins) family members controlling focal adhesion. It also allows us to propose chemicals that can potentially control signaling pathways such as Acetaminophen or 4-hydroxy-oestradiol controlling the MAPK pathway.”*

   Updated manuscript excerpt:* “To find which biological processes are captured in the moon network, we perform a pathway over-representation analysis with sets of nodes down-stream of top deregulated MOON score nodes (absolute MOON score > 1.5) Then, we can represent pathways that are significantly over-represented downstream of given top scoring nodes of the MOON thresholded network in a heatmap (Figure 3E). Reassuringly, this analysis found expected control mechanisms, such as JAK1 or IL6ST very significantly controlling the JAK-STAT signaling pathway or ITGB (integrins) family members controlling focal adhesion. It also allows us to propose chemicals that can potentially control signaling pathways such as Acetaminophen or 4-hydroxy-oestradiol controlling the MAPK pathway. These results appeared overall robust to changes in the initial threshold used to select candidate TFs for network construction (absolute score thresholds: 1.5, 2 [current setting], and 2.5; Supplementary Text S3).”*

   R1.6: Cytosig minimum regulon size

   Reviewer comment: The methods section does not report what minimum regulon size was used for the Cytosig analysis, and no justification is provided for this parameter. It is not clear if it is still 10 or a different one. Cytosig signatures vary widely in how many genes were measured across experiments, meaning a uniform minimum threshold will exclude TFs not because they are inactive but simply because their targets were not measured in a given experiment. The authors should explicitly report the threshold used, assess whether using a lower threshold such as 5 recovers additional scorable signatures without substantially degrading TF activity reliability, and report the distribution of gene coverage across Cytosig signatures to contextualize how many signatures are affected by this limitation.

   Response:

   We thank the reviewer for pointing out that this parameter was not reported clearly enough. The CytoSig benchmark script calls decoupleR::run_ulm() without an explicit minsize (minimal number of measured downstream targets) argument. The previously computed TF-activities independently confirms that the effective setting in the submitted benchmark was minsize = 5: for each of the 1,359 signatures, the number of cached TF activities exactly matched the number of CollecTRI TFs with at least five measured targets (724,571 TF-signature activity estimates in total). Thus, the benchmark already used an effective five-target minimum, rather than a threshold of 10.

   To further explore how much the results would change with a higher cutoff, We have quantified the effect of comparing this setting with a hypothetical minsize = 10 cutoff. The 1,359 filtered CytoSig signatures contained a median of 18,423 measured genes per signature (IQR 16,720–19,073) and a median of 6,005 measured CollecTRI target genes (IQR 5,534–6,085). With the effective five-target minimum, a median of 539 TFs per signature were retained; at 10 targets, this would decrease to 440. Across all TF-signature pairs, the five-target setting retained 724,571 activity estimates, compared with 590,207 under the ten-target cutoff. The five-to-nine-target group therefore contributed 134,364 additional TF-signature estimates (18.5% of the TF-activity input layer).

   The minsize setting filters returned TFs rather than changing the ULM fit for TFs that pass both cutoffs. Scores for TFs with at least 10 measured targets would therefore be identical under the two settings. The added five-to-nine-target estimates had lower median absolute ULM scores than the >=10-target estimates (0.69 versus 0.95) and were less often above an absolute score of 2 (12.7% versus 21.9%). While this highlights an expected association between score magnitude and number of targets, it does not represent an independent validation of the reliability of low-coverage TF activities. We therefore describe the setting as a coverage tradeoff and do not claim that the comparison establishes unchanged reliability.

   We added the effective minimum to the Cytokine scoring Methods and placed the detailed implementation check, coverage distribution, score-magnitude context, and limitation in Supplementary Text S4, “Effect of the minimum CollecTRI target-set size on CytoSig TF-activity coverage.” The supplement also clarifies that this parameter controls the availability of TF inputs; applied-ligand scoring additionally depends on ligand mapping and prior-knowledge-network reachability.

   We have modified the following sections of the manuscript accordingly:

   Location: Methods, Section 4.5.3

   Original manuscript excerpt:* “...referred to as the TF score. Then, for each resulting TF score profile, we filtered out specifically the COSMOS…”*

   Updated manuscript excerpt:* “...referred to as the TF score. Only TFs with at least 5 measured CollecTRI targets were retained for scoring (minsize = 5; see Supplementary Text S4). Then, for each resulting TF score profile, we filtered out specifically the COSMOS…”*

   R1.7: MOON reachability step limit

   Reviewer comment: The maximum number of propagation steps for the reachability filtering is not reported or justified. In the case of the Cytosig analysis, the paper states that 31 ligands could not be scored because they were not reachable upstream of TFs within the allowed number of steps, but never states what that number was. This parameter directly determines which ligands can be benchmarked, but no sensitivity analysis is provided showing whether increasing the step limit recovers additional ligands or changes the benchmark results. The tool might be enhanced with a report of alternative values, discussion of optimal values and discussion of the tradeoff between reachability and the MOON score reliability at greater network distances.

   Response:

   We thank the reviewer for highlighting that the rationale and sensitivity of this parameter were not sufficiently clear. The Methods already stated that the CytoSig reachability filter used ten steps; the CytoSig benchmark script uses the same setting for MOON propagation (n_steps = 10, passed as n_layers = 10) and for the pre-MOON reachability filter. We use this value as a permissive upper bound, rather than as a biologically privileged path length: it avoids considering arbitrarily long prior-knowledge-network paths while allowing multi-step receptor, signaling, and TF routes.

   We assessed the observed depth sensitivity in the submitted benchmark. The 549 cleaned applied-ligand scores across 63 ligands occurred only at MOON levels 1–4 (205, 254, 83, and 7 entries, respectively); no scored applied ligands occurred at levels 5–10. Restricting the benchmark to levels 1, 2, 3 retained 205, 459, and 542, entries across 8, 45 and 60 ligands, respectively. Once levels 1–4 were included, the benchmark summary was unchanged through level 10: 47 of 63 ligands had a positive mean score, while 16 and 4 ligands had nominally significant positive and negative one-sample tests, respectively.

   We also performed a reachability-only check over ten steps for the 52 evaluable (as they are mapable on our network) experiments with missing scores, spanning 11 ligands. After the same expression and prior-knowledge-network membership filtering as in the original CytoSig analysis, none had a directed path from the applied ligand to the filtered TF-input layer at any finite distance. This indicates that simply increasing the propagation limit would not recover the missing applied-ligand scores.

   We have added clarification in the methods (pasted below), and a Supplementary Text S5, “Sensitivity of the CytoSig MOON benchmark to the propagation/reachability step limit.” The supplementary text presents the implementation setting, the level-restricted analysis and its limitation, the extended-horizon reachability-only check, and the limit of what these data can establish about deeper-score reliability. We have modified the following sections of the manuscript accordingly:

   Location: Methods, Section 4.5.3

   Original manuscript excerpt:* “...within ten steps upstream of the TFs.”*

   Updated manuscript excerpt:* “...within ten steps upstream of the TFs. We used this ten-step maximum as a permissive upper bound; in the CytoSig benchmark, applied ligands that were scored by MOON occurred at levels 1–4 (Supplementary Text S5).”*

   R1.8: Benchmark coverage

   Reviewer comment: The benchmark scores only 549 out of 1359 available signatures (40% of the data) due to PKN reachability and TF coverage limitations. Since both limiting parameters are neither reported nor optimized, it is unclear whether the excluded 60% represents a fundamental limitation of the method or an artifact of conservative parameter choices. The authors should explore alternative parameter values and report their effect on both benchmark coverage and performance jointly.

   Response:

   We agree that the parameter choices underlying this benchmark required clearer justification.

   The 549/1,359 value is the coverage of the direct applied-ligand recovery benchmark, rather than the fraction of CytoSig signatures for which MOON can generate a score table. This issue is addressed by the preceding analyses (see response R1.6, R1.7). Supplementary Text S4 documents that the CytoSig TF-activity input layer already used the effective permissive setting minsize = 5, and quantifies the reduction in TF-input coverage that would result from a ten-target minimum. Supplementary Text S5 documents the shared ten-step reachability/propagation setting and its sensitivity analysis: the observed benchmark endpoint was reached by level 4 in a post-hoc level-restricted analysis of the submitted MOON scores, while extending the reachability horizon did not recover scores in the evaluable missing-row subset. The limitations of these analyses—including that they are not full reruns at each alternative setting and do not validate reliability at unobserved deeper levels—are stated explicitly in the corresponding supplementary texts S4 and S5.

   Together, the epxloration of alternative parameters and exploration of filtering steps show that, for 549 signatures (across 61 ligands) both 1) the ligand is clearly identifiable and can be mapped onto the prior knowledge network and 2) the ligand is reachable (within any given number of steps) from downstream TFs. The remaining 810 signatures do not yield an appropriate result to benchmark the ability of MOON scores to capture applied ligands. We therefore do not propose an additional manuscript amendment for R1.8; the R1.6 and R1.7 revisions and their supplementary texts provide the relevant parameter reporting and sensitivity results.

   R1.9: Consistently negative ligands

   Reviewer comment: Four ligands received consistently negative MOON scores across experiments where they were applied (the opposite of the expected direction). The paper identifies and corrects the OSM annotation error but does not investigate or explain the remaining consistently negative ligands. The authors should identify the source of the systematic sign inversion for each of these ligands whether it reflects additional PKN annotation errors, missing interactions, or biological context specificity, and report whether correcting those errors changes the overall benchmark performance.

   Response:

   We thank the reviewer for raising this point. We further explored the four ligands with nominally (i.e. before multiple-hypothesis correction) significant negative mean applied-ligand MOON scores in the submitted CytoSig benchmark: OSM (n = 5, mean = -3.25), FGF10 (n = 9, mean = -1.04), IL6 (n = 39, mean = -0.68), and IGF1 (n = 6, mean = -0.43). These results were identified with the exploratory, unadjusted per-ligand one-sample test used in the benchmark summary. OSM, FGF10, and IGF1 were negative in every corresponding signature; IL6 was negative in 35 of 39 signatures. All four perturbations were annotated as activating treatments. This does not exclude experiment-specific context or data-quality effects, but it does not support a common reversal of the CytoSig treatment-label direction as the explanation for these cases.

   OSM remains the clearest localized PKN error: as reported in the manuscript, the erroneous inhibitory IL6ST annotation was corrected in OmniPath version 2024.03.19. A separate, broad endpoint/bridge PKN diagnostic also changed all five OSM scores from negative to positive (mean -3.25 to 3.63), but this is not the same network comparison as the 2024.03.19 update reported in the manuscript and is not presented as a replacement benchmark. IL6 was strongly rescued in this diagnostic (mean -0.68 to 1.89; negative scores 35/39 to 4/39). Its submitted local network included direct edges to the response genes CRP, A2M, and PTHLH that were absent from the updated representation. Since the diagnostic changes a broad PKN, this supports sensitivity to PKN representation or versioning without assigning the IL6 shift to one interaction.

   FGF10 and IGF1 did not reveal a similarly specific local annotation error. FGF10 remained weakly negative in the diagnostic (mean -1.04 to -0.15; 6/9 scores still negative), and its immediate FGF10–FGFR2 neighborhood was unchanged. All nine submitted FGF10 scores were computed at level 3 through the single preceding-layer node FGFR2. IGF1 was partly rescued (mean -0.43 to 0.49; 2/6 scores still negative), and each submitted score was computed at level 2 through INSR. These cases are compatible with sensitivity to the layer-wise MOON heuristic: an upstream node is scored when first reached through the immediately preceding layer, so evidence available only through more distant downstream paths does not directly re-enter that score. This is a plausible mechanism rather than proof that it is the sole cause of either result. The bounded analysis did not directly test cell-type, protocol, or other biological-context explanations.

   The endpoint/bridge comparison was intentionally treated as a diagnostic rather than a corrected replacement benchmark because it had slightly different coverage (547 scored entries across 61 ligands, compared with 549 entries across 63 ligands in the submitted cache). It contained one nominally significant negative ligand, CD40LG, while the number of nominally significant positive ligands remained 16. Since coverage and the identity of the nominally positive and negative ligands changed, we do not interpret this comparison as demonstrating a net improvement in overall benchmark performance. Instead, the analysis of these four representative ligands supports the conclusion that the negative cases have distinct potential reasons for failure, including PKN representation and layer-wise propagation sensitivity.

   The complete analysis and exploration of numerical results, and limitations of the diagnostic comparison are presented in proposed Supplementary Text S6. To keep the Results section concise, we propose one cross-reference immediately after the initial OSM/IL6 example.

   We have modified the following sections of the manuscript accordingly:

   Location: Results section 2.2

   Original manuscript excerpt:* “...This erroneous annotation also explained the seemingly poor MOON score estimation of the IL6 ligand (Figure 2D), which has an average score of 2.0 in the newer version.”*

   Updated manuscript excerpt:* “...This erroneous annotation also explained the seemingly poor MOON score estimation of the IL6 ligand (Figure 2D), which has an average score of 2.0 in the newer version. Further analyses of the four ligands with nominally significant negative mean MOON scores are presented in Supplementary Text S6.”*

   R1.10: Low positive-ligand reliability

   Reviewer comment: The paper reports only 16 out of 63 ligands (25%) show statistically significant positive MOON scores across experiments. The authors should investigate and discuss what drives this low reliability. Is the variability in scores across experiments for the same ligand explained by cell type specificity, experimental protocol differences (if there are any), PKN incompleteness, or limitations of the MOON scoring procedure itself? The use of a generic PKN like OmniPath, which contains interactions derived from many different cell types and conditions, may introduce context-irrelevant edges that add noise or produce sign inversions when scoring ligands in specific experimental contexts. Context-specific network inference approaches such as ARACNE, GENIE3, SCENIC applied directly to the transcriptomic data of each experiment or to other datasets in similar cell types/conditions, could potentially improve MOON's reliability by restricting propagation to interactions that are actually active in the biological context being analyzed.

   Without understanding the sources of failure it is unclear whether performance can be improved through such strategies or whether the low reliability reflects a more fundamental limitation of the prior knowledge network approach in diverse biological contexts.

   Response:

   The extended exploration reported in Supplementary Text S6 addresses the consistently negative mean-score cases individually: it identifies a localized OSM PKN error, supports sensitivity of IL6 to PKN representation or versioning, and identifies layer-wise propagation bottlenecks as plausible contributors for the unresolved FGF10 and partly IGF1 cases. That analysis does not directly test biological-context or protocol effects, but it shows that the negative-score cases are heterogeneous rather than attributable to a simple common sign reversal.

   To address the reviewer's suggestion of an inferred regulatory network, we performed a separate, bounded IFNA1 diagnostic with GENIE3. We used 104 CytoSig IFNA contrast-statistic profiles to infer a pooled IFNA response-context TF-target layer, while retaining the curated ligand/receptor/signaling PKN. This is not an experiment-specific or cell-type-specific network, and the same IFNA profiles were used for network inference and for evaluation; it is therefore a feasibility diagnostic rather than an independent validation.

   Replacing the curated CollecTRI TF-target layer with GENIE3 did not improve raw applied-ligand IFNA1 MOON scores. IFNA1 was scored in 76 of 104 contrasts in both the submitted CollecTRI/COSMOS cache and the GENIE3-only analysis, but the mean raw score decreased from 7.42 to 4.56 (only 10 of 76 paired scores were higher with GENIE3). A CollecTRI-priority merged layer partly recovered the GENIE3-only score reduction (mean 5.42), but remained below the submitted baseline and improved raw scores in only 18 of 76 paired contrasts. Mean within-experiment rank quantiles were slightly higher for the inferred and merged layers, but this observation is not independent evidence of improved prioritization because the same profiles were used to infer and assess the GENIE3 network.

   This analysis therefore does not show that altering the TF-target layer alone improves ligand-score recovery. It also cannot resolve the R1.9 reasons for failure directly, because the upstream PKN and the layer-wise MOON scoring heuristic were unchanged. Rather, the two analyses support a more specific interpretation: heterogeneous ligand scores can reflect several components of the workflow, including PKN representation, propagation behavior, and the regulatory layer.

   The complete GENIE3 comparison, its numerical results, and its limitations are presented in proposed Supplementary Text S7. To keep the Results section concise, we propose one cross-reference at the end of the existing IFNA1 example, alongside the proposed R1.9 failure-mode analysis.

   We have modified the following sections of the manuscript accordingly:

   Location: Results section 2.2

   Original manuscript excerpt:* “The high MOON scores indicate that treating cell lines with IFNA1 will consistently lead to the activation of IRF9.”*

   Updated manuscript excerpt:* “The high MOON scores indicate that treating cell lines with IFNA1 will consistently lead to the activation of IRF9. An analysis of the impact of a GENIE3-inferred regulatory network on IFNA1 scoring is detailed in Supplementary Text S7, alongside the failure-mode analysis in Supplementary Text S6.”*

   R1 minor comments

   A transcriptional consistency check can be performed that removes any

   Interaction

   number of downstream TF(s) participating

   Check for contracted forms , typically not accepted in written text

   Figure 3 these should be D and E, there are two Cs

   1. C) MOON Network connecting the top deregulated TFs and LR interactions of factor 4 based on a signed directed prior knowledge network. D) Heatmap of the top results of the Pathway control analysis. We represent pathways that are significantly over-represented downstream of given nodes of the MOON thresholded network and merely contextualize(s) a part of it. Check MOON caps

   The precision is again high among the top PC1

   Loading (s)

   the complementarity of MOON scores with

   expression value(s)

   Comma missing

   However, there are many more types of domain knowledge that can potentially be used to interpret feature weights of factors beyond pathway ontologies, such as prior knowledge in the form of footprints and signed-directed networks can help to provide interpretable insights from factor weights.

   Unclear:

   We saw that it was able to recover expected regulation mechanisms, as some of the top gene-pathway interactions were found to be e.g. JAK regulates the JAK-STAT pathway.

   Methods:

   The ULM method of decoupleR was used with CollecTRI to estimate patient specific transcription factor activity signatures (XXX min target per TF).

   Response:

   We thank the reviewer for raising these points. We have corrected the manuscript accordingly.

   Reviewer #1: Significance

   General significance comment

   Reviewer comment: This paper addresses two main issues in the field: 1) going beyond correlation and towards mechanistic and causal explanations in biomedically relevant regulation processes and 2) using prior knowledge while also checking its consistency with data. The approach proposed is likely to provide extremely useful insight, as shown by applications both in-vitro and in patient cohorts.

   The only limitation, which cannot easily be addressed but it is generally an issue in the field, is our lack of certainty about the ground truth which makes it very difficult to assess the relevance of the hypotheses provided by the computational approach.

   The paper will be of broad interest to the computational biology community, especially in oncology and drug discovery.

   We are researchers in the same field, we have tested several other approaches to perform similar tasks.

   Response:

   We thank the reviewer for this positive assessment and agree that uncertainty about biological ground truth is a central limitation of computational-mechanistic inference. COSMOS+ uses signed prior knowledge and multi-omic observations to prioritize mechanistically coherent explanations, but these explanations remain hypotheses whose biological relevance and causality require context-matched validation.

   With the revisions listed above, we aimed to make these limitations more explicit and, where possible, to quantify their impact. These complementary analyses support selected aspects of the framework, but none establishes every inferred regulator or network edge as a causal driver, therapeutic target, or clinically validated biomarker. The scope of the functional and clinical evidence is further detailed under R2.1 and R2.2. We therefore narrowed the corresponding claims and identified targeted perturbation experiments and independent clinical validation as necessary next steps.

   Reviewer #2: Evidence, Reproducibility And Clarity

   R2.1: Computational validation and functional-validation scope

   Reviewer comment: This manuscript presents COSMOS+, an extension of the COSMOS framework that integrates multi-omics factor analysis with prior-knowledge signaling and metabolic networks through the newly developed MOON algorithm.

   The authors demonstrate the approach using the NCI60 dataset, breast cancer resistance models, and a breast cancer patient cohort. The computational framework is technically sophisticated and addresses an important challenge in systems biology, the principal novelty resides in the development of the MOON/COSMOS+ computational framework itself. The biological applications presented throughout the manuscript function largely as case studies illustrating the algorithm rather than generating fundamentally new biological insights. The analyses of breast cancer resistance and the NCI60 dataset are interesting examples, but most conclusions remain computationally inferred and are not experimentally demonstrated.

   The manuscript proposes multiple signaling regulators, resistance-associated pathways, and mechanistic hypotheses, yet none of these predictions are directly tested. Given the emphasis on identifying resistance drivers and actionable biological mechanisms, additional experimental evidence would substantially strengthen the work. Validation through CRISPR-mediated perturbation, knockdown experiments, or other functional assays would help establish whether the inferred network regulators truly contribute to the phenotypes described.

   Response:

   We thank the reviewer for this comment and we agree that the primary contribution of this study is methodological and that COSMOS+/MOON outputs should be interpreted as mechanistically informed, testable hypotheses rather than direct evidence of causality. We have therefore revised the manuscript throughout to avoid describing inferred regulators or pathways as established resistance drivers, causal mechanisms, or clinically actionable targets.

   We also clarify the scope of the existing evaluations. The CytoSig analysis benchmarks recovery of known upstream perturbations, whereas comparison with an orthogonal CRISPR screen in the breast-cancer models provides limited, context-specific support for the prioritised genes. Neither analysis constitutes functional validation of individual network regulators or inferred mechanisms.

   We now explicitly acknowledge that targeted perturbation and, where appropriate, rescue experiments in the relevant model systems will be needed to establish causal roles for specific candidates. We have revised the Abstract, Introduction, Results, Methods, and Discussion accordingly.

   We have modified the following sections of the manuscript accordingly:

   Location: Abstract

   Original manuscript excerpt:* “We apply this approach on a novel multi-omics dataset of cell line models of breast cancer resistance to evaluate the ability of such mechanistic hypotheses to identify resistance drivers, as well as a breast cancer patient cohort. Our approach offers an interpretable framework to generate actionable insights from multi-omic data particularly suited for high dimensional datasets.”*

   Updated manuscript excerpt:* “We apply this approach on a novel multi-omics dataset of cell line models of breast cancer resistance to explore resistance-associated mechanistic hypotheses, as well as a breast cancer patient cohort. Our approach offers an interpretable framework to generate mechanistically informed hypotheses from multi-omic data particularly suited for high dimensional datasets.”*

   Location: Introduction

   Original manuscript excerpt:* “We show that mechanistic hypotheses of signaling deregulation in resistant and sensitive cell lines treated with CDK inhibitors are correlated with resistance markers identified through knock-out screening.”*

   Updated manuscript excerpt:* “We show that mechanistic hypotheses of signaling deregulation in resistant and sensitive cell lines treated with CDK inhibitors are modestly correlated with gene sensitization profiles obtained through knock-out screening.”*

   Location: Results, Section 2.3

   Original manuscript excerpt:* “We can interpret such scores as TFs that are responsible for the transcriptional programs that are captured by given factors.”*

   Updated manuscript excerpt:* “We can interpret such scores as TF activities that are consistent with the transcriptional programs captured by given factors.”*

   Location: Results, Section 2.4

   Original manuscript excerpt:* “Therefore, the network associated with MCF7SYL recapitulates well the CDK2 mediated acquired resistance to treatment while highlighting a potential positive feedback loop through MYC and the NOTCH pathway.”*

   Updated manuscript excerpt:* “Therefore, the network associated with MCF7SYL is consistent with the hypothesis that CDK2 contributes to acquired resistance to treatment while highlighting a potential positive feedback loop through MYC and the NOTCH pathway.”*

   Location: Results, Section 2.5

   Original manuscript excerpt:* “We sought to evaluate their accuracy and functional relevance with respect to resistance to treatment by investigating whether feature weights of the PCA on MOON scores were associated with mechanisms of resistance to the CDK inhibitors CDK4/6i (CDK4/6) and CDK2/4/6i (CDK2/4/6).”*

   Updated manuscript excerpt:* “We assessed the concordance between feature weights of the PCA on MOON scores and gene sensitization profiles obtained in the presence of the CDK inhibitors CDK4/6i (CDK4/6) and CDK2/4/6i (CDK2/4/6).”*

   Location: Results, Section 2.5

   Original manuscript excerpt:* “Thus, while small, the correlation patterns are significantly consistent across every factor for the resistant cell line but not for the sensitive one. This agrees with the expectation that PCA factors that are specifically associated with resistant and sensitive cell separation would potentially be associated with resistance mechanisms that can be sensitization targets, but only in resistant cells (HCC1806) and not cells that are already sensitive (MCF7).”*

   Updated manuscript excerpt:* “Thus, while small, the relationship between the factor-specific MOON–crispR correlations and the separation between resistant and sensitive cells was significant for HCC1806 but not for MCF7. This indicates concordance between MOON-based prioritization and crispR sensitization profiles in this comparison, but does not establish that individual prioritized genes or inferred network edges mediate resistance.”*

   Location: Results, Section 2.5

   Original manuscript excerpt:* “...we first considered genes that had a KO sensitization score of -2 at least as true positives, that is genes that are driving resistance to CDK4/6i or CDK2/4/6i. We could then compute the area under the precision curve (AUPRC) for genes that have negative weights in PC1 (that is, genes that are more active in the resistant cell line after treatment with CDK4/6i and CDK2/4/6i). The AUPRC is favored in this case to the area under the receiving operator curve (AUROC) due to 1) the imbalance of the ratio of true positive and true negative (only 7% of true positive) and 2) precision is a usual metric for a model that we assume does not inherently capture the full complexity of the underlying biological mechanism, and merely contextualize a part of it.[...] Therefore, while the AUPRC is higher than a random baseline, this shows that the precision of the PC1 weight is better than the recall to capture sensitization drivers.”*

   Updated manuscript excerpt:* “...we considered genes with a KO sensitization score of -2 or lower as positive instances for the precision-recall analysis. We then computed the area under the precision-recall curve (AUPRC) for genes that have negative weights in PC1 (that is, genes that are more active in the resistant cell line after treatment with CDK4/6i and CDK2/4/6i). The AUPRC is favored over the area under the receiver operating characteristic curve (AUROC) because the positive instances are imbalanced (7% of genes) and because the MOON score PCA is not expected to capture the full complexity of the underlying biological mechanisms.[...] Therefore, while the AUPRC is higher than a random baseline, the PC1 weights had higher precision than recall for identifying genes with strong crispR sensitization scores.”*

   Location: Methods, Section 4.10

   Original manuscript excerpt:* “4.10 Analysis of crispR KO data and validation of resistance mechanisms”*

   Updated manuscript excerpt:* “4.10 Analysis of crispR KO data and comparison with resistance-associated MOON features”*

   Location: Methods, Section 4.10

   Original manuscript excerpt:* “It can then be interpreted as e.g. an indication of how much does a gene that is more active in a resistant cell line treated with a CDK inhibitor such as CDK4/6i actually is directly responsible for the resistance to the treatment.”*

   Updated manuscript excerpt:* “It can then be interpreted as an indication of the concordance between resistance-associated MOON features and crispR sensitization scores, rather than as evidence that the corresponding genes are directly responsible for resistance to treatment.”*

   Location: Discussion

   Original manuscript excerpt:* “We show how MOON can identify biological mechanisms underlying cancer treatment resistance.”*

   Updated manuscript excerpt:* “We show how MOON can be used to explore biological mechanisms associated with cancer treatment resistance.”*

   Location: Discussion

   Original manuscript excerpt:* “Therefore, we also assessed the ability of mechanistic hypotheses generated by COSMOS+ to support the identification of drivers of treatment resistance by applying it on a novel combined transcriptomic and phospho-proteomic dataset of breast cancer cell lines with different resistance profiles exposed to CDK inhibitor drugs measured at early (2-4 hours) and late time points (72-96 hours).”*

   Updated manuscript excerpt:* “Therefore, we applied COSMOS+ to a novel combined transcriptomic and phospho-proteomic dataset of breast cancer cell lines with different resistance profiles exposed to CDK inhibitor drugs measured at early (2-4 hours) and late time points (72-96 hours), and compared the resulting resistance-associated MOON prioritizations with crispR sensitization profiles.”*

   Location: Discussion,

   Original manuscript excerpt:* “In this context, we sought to use COSMOS+ to examine both known and potentially novel mechanistic hypotheses mediating response and resistance to CDK4/6i. COSMOS+ recapitulated known signaling and transcriptional regulation components following CDK inhibition such as the CDK2,4,6, RB1 and E2Fs crosstalk, while also highlighting a wide set of other potentially important mechanisms that were differentially regulated between the sensitive and resistant cell line, such as a YWHAQ/E2F1/FOXO3 and a CDK2/NOTCH crosstalk.”*

   Updated manuscript excerpt:* “In this context, we sought to use COSMOS+ to examine both known and potentially novel mechanistic hypotheses associated with response and resistance to CDK4/6i. COSMOS+ highlighted expected signaling and transcriptional regulation components following CDK inhibition such as the CDK2,4,6, RB1 and E2Fs crosstalk, while also highlighting a wide set of other candidate mechanisms that were differentially regulated between the sensitive and resistant cell line, such as a YWHAQ/E2F1/FOXO3 and a CDK2/NOTCH crosstalk.”*

   Location: Discussion

   Original manuscript excerpt:* “We then compared the mechanistic hypothesis identified [...]. The poor recall suggests that there is a large gap between the ability to identify deregulated signaling and gene regulation between different cell lines treated with a drug and the actual identification of direct drivers of drug resistance. [...] We also showed how pathway control analysis could correctly recapitulate the control of E2Fs transcription factors over cell cycle processes in a breast cancer cell line dataset and how it is differentially regulated between sensitive and resistant cell lines treated with CDK inhibitors.”*

   Updated manuscript excerpt:* “We then compared the mechanistic hypotheses generated [...]. The poor recall highlights the gap between identifying deregulated signaling and gene regulation between different cell lines treated with a drug and identifying genes with drug-specific sensitization in this assay. [...] We also showed how pathway control analysis highlighted the expected control of E2F transcription factors over cell cycle processes in a breast cancer cell line dataset and how it is differentially regulated between sensitive and resistant cell lines treated with CDK inhibitors.”*

   Location: Discussion

   Original manuscript excerpt:* “... this hypothesis can be validated and translated into actionable insights. A common example of such actionable insight in pharmacological research is the identification of new drug targets, especially for cells and tissues that are resistant to existing treatments. [...] In this study, we show that there is partial consistency between differentially regulated molecular drug responses and molecular drivers of treatment sensitization. Furthermore, the good precision of COSMOS+ to capture sensitization targets but poor recall is consistent with the idea that such a model, while being able to accurately capture important biological mechanisms, is not yet able to generate predictive mechanistic models of cell signaling.”*

   Updated manuscript excerpt:* “... it requires targeted experimental validation in the relevant biological context before it can be considered a causal or actionable mechanism. A common potential downstream application in pharmacological research is the identification of new drug targets, especially for cells and tissues that are resistant to existing treatments. [...] In this study, we show that there is partial consistency between differentially regulated molecular drug responses and genes with crispR sensitization scores. Furthermore, the higher precision than recall of COSMOS+ for genes with crispR sensitization scores is consistent with the idea that such a model can prioritize candidate mechanisms for experimental testing, but is not yet able to generate predictive mechanistic models of cell signaling.”*

   R2.2: Clinical and translational scope

   Reviewer comment: I also have reservations regarding the translational implications proposed by the authors. The manuscript suggests potential utility for patient outcome prediction and resistance mechanism discovery. However, the patient cohort analyzed is relatively limited in size and, importantly, no independent external validation cohort is provided. As a result, it remains unclear how robust and generalizable these predictive signatures would be in broader clinical settings. This aspect should be discussed more cautiously.

   Response:

   We thank the reviewer for this important comment and agree that the PALOMA3 analysis requires a more cautious clinical interpretation. The analysed PALOMA3 RNA cohort is relatively limited in size (302 patients) and no independent external validation cohort was available. This analysis was not intended to develop or validate a clinically deployable progression-free-survival predictor. Its purpose was methodological: to examine, within one cohort, whether MOON-derived activity estimates capture outcome associations that are complementary to those obtained from RNA abundance alone.

   In particular, the Cox coefficients used to construct the RNA, MOON, and hybrid scores, and the correlations of those scores with observed PFS, were derived and assessed in the same PALOMA3 cohort. The higher correlation of the hybrid score is therefore an in-cohort illustration of potential complementarity between RNA and MOON representations. It does not quantify out-of-sample predictive performance, incremental clinical value, or generalisability to other breast-cancer cohorts.

   As detailed under R1.4, no individual MOON score remained significant after Benjamini-Hochberg correction in the separately evaluated treatment-arm node-wise Cox scan. The RB1 observation is accordingly presented as a nominal, hypothesis-generating illustration rather than an individually validated prognostic or predictive biomarker. Independent validation in larger, clinically comparable cohorts would be required before any RNA, MOON, or hybrid signature could be considered for outcome prediction, treatment selection, or clinical biomarker use.

   We therefore amended the Introduction, Results, and Discussion. They retain the PALOMA3 analysis as an exploratory within-cohort comparison of molecular representations, while removing language that implies validated clinical prediction, clinical biomarker identification, or generalisable treatment-efficacy prediction. The Figure 6 legend is retained unchanged.

   We have modified the following sections of the manuscript accordingly:

   Location: Introduction

   Original manuscript excerpt:* “Finally, we apply the COSMOS+ framework to a breast cancer patient cohort, and we demonstrate the complementarity of the resulting features with omic data to predict patient outcomes.”*

   Updated manuscript excerpt:* “Finally, we apply the COSMOS+ framework to a breast cancer patient cohort to examine the complementarity of the resulting features with omic data in associations with patient outcomes.”*

   Location: Results, Section 2.6

   Original manuscript excerpt:* “We also evaluated if the PC1 loadings from the cell line data could be used as a signature to predict poor patient response.”*

   Updated manuscript excerpt:* “We also explored whether the PC1 loadings from the cell line data were associated with patient response in the PALOMA3 treatment arm.”*

   Original manuscript excerpt:* “This showed that there is a significant amount of information that can be extracted from a sensitive/resistance treatment response cell line model to inform patient response.”*

   Updated manuscript excerpt:* “Because the signature size and its association with PFS were explored in the same cohort, this result is descriptive and does not establish a clinically generalisable prediction signature.”*

   Location: Results, Section 2.6

   Original manuscript excerpt:* “Finally, to further explore the complementarity of MOON scores with expression value, we created a hybrid Cox coefficient signature by combining the most extreme cox coefficients computed from the RNA data on the one hand and from MOON scores on the other hand (Figure 6G). That is, for each gene, the cox coefficient that was the most extreme between MOON and RNA was included. We then evaluated the ability of such a signature to estimate patient outcome by multiplying the hybrid Cox coefficient signature with the corresponding hybrid patient cohort of MOON scores and RNA measurements. The resulting vector essentially corresponds to scaled PFS predictions. The scaled PFS prediction was compared with the actual PFS for the hybrid signature, as well as with a signature based on RNA or MOON scores alone. The correlation for the hybrid signature PFS prediction was 0.45, while the RNA alone and MOON scores alone were 0.4 and 0.35, respectively (Figure 6H). This further supports that both MOON score and RNA values can bring complementary information to inform us about gene regulation events and processes that may be associated with better or worse outcomes in patients.”*

   Updated manuscript excerpt:* “Finally, to further explore the complementarity of MOON scores with expression values, we created a hybrid Cox coefficient signature by combining the most extreme Cox coefficients computed from RNA data and MOON scores (Figure 6G). For each gene, the most extreme coefficient between the MOON and RNA analyses was included. We then compared the association of this hybrid score with observed patient outcomes, using the same PALOMA3 cohort to derive the Cox coefficients and assess the resulting scores. The correlation with observed PFS was 0.45 for the hybrid score, compared with 0.40 and 0.35 for RNA-only and MOON-only scores, respectively (Figure 6H). This within-cohort comparison is intended to illustrate the potential complementarity of RNA and MOON features; it does not assess out-of-sample predictive performance or establish a generalisable model for PFS prediction.”*

   Location: Discussion

   Original manuscript excerpt:* “These results also illustrate how mechanistic hypotheses generated in the context of pre-clinical cell-line models can translate into a clinical setting and potentially help predict treatment efficacy as well as identify underlying molecular causes of treatment resistance.”*

   Updated manuscript excerpt:* “These results also illustrate how mechanistic hypotheses generated in the context of pre-clinical cell-line models can be compared with outcome associations in a clinical cohort, but do not establish generalisable prediction of treatment efficacy or the underlying molecular causes of treatment resistance.”*

   Location: Discussion

   Original manuscript excerpt:* “Finally, we demonstrated how COSMOS+ can be used with clinical baseline patient cohort data to complement biomarker predictions based on transcriptomics data alone. We showed how a Cox survival analysis performed with MOON scores instead of transcriptomic data allowed for the identification of expected markers that were missed by transcriptomic data alone. We then showed that the Cox coefficients based on MOON scores were correlated with MOON scores estimated from a cell line model of resistant and sensitive cell lines treated with similar drug regimens, suggesting that mechanisms identified in such a cell line model could support the identification of clinical biomarkers. While small (Pearson correlation coefficients = -0.12), the correlation was highly significant (p-value *

   Updated manuscript excerpt:* “Finally, we used baseline data from the PALOMA3 patient cohort to explore whether MOON-derived activity estimates captured outcome associations complementary to transcriptomic data alone. The correlation between the cell-line PC1 loadings and MOON-score Cox coefficients (Pearson correlation coefficient = -0.12, p-value *

   Location: Discussion

   Original manuscript excerpt:* “Nonetheless, the mechanisms hypothesized by COSMOS+ seem to complement the information in the omic dataset alone. Indeed, when COSMOS+ was applied at the level of the transcriptome measurements of the breast cancer cohort Paloma3 patients, it revealed markers of treatment outcomes that could not have been captured at the level of their expression alone. For example, RB1 activity estimated from its expected downstream regulated expression target was significantly associated with the worst outcome in the treatment arm of the Paloma3 cohort, while its own expression was marginally associated with the outcome in the wrong direction (that is, high expression of the RB1 tumor suppressor associated with worst outcome).”*

   Updated manuscript excerpt:* “Nonetheless, the analysis of the PALOMA3 cohort suggests that MOON scores may capture outcome associations complementary to RNA expression. For example, the RB1 MOON score was nominally associated with worse outcome in the treatment arm, whereas its expression was not associated with outcome in the same direction (see Results); this is presented as a hypothesis-generating illustration rather than an individual prognostic biomarker.”*

   Location: Discussion

   Original manuscript excerpt:* “This complementarity is further illustrated by the increased correlation between a PFS prediction performed with RNA or MOON scores alone compared to a hybrid signature. Therefore, it will be interesting in the future to assess if COSMOS+ could be used in combination with omic features to improve the performance of predictive approaches.”*

   Updated manuscript excerpt:* “This complementarity is further illustrated by the higher within-cohort correlation with observed PFS of the hybrid RNA/MOON score than of the RNA- or MOON-only scores. Because the Cox coefficients and scores were derived and assessed in the same relatively limited PALOMA3 cohort, this comparison does not establish the performance or generalisability of a PFS prediction model; independent validation in larger, clinically comparable cohorts will be required.”*

   R2.3: Transcriptomic, prior-knowledge, and interpretation biases

   Reviewer comment: Suggest to have a look at paper such as "Technical and Biological Biases in Bulk Transcriptomic Data Mining for Cancer Research" which discuss database bias.

   The authors themselves provide a compelling example in which an erroneous IL6ST annotation leads to incorrect MOON predictions. This example highlights a broader issue: the quality of the inferred mechanistic networks is inherently constrained by the accuracy and completeness of resources such as OmniPath, STITCH, and related databases. Consequently, false or incomplete annotations may propagate through the analytical pipeline and influence biological interpretation.

   Response:

   We thank the reviewer for highlighting this important limitation. We agree that the accuracy and completeness of the prior-knowledge resources constrain the mechanistic hypotheses generated by COSMOS+. The OSM–IL6ST example is a clear illustration that an incorrect interaction can propagate through the network and affect the resulting scores and interpretation.

   As detailed under R1.9, we explored the four ligands with consistently negative MOON scores. This distinguishes the localized OSM–IL6ST annotation error from broader sensitivity of the IL6 result to network representation, while no comparably specific annotation error was identified for FGF10 or IGF1. The diagnostic comparison is not presented as a replacement benchmark and does not exclude biological-context or other experimental explanations.

   Consistent with the claim-scope revisions under R2.1 and R2.2, we therefore interpret MOON scores and subnetworks as prior-knowledge-dependent, mechanistically informed hypotheses rather than definitive mechanisms, causal drivers, or clinical biomarkers. Furthermore, we believe such an approach can in fact help pin-pointing problems in prior knowledge and help fix them, as in the OSM case. No additional manuscript-text amendment is proposed for R2.3: the detailed ligand-level analysis is addressed under R1.9, and the applied R2.1/R2.2 amendments already make this interpretation boundary explicit.

   R2.4: Positioning relative to integrative oncology studies

   Reviewer comment: The authors may wish to discuss how COSMOS+ compares conceptually with recent efforts aimed at identifying clinically relevant oncogenic determinants through integrative multi-omics analyses. The article "Comprehensive analysis of oncogenic determinants across tumor types via multi-omics integration" would be a valuable addition and would help position the present framework within the broader landscape of cancer systems biology and biomarker discovery.

   Response:

   We thank the reviewer for this helpful suggestion. We added a Discussion paragraph that positions COSMOS+ relative to complementary integrative oncology approaches. Pan-cancer studies combine multiple molecular layers to characterize recurrent oncogenic alterations and molecular subtypes, whereas COSMOS+ combines activity estimates with a signed prior knowledge network to formulate context-specific, testable mechanistic hypotheses. We make explicit that COSMOS+ is not intended to establish recurrent genomic drivers or validated clinical biomarkers.

   We cite the suggested article as a recent broad review of this landscape, together with two primary TCGA studies that directly support the respective recurrent-alteration and molecular-subtype statements. The functional-validation and clinical-interpretation boundaries are already specified in the accepted R2.1 and R2.2 amendments.

   Citation additions: Sanchez-Vega et al (2018), Cell, doi:10.1016/j.cell.2018.03.035; and Ubaid et al (2025), Cancer Genetics, doi:10.1016/j.cancergen.2025.08.010.

   We have modified the following sections of the manuscript accordingly:

   Location: Discussion

   Additional text in manuscript:* “These analyses demonstrate the flexibility and efficiency of COSMOS+ to generate mechanistic hypotheses across multi-omic layers and in diverse contexts. These results also illustrate how mechanistic hypotheses generated in the context of pre-clinical cell-line models can be compared with outcome associations in a clinical cohort, but do not establish generalisable prediction of treatment efficacy or the underlying molecular causes of treatment resistance.*** ** Pan-cancer studies have integrated multiple molecular layers to characterize recurrent oncogenic alterations and molecular subtypes (Sanchez-Vega et al, 2018; Hoadley et al, 2018) Ubaid et al (2025). COSMOS+ is orthogonal to these approaches, as it combines activity estimates with a signed prior knowledge network to formulate context-specific, testable mechanistic hypotheses rather than to establish recurrent genomic drivers or validated clinical biomarkers. This opens the possibility to combine both approaches, using functional mutations as inputs for COSMOS+.”

   Reviewer #2: Significance

   R2.S1: Validation scope and significance

   Reviewer comment: I find the methodological contribution potentially valuable, but the current manuscript remains largely computational and exploratory. Additional functional validation and stronger clinical validation would considerably enhance the biological significance and translational impact of the study.

   Response:

   We agree that additional targeted, context-matched perturbation or rescue experiments and independent clinical validation would strengthen any individual mechanistic or translational claim. The primary contribution of this study is indeed computational, and in particular methodological. The relevant clarifications of scope are detailed above (R2.1 and R2.2). These amendments make the functional and clinical validation boundaries explicit. We did not add further wet-lab validations or an independent clinical cohort in this revision; such validation would need to be designed separately for each prioritized mechanism and biological context and we estimated that they would fall outside of the clarified scope of this manuscript.

   Reviewer #3: Evidence, Reproducibility And Clarity

   General comment

   Reviewer comment: The model and datasets are clearly presented for the reproducibility and clarity.

   Response:

   We thank Reviewer 3 for recognizing the clear presentation of the model and datasets.

   Reviewer #3: Significance

   R3.1: Sparse biologically meaningful subnetworks

   Reviewer comment: It is a good idea to integrate multi-omics data to mine the potential signaling pathways. One challenge is that the network is complex and it is hard to identify a small scale, sparse but biologically meaningful signaling pathways/networks. It might be helpful to provide solid results or discussion about this challenge.

   Response:

   We thank the reviewer for highlighting this important challenge. We agree that reducing complex prior-knowledge networks to small, sparse, and biologically interpretable subnetworks remains difficult, particularly because alternative network paths can explain the same multi-omic observations.

   COSMOS+ already provides a set of complementary, user-accessible functions intended to improve interpretation at different scales. For a selected upstream node, get_moon_scoring_network extracts a focused scoring network containing the downstream footprint used to compute its MOON score, allowing the local evidence underlying an individual prioritization to be examined. For a compact solution-level view, reduce_solution_network applies a user-chosen absolute-score cutoff, retains only sign-consistent interactions, restricts the network to nodes reachable from the specified upstream inputs, and prunes unsupported sources and sinks. Pathway-control analysis provides a complementary pathway-level summary: it tests which pathways are over-represented among the downstream neighbourhoods of high-scoring nodes and represents the results in a node--pathway matrix. In the manuscript, these approaches are illustrated through score-specific ligand networks and through the thresholded factor-4 MOON network together with its pathway-control analysis.

   We further provide additional functions in the COSMOS NCI60 tutorial; here a dual-threshold network-extraction procedure. A stringent primary threshold defines a high-scoring core, while a lower secondary threshold retains score-supported neighbouring nodes of that core; sign-consistency filtering is then applied. Although this procedure is not used for the networks shown in the main manuscript, it provides a practical analysis option for obtaining compact, interpretable subnetworks while preserving locally supported context.

   These functions improve the interpretability of complex network results but do not provide a general solution to sparse-network inference. Developing and evaluating more robust strategies to produce sparse, biologically meaningful network summaries remains an active area of our ongoing work.

   We have added the following sections to the manuscript accordingly:

   Location: Discussion

   Additional text in manuscript:* “In order to pin-pointing such errors more easily, we have developed a range of functions to help extract and interpret sparse subnetworks for the contextualised prior knowledge. For example, some functions help extract sub-networks based on minimum activity scores of nodes, other focus on connecting a specific node to downstream measurements that it relies on for its score, and a function to explore over-representation of biological terms and pathways downstream of a given node.”*

   Read the original source
2. 

   Review Commons

   Aug 17, 2026

   Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

   Learn more at Review Commons

   ---

   Referee #3

   Evidence, reproducibility and clarity

   The model and datasets are clearly presented for the reproducibility and clarity.

   Significance

   It is a good idea to integrate multi-omics data to mine the potential signaling pathways. One challenge is that the network is complex and it is hard to identify a small scale, sparse but biologically meaningful signaling pathways/networks. It might be helpful to provide solid results or discussion about this challenge.

   Read the original source
3. 

   Review Commons

   Aug 17, 2026

   Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

   Learn more at Review Commons

   ---

   Referee #2

   Evidence, reproducibility and clarity

   This manuscript presents COSMOS+, an extension of the COSMOS framework that integrates multi-omics factor analysis with prior-knowledge signaling and metabolic networks through the newly developed MOON algorithm.

   The authors demonstrate the approach using the NCI60 dataset, breast cancer resistance models, and a breast cancer patient cohort. The computational framework is technically sophisticated and addresses an important challenge in systems biology，the principal novelty resides in the development of the MOON/COSMOS+ computational framework itself. The biological applications presented throughout the manuscript function largely …

   Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

   Learn more at Review Commons

   ---

   Referee #2

   Evidence, reproducibility and clarity

   This manuscript presents COSMOS+, an extension of the COSMOS framework that integrates multi-omics factor analysis with prior-knowledge signaling and metabolic networks through the newly developed MOON algorithm.

   The authors demonstrate the approach using the NCI60 dataset, breast cancer resistance models, and a breast cancer patient cohort. The computational framework is technically sophisticated and addresses an important challenge in systems biology，the principal novelty resides in the development of the MOON/COSMOS+ computational framework itself. The biological applications presented throughout the manuscript function largely as case studies illustrating the algorithm rather than generating fundamentally new biological insights. The analyses of breast cancer resistance and the NCI60 dataset are interesting examples, but most conclusions remain computationally inferred and are not experimentally demonstrated.

   The manuscript proposes multiple signaling regulators, resistance-associated pathways, and mechanistic hypotheses, yet none of these predictions are directly tested. Given the emphasis on identifying resistance drivers and actionable biological mechanisms, additional experimental evidence would substantially strengthen the work. Validation through CRISPR-mediated perturbation, knockdown experiments, or other functional assays would help establish whether the inferred network regulators truly contribute to the phenotypes described.

   I also have reservations regarding the translational implications proposed by the authors. The manuscript suggests potential utility for patient outcome prediction and resistance mechanism discovery. However, the patient cohort analyzed is relatively limited in size and, importantly, no independent external validation cohort is provided. As a result, it remains unclear how robust and generalizable these predictive signatures would be in broader clinical settings. This aspect should be discussed more cautiously. Suggest to have a look at paper such as ""Technical and Biological Biases in Bulk Transcriptomic Data Mining for Cancer Research" which discuss database bias.

   The authors themselves provide a compelling example in which an erroneous IL6ST annotation leads to incorrect MOON predictions. This example highlights a broader issue: the quality of the inferred mechanistic networks is inherently constrained by the accuracy and completeness of resources such as OmniPath, STITCH, and related databases. Consequently, false or incomplete annotations may propagate through the analytical pipeline and influence biological interpretation.the authors may wish to discuss how COSMOS+ compares conceptually with recent efforts aimed at identifying clinically relevant oncogenic determinants through integrative multi-omics analyses. The article "Comprehensive analysis of oncogenic determinants across tumor types via multi-omics integration" would be a valuable addition and would help position the present framework within the broader landscape of cancer systems biology and biomarker discovery.

   Significance

   I find the methodological contribution potentially valuable, but the current manuscript remains largely computational and exploratory. Additional functional validation and stronger clinical validation would considerably enhance the biological significance and translational impact of the study.

   Read the original source
4. 

   Review Commons

   Aug 17, 2026

   Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

   Learn more at Review Commons

   ---

   Referee #1

   Evidence, reproducibility and clarity

   Dugourd et al. present COSMOS+, a computational framework that integrates multi-omics data with prior knowledge networks to generate mechanistic hypotheses connecting signaling, transcriptional regulation, and metabolism. The study addresses a relevant challenge in the field: the difficulty of moving beyond purely data-driven factor analysis toward biologically interpretable and causally grounded insights, introducing MOON, a scalable iterative network scoring algorithm as a practical alternative to computationally expensive optimization-based approaches such as CARNIVAL. Importantly, the method has the potential to highlight errors …

   Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

   Learn more at Review Commons

   ---

   Referee #1

   Evidence, reproducibility and clarity

   Dugourd et al. present COSMOS+, a computational framework that integrates multi-omics data with prior knowledge networks to generate mechanistic hypotheses connecting signaling, transcriptional regulation, and metabolism. The study addresses a relevant challenge in the field: the difficulty of moving beyond purely data-driven factor analysis toward biologically interpretable and causally grounded insights, introducing MOON, a scalable iterative network scoring algorithm as a practical alternative to computationally expensive optimization-based approaches such as CARNIVAL. Importantly, the method has the potential to highlight errors in prior knowledge networks. This will not palliate the incompleteness of the existing prior knowledge but, at least, it can identify inconsistencies and help 'correct' the databases. It is not entirely clear whether inconsistencies stem from specificities of cell lines/samples that are not reflected in the general databases, or simply mistakes in the databases. In any case, reconstructing mechanistic hypotheses based on (multi)omics datasets remains an important endeavour, necessary to understand biological processes and also with clear applications in medicine. The code is available and clearly documented ensuring reproducibility. The literature survey is comprehensive and useful to understand the need for developments in this field.

   While the framework is flexible and the applications span a diverse range of contexts from cell line collections to clinical cohorts, some aspects of the work require stronger justification and additional evaluation tests might increase the usefulness of the methods for the community.

   1. MOON scores are ULM t-values at layer 1, but at deeper layers they are t-values computed from previous t-values. It is unclear whether the propagated scores follow the same distribution across layers, which is a requirement to apply the same threshold of |score| > 1.5 uniformly. The authors should either provide a justification for why the threshold remains meaningful at deeper layers, or characterize through simulation how score distributions change across layers and propose a layer-specific thresholding strategy.
   2. The transcriptional consistency check removes TF -> target edges where TF score is incoherent with target expression, but no analogous check is applied to kinase -> TF or receptor -> kinase edges. The upstream signaling layer is therefore only constrained by the upstream anchor comparison while the transcriptional layer (TF -> target gene) undergoes aggressive pruning. The authors should justify this asymmetry explicitly or extend the consistency check to upstream layers where phosphoproteomic data is available. We are aware that prior knowledge on signs of activation might also be lacking, often it will be hard to know what is the real ground truth without doing specific experiments in the right cell types/conditions. Some more discussion about this would clarify the difficulty of this problem.
   3. For the MOFA section, the authors state they used the parameter viewsscale_views = False. This means the three omics layers were not variance-normalized before input to MOFA. Since transcriptomics, proteomics, and metabolomics have different features ranges, the omic layer with higher absolute variance will dominate the factor structure. The authors should either justify this choice explicitly or show that the factor structure is not dominated by a single omic layer.
   4. There is no multiple testing correction in the clinical association analysis. The ULM-based clinical metadata association tests each factor against each clinical category independently. With 9 factors and the number of clinical categories available in NCI60 (tissue of origin alone has ~10 categories, plus age, pathology, and other variables), the number of simultaneous tests is large. The same issue applies to the Cox survival analysis in the Paloma3 section, where a separate Cox model is fitted for every node in the MOON network and results are reported without any correction for multiple comparisons. The authors should apply FDR correction at both stages to confirm that reported associations are not false positives arising from the large number of simultaneous tests performed.
   5. The threshold of |t| > 2 is applied uniformly across TF activity scoring, kinase scoring, LR scoring, and clinical associations without justification. A t-value of 2 corresponds approximately to p < 0.05 only for specific sample sizes, and the implied p-value therefore differs across the three datasets analyzed in this paper which have substantially different sample sizes. The authors should justify this threshold, perform a sensitivity analysis, or demonstrate that conclusions are robust to alternative thresholds.
   6. The methods section does not report what minimum regulon size was used for the Cytosig analysis, and no justification is provided for this parameter. It is not clear if it is still 10 or a different one. Cytosig signatures vary widely in how many genes were measured across experiments, meaning a uniform minimum threshold will exclude TFs not because they are inactive but simply because their targets were not measured in a given experiment. The authors should explicitly report the threshold used, assess whether using a lower threshold such as 5 recovers additional scorable signatures without substantially degrading TF activity reliability, and report the distribution of gene coverage across Cytosig signatures to contextualize how many signatures are affected by this limitation.
   7. The maximum number of propagation steps for the reachability filtering is not reported or justified. In the case of the Cytosig analysis, the paper states that 31 ligands could not be scored because they were not reachable upstream of TFs within the allowed number of steps, but never states what that number was. This parameter directly determines which ligands can be benchmarked, but no sensitivity analysis is provided showing whether increasing the step limit recovers additional ligands or changes the benchmark results. The tool might be enhanced with a report of alternative values, discussion of optimal values and discussion of the tradeoff between reachability and the MOON score reliability at greater network distances.
   8. The benchmark scores only 549 out of 1359 available signatures (40% of the data) due to PKN reachability and TF coverage limitations. Since both limiting parameters are neither reported nor optimized, it is unclear whether the excluded 60% represents a fundamental limitation of the method or an artifact of conservative parameter choices. The authors should explore alternative parameter values and report their effect on both benchmark coverage and performance jointly.
   9. Four ligands received consistently negative MOON scores across experiments where they were applied (the opposite of the expected direction). The paper identifies and corrects the OSM annotation error but does not investigate or explain the remaining consistently negative ligands. The authors should identify the source of the systematic sign inversion for each of these ligands whether it reflects additional PKN annotation errors, missing interactions, or biological context specificity, and report whether correcting those errors changes the overall benchmark performance.
   10. The paper reports only 16 out of 63 ligands (25%) show statistically significant positive MOON scores across experiments. The authors should investigate and discuss what drives this low reliability. Is the variability in scores across experiments for the same ligand explained by cell type specificity, experimental protocol differences (if there are any), PKN incompleteness, or limitations of the MOON scoring procedure itself? The use of a generic PKN like OmniPath, which contains interactions derived from many different cell types and conditions, may introduce context-irrelevant edges that add noise or produce sign inversions when scoring ligands in specific experimental contexts. Context-specific network inference approaches such as ARACNE, GENIE3, SCENIC applied directly to the transcriptomic data of each experiment or to other datasets in similar cell types/conditions, could potentially improve MOON's reliability by restricting propagation to interactions that are actually active in the biological context being analyzed. Without understanding the sources of failure it is unclear whether performance can be improved through such strategies or whether the low reliability reflects a more fundamental limitation of the prior knowledge network approach in diverse biological contexts.

   Minor comments:

   A transcriptional consistency check can be performed that removes any Interaction number of downstream TF(s) participating

   Check for contracted forms , typically not accepted in written text

   Figure 3 these should be D and E, there are two Cs C) MOON Network connecting the top deregulated TFs and LR interactions of factor 4 based on a signed directed prior knowledge network. D) Heatmap of the top results of the Pathway control analysis. We represent pathways that are significantly over-represented downstream of given nodes of the MOON thresholded network

   and merely contextualize(s) a part of it.

   Check MOON caps The precision is again high among the top PC1 Loading (s)

   the complementarity of MOON scores with expression value(s)

   Comma missing However, there are many more types of domain knowledge that can potentially be used to interpret feature weights of factors beyond pathway ontologies, such as prior knowledge in the form of footprints and signed-directed networks can help to provide interpretable insights from factor weights.

   Unclear: We saw that it was able to recover expected regulation mechanisms, as some of the top gene-pathway interactions were found to be e.g. JAK regulates the JAK-STAT pathway.

   Methods: The ULM method of decoupleR was used with CollecTRI to estimate patient specific transcription factor activity signatures (XXX min target per TF).

   Significance

   This paper addresses two main issues in the field: 1) going beyond correlation and towards mechanistic and causal explanations in biomedically relevant regulation processes and 2) using prior knowledge while also checking its consistency with data. The approach proposed is likely to provide extremely useful insight, as shown by applications both in-vitro and in patient cohorts. The only limitation, which cannot easily be addressed but it is generally an issue in the field, is our lack of certainty about the ground truth which makes it very difficult to assess the relevance of the hypotheses provided by the computational approach. The paper will be of broad interest to the computational biology community, especially in oncology and drug discovery.

   We are researchers in the same field, we have tested several other approaches to perform similar tasks.

   Read the original source
5. 

   Version published to 10.1101/2024.07.15.603538 on bioRxiv

   Jul 17, 2024