diff --git a/scripts/plotting/chimerax_plot.py b/scripts/plotting/chimerax_plot.py index b63ba41..193d2d9 100644 --- a/scripts/plotting/chimerax_plot.py +++ b/scripts/plotting/chimerax_plot.py @@ -2,7 +2,8 @@ Use ChimeraX to generate 3D structures colored by metrics """ -import os +import os +import shlex import subprocess import math @@ -11,12 +12,41 @@ import daiquiri from scripts import __logger_name__ -from scripts.plotting.utils import detect_af_version +from scripts.plotting.utils import cap_inf_scores, detect_af_version logger = daiquiri.getLogger(__logger_name__ + ".plotting.chimerax_plot") -def create_attribute_file(path_to_file, +def _run_chimerax(command, label): + """Run a ChimeraX subprocess. Surface failures as warnings so one bad image + doesn't abort the whole loop, but the user actually sees something went wrong. + + `command` is an argv list (no shell interpretation), so paths and labels + with spaces or shell metacharacters pass through verbatim. + """ + # Discard stdout (never read) but keep stderr for the failure warning; + # capture_output=True would buffer both into memory across many genes. + try: + result = subprocess.run( + command, + stdout=subprocess.DEVNULL, + stderr=subprocess.PIPE, + text=True, + ) + except OSError as exc: + # E.g. FileNotFoundError when chimerax_bin is missing or not executable. + # Treat it like a failed run so the loop continues across other genes. + logger.warning(f"ChimeraX could not be executed for {label}: {exc}") + return None + if result.returncode != 0: + logger.warning( + f"ChimeraX failed for {label} " + f"(rc={result.returncode}): {result.stderr.strip()}" + ) + return result + + +def create_attribute_file(path_to_file, df, attribute_col, pos_col="Pos", @@ -100,11 +130,10 @@ def get_chimerax_command(chimerax_bin, transparent_bg=False): palette = get_palette(intervals, type="diverging") if attribute == "logscore" else get_palette(intervals, type="sequential") - transparent_bg = " transparentBackground true" if transparent_bg else "" - - chimerax_command = ( - f"{chimerax_bin} --nogui --offscreen --silent --cmd " - f"\"open {pdb_path}; " + transparent_bg_suffix = " transparentBackground true" if transparent_bg else "" + + chimerax_script = ( + f"open {pdb_path}; " "set bgColor white; " "color lightgray; " f"open {attr_file_path}; " @@ -118,18 +147,20 @@ def get_chimerax_command(chimerax_bin, "graphics silhouettes true width 1.3;" "zoom;" ) - + if clusters is not None and len(clusters) > 0: for pos in clusters: - chimerax_command += f"marker #10 position :{pos} color #dacae961 radius 5.919;" + chimerax_script += f"marker #10 position :{pos} color #dacae961 radius 5.919;" cluster_tag = "_clusters" else: cluster_tag = "" - + output_path = os.path.join(chimera_output_path, f"{cohort}_{i}_{gene}_{attribute}{cluster_tag}.png") - chimerax_command += f"save {output_path} pixelSize {pixelsize} supersample 3{transparent_bg};exit\"" - - return chimerax_command + chimerax_script += f"save {output_path} pixelSize {pixelsize} supersample 3{transparent_bg_suffix};exit" + + # Return as an argv list so the caller can use subprocess.run(..., shell=False) + # and avoid shell interpretation of paths/labels with spaces or metacharacters. + return [chimerax_bin, "--nogui", "--offscreen", "--silent", "--cmd", chimerax_script] def generate_chimerax_plot(output_dir, @@ -151,6 +182,7 @@ def generate_chimerax_plot(output_dir, result = pd.read_csv(pos_result_path) if "Ratio_obs_sim" in result.columns: result = result.rename(columns={"Ratio_obs_sim" : "Score_obs_sim"}) + result = cap_inf_scores(result) result["Logscore_obs_sim"] = np.log(result["Score_obs_sim"]) # Detect the AlphaFold version actually present in the dataset; the @@ -233,9 +265,9 @@ def generate_chimerax_plot(output_dir, cohort, pixelsize=pixel_size, transparent_bg=transparent_bg) - subprocess.run(chimerax_command, shell=True) - logger.debug(chimerax_command) - + _run_chimerax(chimerax_command, f"{gene} ({attribute})") + logger.debug(shlex.join(chimerax_command)) + if attribute == "score" or attribute == "logscore": chimerax_command = get_chimerax_command(chimerax_bin, pdb_path, @@ -252,8 +284,8 @@ def generate_chimerax_plot(output_dir, clusters=clusters, pixelsize=pixel_size, transparent_bg=transparent_bg) - subprocess.run(chimerax_command, shell=True) - logger.debug(chimerax_command) + _run_chimerax(chimerax_command, f"{gene} ({attribute}, clusters)") + logger.debug(shlex.join(chimerax_command)) else: tried_files = ', '.join(os.path.basename(path) for path in pdb_candidates) diff --git a/scripts/plotting/plot.py b/scripts/plotting/plot.py index 4e809da..6face81 100644 --- a/scripts/plotting/plot.py +++ b/scripts/plotting/plot.py @@ -17,6 +17,7 @@ import daiquiri from scripts.plotting.utils import ( + cap_inf_scores, get_broad_consequence, save_annotated_result, get_enriched_result, @@ -2308,6 +2309,7 @@ def generate_plots(gene_result_path, gene_result = pd.read_csv(gene_result_path) pos_result = pd.read_csv(pos_result_path) + pos_result = cap_inf_scores(pos_result) maf = pd.read_csv(maf_path, sep="\t") miss_prob_dict = json.load(open(miss_prob_path)) seq_df = pd.read_csv(seq_df_path, sep="\t") diff --git a/scripts/plotting/utils.py b/scripts/plotting/utils.py index 35bca05..4951e71 100644 --- a/scripts/plotting/utils.py +++ b/scripts/plotting/utils.py @@ -17,6 +17,50 @@ logging.getLogger('urllib3.connectionpool').setLevel(logging.WARNING) +def cap_inf_scores(pos_result, col="Score_obs_sim"): + """ + Replace +inf values in a score column with a saturated finite value. + + `Score_obs_sim` is +inf for extreme hotspot positions where the observed + anomaly score is at the mathematical limit (the high-precision Decimal + fallback in score_and_simulations.get_dcm_anomaly_score returns +inf when + even 600-digit precision underflows). +inf breaks downstream plotting: + normalization (sum=inf), log transforms, and axis auto-scaling. + + We cap +inf at 1.5 * max(finite) so the point stays visible at the top of + the y-axis without distorting the rest of the track. + """ + if col not in pos_result.columns: + return pos_result + + # Only +inf is expected here (the score is bounded below by 0); use the + # positive-only check so any future -inf sentinel is preserved instead of + # silently rewritten. + inf_mask = np.isposinf(pos_result[col]) + if not inf_mask.any(): + return pos_result + + finite_values = pos_result.loc[~inf_mask, col] + finite_max = finite_values.max() if not finite_values.empty else np.nan + if pd.notna(finite_max) and finite_max > 0: + cap = finite_max * 1.5 + # Guard against finite_max being so large that *1.5 overflows back to inf + # (would defeat the purpose of capping). Fall back to finite_max itself. + if not np.isfinite(cap): + cap = finite_max + else: + cap = 1.0 + + n_inf = int(inf_mask.sum()) + logger.warning( + f"Capping {n_inf} `{col}` value(s) at {cap:.4f} for plotting " + f"(originally +inf — extreme hotspot positions)" + ) + pos_result = pos_result.copy() + pos_result.loc[inf_mask, col] = cap + return pos_result + + _AF_VERSION_RE = re.compile(r"-model_v(\d+)\.pdb(?:\.gz)?$")