{"id":"7aaf3dce-a48b-44cf-be39-6d490396857c","arxiv_id":"2607.13994","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"Causal graphs learned from 108 low-dose-radiation RNA-seq samples yield gene sets more enriched for radiation-response pathways than differential-expression baselines and show transcription-factor hubs and housekeeping sinks.","lead":"This study applies causal-discovery algorithms to RNA-sequencing data from human cells exposed to low-dose radiation, reporting that the inferred gene networks pick out known radiation-response genes better than standard differential-expression analysis. It is a use case for mining mechanism hypotheses from small, perturbed transcriptomics cohorts, with several statistical weaknesses that the authors partly acknowledge.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"DESeq2-on-LFC prefilter is statistically invalid and seeds the gene universe for every downstream result; the central claim is untested until it is rerun on counts.","rationale":"The central claim is that causal discovery with a joint radiation variable identifies radiation-response gene sets better than DE or supervised ML. For that to be true, the gene universe fed to DAG-GNN must at least be a defensible candidate set. The paper does not establish this: DESeq2 is applied to an LFC matrix, violating the method's count-data assumptions. This is not a disagreement with consensus; it is a statistical misapplication. The concern is testable by re-running the prefilter on counts. The reader's weakest assumption identifies the same point; I agree. A secondary issue — small n, 10 bootstraps, and STRING-based partition/merge — could also threaten edge-level claims, but the prefilter is more fundamental because it precedes those steps. The reader's CONDITIONAL verdict remains appropriate: the findings are biologically plausible, the authors explicitly label the analysis as primarily qualitative, and the issue is addressable. The paper as written cannot be accepted until the prefilter is corrected and the main comparisons are recomputed on a valid gene pool.","tokens_in":13398,"tokens_out":5577,"duration_ms":64203,"concrete_test":"Require the raw per-sample count matrix (from the authors or the Jantre et al. data release) and rerun the §4.1.1 prefilter correctly: pyDESeq2 on counts with a design including dose rate and week, Wald test, BH FDR; define the per-dose union with FDR<0.05 and |log2FC|>1. Compare this corrected pool with the p≤0.05 pool used in the paper, in both overlap and pathway composition. Then rerun DAG-GNN and the enrichment analysis on the corrected pool, or at least recompute Fig. 4 and Table 3 with corrected DE gene sets. If the corrected causal invariant gene set no longer shows the reported radiation-pathway enrichment, or if the causal-vs-DE advantage disappears, the central claim fails. If corrected and original pools overlap by >90% and the enrichment top hits are essentially unchanged, the concern is settled.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 4.1.1 states that differential expression was \"performed using pyDESeq2 v0.4.9 with dataset X_LFC\" — a continuous matrix of log2 fold changes, not counts. DESeq2's negative-binomial model and Wald test are defined for count data; applied to an LFC matrix, the p-values have no defined statistical meaning. These p-values are the gate for the entire causal pipeline: §4.1.3 filters out genes with p>0.05 before DAG-GNN, yielding per-dose gene pools of 6,274–10,603 genes. Every downstream claim — the causal gene sets, the invariant 438-gene set, the enrichment comparisons in Fig. 4/Table 3, the TF/housekeeping hub analysis, and the stress-versus-apoptosis branch — is computed inside this invalidly chosen gene universe. The comparison to \"differential expression\" baselines is doubly affected because those baselines use the same invalid p-values. The random-gene control in Appendix D does not fix this: random sets are drawn from the full background, not from the potentially arbitrary filtered pool. This is the load-bearing entry point of the analysis.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper applies causal discovery (DAG-GNN) to bulk RNA-seq data from RPE1 cells exposed to five radiation dose rates, jointly modeling a cumulative radiation-dose variable with gene expression. It compares the resulting causal gene sets against differential expression (DESeq2) and random-forest baselines via pathway enrichment, and analyzes the learned directed graphs for known regulatory structure (transcription-factor hubs, housekeeping sinks, invariant subgraphs). The central claim is that causal discovery achieves a higher enrichment for radiation-response pathways and provides structured hypotheses about stress-versus-apoptosis branching. The paper includes code, a real dataset, and detailed appendices with additional baselines and bootstrap analyses.","tokens_in":13647,"tokens_out":3388,"duration_ms":38527,"significance":"If the methodological concerns were resolved, this would be a valuable demonstration of causal discovery on low-sample, perturbed transcriptomics—a setting where exact graph recovery is known to be hard but where structured hypothesis generation could still be useful. The graph-structural findings (TF hubs, housekeeping sinks, a stress/apoptosis branch) are concrete and experimentally testable, and the paper provides a public implementation. However, the statistical validity of the prefilter and the circularity in graph validation currently prevent full confidence in the conclusions.","major_comments":[{"comment":"DESeq2 is applied to X_LFC, a matrix of log2 ratios of TPM values, not raw counts. DESeq2's negative-binomial model and Wald test are defined for count data; applying them to continuous log-ratios gives p-values with no defined statistical meaning. These p-values are the gate for the entire causal pipeline (genes with p>0.05 are removed before DAG-GNN), so every downstream result—causal gene sets, invariant 438-gene set, enrichment comparisons in Fig. 4/Table 3, and structural analyses—is computed over an arbitrarily filtered gene universe. This is a load-bearing error that requires rerunning the DE analysis on raw counts or an appropriate continuous-data method and repeating the full pipeline.","section":"§4.1.1 and Table 1"},{"comment":"The STRING PPI network is used as the initial undirected structure to define the causal partition (Appendix A), constraining the search space for DAG-GNN. The same STRING network is then used as a validation target for edge overlap in Table 2. This creates partial circularity: edges consistent with STRING are more likely to be discovered because the partition restricts possible edges. The reported F1 values (max 0.031) are low, but they are not calibrated against the constrained search space. The random-gene control in Appendix D draws from the full background, not from the STRING-partitioned pool, so it does not correct this bias.","section":"§4.1.3, Appendix A, and Table 2"},{"comment":"With n=36 samples per dose and up to 2,866 variables per DAG-GNN run, the problem is severely underdetermined. The consensus graph is built from only 10 bootstrap resamples with an arbitrary 50% edge-presence threshold, and the adjacency threshold is fixed at the DAG-GNN default 0.3. Appendix G shows very high variance across individual bootstrap runs; e.g., edge overlap with knowledge bases varies widely. The authors note the analysis is 'primarily qualitative,' but the central claim of 'higher hit rate' rests on consensus graphs derived from these fragile thresholds. A sensitivity analysis (e.g., varying the consensus threshold or adjacency threshold) or a more robust statistical evaluation is needed.","section":"§4.1.3 and Appendix G"},{"comment":"The comparison of pathway enrichment across methods is not statistically formalized. Fig. 4 reports -log10 p-values without confidence intervals or a formal test for whether causal gene sets are significantly more enriched than DE or random-forest sets. The causal gene sets are generally smaller than DE sets, and the figures do not control for gene-set size or the DE-p-value prefilter. The random-gene control in Appendix D (drawing from the full background) is informative but does not address the specific bias induced by the STRING partition or the DE prefilter. A permutation test that samples gene sets of equal size from the same filtered pool would strengthen the claim.","section":"Fig. 4 and Table 3"}],"minor_comments":[{"comment":"The training-curve figure is referenced as 'Fig.??'—the placeholder remains in the manuscript.","section":"Appendix A"},{"comment":"The caption says 'dose rate I (6.66 mGy/min)' but the unit should be mGy/hr, consistent with the rest of the paper.","section":"Table 6"},{"comment":"The random-forest feature-selection threshold (20 of 25 folds) and the correlation baseline threshold (>50% folds) are arbitrary; please justify or perform sensitivity checks.","section":"§4.1.2"},{"comment":"The caption mentions '5 repeated runs' but the text does not specify how many random gene sets were drawn; clarify the procedure.","section":"Fig. 15"},{"comment":"The text states 'Differential expression identifies larger sets' but does not give the exact sizes of the causal gene sets at each dose rate; reporting these sizes (and the number of nodes in each consensus graph) would help interpret Venn diagrams and enrichment comparisons.","section":"§5"}],"recommendation":"major_revision","confidential_remarks":"The DESeq2-on-LFC issue is indeed the load-bearing weakness identified by the stress-test. It is not merely a presentation issue—it invalidates the prefilter and thus the entire gene universe for the causal analysis. The STRING circularity is also real, though its impact on the main enrichment claim is indirect (it affects which nodes have edges and thus which genes are selected). The authors should be given the opportunity to rerun the pipeline on raw counts and to recalibrate the STRING validation. I expect that after these fixes the paper could be suitable for publication, but in its current form the central claim is not supported."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"You should know two things about this paper. First, it is not a routine X-applied-to-Y. Jointly modeling a continuous dose-rate perturbation with gene expression via DAG-GNN on bulk RNA-seq, then studying hub/sink structure and an invariant causal gene set, is a real methodological step forward for low-sample perturbation studies. Second, the entry point of the analysis is statistically broken, and the break propagates into almost everything they report.\n\nThe good parts: the idea is sensible, the paper is plainly written, and the authors include controls that many applied papers skip — random-gene enrichment checks, a correlation baseline, and honest discussion that the analysis is qualitative. The structural findings (TF hubs, housekeeping sinks, stress-versus-apoptosis branch) are biologically credible and worth testing. The scaling trick from Shah et al. is interesting and appears to work.\n\nThe soft spot is not minor. Section 4.1.1 says DESeq2 was run on X_LFC, a matrix of log2 fold changes, not counts. DESeq2's negative-binomial model is defined for raw counts; feeding it ratio-transformed TPMs produces p-values with no defined meaning. Those p-values are the gate for the entire causal pipeline: genes with p>0.05 are discarded before DAG-GNN, producing the dose-specific gene pools of 6,000–10,000 genes. Every downstream claim — the causal gene sets, the invariant 438 genes, the enrichment comparisons, the hub/sink analysis — lives inside that arbitrarily chosen universe. The comparison to differential expression baselines is doubly confounded because those baselines are computed from the same invalid p-values. The random-gene control does not fix it: random sets are drawn from the full background, not from the filtered pool. This is load-bearing, and the stress-test note is right.\n\nA secondary but real issue is the STRING circularity: STRING is used as the skeleton for the causal partition and then as a validation target for edge overlap. That inflates the reported overlap at least somewhat, and the edge validation numbers are anyway close to zero (F1 ≤ 0.03). The bootstrap consensus with only 10 resamples and a 50% threshold is thin, and the authors themselves flag the lack of error bars.\n\nWho is this for? Anyone working on causal discovery for healthcare transcriptomics, or on gene prioritization in low-sample perturbation studies. They will see a useful framework but will need to treat the results as hypotheses. The paper deserves a serious referee — the ideas are worth the effort of a careful review and the fix is straightforward: rerun the DE prefilter on count data, redo the gene set comparisons with matched sizes and bootstrap intervals, and separate STRING-as-skeleton from STRING-as-validation. If the enrichment advantage survives that, it becomes a solid paper. My recommendation to an editor: send it to review with a request for major revision, not a desk reject.","headline":"A genuinely novel application of causal discovery to perturbed transcriptomics, undermined by a statistically invalid DESeq2-on-LFC gene filter that seeds every downstream result; the paper deserves a serious referee but the central claim is untested until that filter is fixed.","tokens_in":14195,"tokens_out":1443,"would_cite":false,"duration_ms":15978,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Jointly modeling radiation exposure with gene expression, the authors aim to show that causal discovery finds perturbation-response gene sets and regulatory structure that differential expression misses.","keywords":["causal discovery","gene regulatory network","DAG-GNN","radiation response","RNA-seq","pathway enrichment","dose-rate perturbation","transcription factor hubs"],"falsifier":"Re-run the entire pipeline with the gene prefilter performed by DESeq2 (or an equivalent count-based method) on the raw count matrix rather than on log2-fold-change values, and check whether the causal gene sets still show stronger radiation-pathway enrichment than the differential-expression or random-forest baselines; if the 438-gene invariant set shifts and the enrichment advantage disappears, the claim collapses.","tokens_in":13182,"feed_emoji":"🧬","tokens_out":5702,"duration_ms":54788,"temperature":0.7,"pith_summary":"This paper tries to establish that causal discovery—not just correlation or per-gene fold changes—can identify the gene programs behind radiation response in human cells. By feeding a graph-learning model both the measured gene expression and the cumulative radiation dose, the authors learn directed gene networks for each dose rate. They report that these networks flag smaller gene sets that are more strongly enriched for known radiation-response pathways than gene sets from differential expression, random forests, or linear correlation; that the network hubs are transcription factors and the sinks are housekeeping genes; and that an invariant 438-gene core across dose rates encodes the expected radiation-response pathways. If true, this would make causal discovery a practical complement to standard bioinformatics for turning small perturbation studies into structural mechanism hypotheses.","feed_headline":"Causal graphs find radiation-response genes that DE misses","feed_subtitle":"Joint modeling of dose and expression yields smaller gene sets with stronger pathway enrichment, plus testable network hypotheses.","key_machinery":"The central mechanism is DAG-GNN, a variational autoencoder that learns a directed adjacency matrix under a differentiable acyclicity constraint, here trained to reconstruct the joint distribution of gene expression and cumulative radiation dose. Because the full 15,694-gene matrix is intractable, the authors pre-filter genes by differential-expression p-values, partition the remaining genes into overlapping subsets using the STRING protein interaction network as an undirected skeleton, learn a DAG per subset, and merge the subgraphs. Ten bootstrap resamples produce a consensus graph (edges present in at least 50% of runs) at each dose rate; 'perfect' edges present in all ten runs are treate","core_discovery":"On its own terms, the paper claims that causal discovery over the joint distribution of radiation dose and gene expression outperforms standard gene-selection pipelines for perturbation response. Using RPE1 cells exposed to five dose rates over nine weeks, the authors build dose-specific directed acyclic graphs with a graph variational autoencoder, include the cumulative radiation dose as an explicit variable, and retain edges that appear in at least half of ten bootstrap resamples. The resulting gene sets are smaller than those from differential expression yet are more strongly enriched for curated radiation-response pathways, and the invariant 438-gene intersection is enriched almost entir","pith_inferences":["If the approach holds, the same joint-perturbation causal-discovery recipe could transfer to drug-response or disease-progression transcriptomics, where structured perturbations and few samples are the norm.","Because the graph is built over a p-value-filtered gene pool, the claim's scope is restricted to genes passing that prefilter; the method's generalizability depends on the prefilter being statistically valid.","The invariant-gene intersection behaves like an environment-invariant causal signature, so the dose-rate-stable edges might be testable mechanistic predictions: silencing the upstream transcription-factor hubs should shift the downstream housekeeping-ancestor enrichment profile.","A cheap additional robustness check the authors did not run would be to compare the consensus graphs against graphs learned from shuffled radiation-dose labels, to ensure the perturbation variable is doing the work and not the expression covariance alone."],"forward_implications":["A joint model of radiation exposure and multivariate gene expression produces gene sets with stronger radiation-specific pathway enrichment than differential expression, random forest feature importance, or a linear correlation baseline, even though these sets are smaller.","The learned directed graphs consistently place transcription factors at high out-degree hubs and housekeeping genes at high in-degree sinks, and this assignment is not driven by expression variance (maximum Spearman rho = 0.119).","The 438 genes present in every dose-rate graph are enriched almost exclusively for known radiation-response pathways (cell cycle, p53 signaling, DNA damage response, senescence), supporting the intersection as a core response set.","Ancestors of housekeeping genes are enriched for oxidative-stress response while non-ancestors are enriched for membrane-permeability and apoptosis processes, suggesting a stress-versus-death branch in the graph.","Edges that appear in 100% of bootstrap runs, though mostly absent from current knowledge bases, are candidates for novel radiation-response mechanisms that could be tested by CRISPR knockout."],"fun_headline_variants":["Causal discovery beats DE for radiation gene networks","Dose-aware causal graphs reveal radiation response genes","Causal graphs uncover radiation genes DE can't see","Gene networks from causal discovery outdo differential expression","Causal modeling finds stronger radiation response signals"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"The whole pipeline depends on the statistical validity of the gene prefilter, which computes p-values with a negative-binomial test designed for count data but applies it to a log2-fold-change matrix—so if those p-values are meaningless, the causal graphs are learned over an arbitrary gene subset and the enrichment comparisons lose their footing.","fun_headline_variants_meta":{"raw":{"variants":["Causal discovery beats DE for radiation gene networks","Dose-aware causal graphs reveal radiation response genes","Causal graphs uncover radiation genes DE can't see","Gene networks from causal discovery outdo differential expression","Causal modeling finds stronger radiation response signals"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00015,"raw_usage":{"total_tokens":998,"prompt_tokens":676,"completion_tokens":322,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":420,"completion_tokens_details":{"reasoning_tokens":251}},"tokens_in":420,"tokens_out":322,"duration_ms":3201,"temperature":1.0,"reasoning_tokens":251,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-02T03:05:21.034103+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Re-run the entire pipeline with the gene prefilter performed by DESeq2 (or an equivalent count-based method) on the raw count matrix rather than on log2-fold-change values, and check whether the causal gene sets still show stronger radiation-pathway enrichment than the differential-expression or random-forest baselines; if the 438-gene invariant set shifts and the enrichment advantage disappears, the claim collapses.","supporting_citations":[],"review_version":1}