#!/usr/bin/env python3 """ER RQmodified calculation logic (self-contained, no external project dependencies). This module is a pure calculation copy synced from the iHARP project root `ER_RQmodified.py`, for standalone execution by the skill. It keeps the algorithm identical to the project core engine: - single-endpoint score ScoreSingle = 100 * reference_AC50 / chemical_AC50 - CSraw_agonist / CSraw_antagonist = harmonic mean(valid ScoreSingle) - ch1/ch2, viability are not scored directly - viability false positive: endpoint_AC50 >= paired_viability_AC50 / 3 -> excluded - LUC_VM7 special QC - concentration conversion ng/L = ng/mL * 1000 Input columns: CompoundName, Endpoint, AC50, HitCall, Concentration(ng/mL)_adjusted References: agonist 17beta-Estradiol; antagonist 4-Hydroxytamoxifen Default PNEC: E2 = 2.0, 4OHT = 30.0 Note: if the ER algorithm of the project core engine is modified, this file must be updated in sync. """ from __future__ import annotations import csv from collections import defaultdict from pathlib import Path from typing import Iterable DEFAULT_PNEC_E2 = 2.0 DEFAULT_PNEC_4OHT = 30.0 INPUT_CONCENTRATION_COLUMN = 'Concentration(ng/mL)_adjusted' OUTPUT_CONCENTRATION_COLUMN = 'Concentration(ng/L)_adjusted' AGONIST_ENDPOINTS = { 'OT_ER_ERaERa_0480', 'OT_ER_ERaERa_1440', 'OT_ER_ERaERb_0480', 'OT_ER_ERaERb_1440', 'OT_ER_ERbERb_0480', 'OT_ER_ERbERb_1440', 'OT_ERa_EREGFP_0120', 'OT_ERa_EREGFP_0480', 'ATG_chERa_XSP1_up', 'ATG_chERa_XSP2_up', 'ATG_ERa_TRANS_up', 'ATG_ERb_TRANS2_up', 'ATG_frER1_XSP1_up', 'ATG_frER1_XSP2_up', 'ATG_frER2_XSP1_up', 'ATG_frER2_XSP2_up', 'ATG_hERa_XSP1_up', 'ATG_hERa_XSP2_up', 'ATG_hERb_XSP1_up', 'ATG_hERb_XSP2_up', 'TOX21_ERa_BLA_Agonist_ratio', 'TOX21_ERa_LUC_VM7_Agonist', 'TOX21_ERa_LUC_VM7_Agonist_10nM_ICI182780', 'TOX21_ERb_BLA_Agonist_ratio', 'ACEA_ER_80hr', } ANTAGONIST_ENDPOINTS = { 'TOX21_ERa_BLA_Antagonist_ratio', 'TOX21_ERa_LUC_VM7_Antagonist_0.5nM_E2', 'TOX21_ERb_BLA_Antagonist_ratio', } # ER rule: ch1/ch2 and viability endpoints only participate in QC or retain details; they are not scored directly. NON_SCORING_ENDPOINTS = { 'TOX21_ERa_BLA_Agonist_ch1', 'TOX21_ERa_BLA_Agonist_ch2', 'TOX21_ERa_BLA_Antagonist_ch1', 'TOX21_ERa_BLA_Antagonist_ch2', 'TOX21_ERa_BLA_Antagonist_viability', 'TOX21_ERa_LUC_VM7_Agonist_10nM_ICI182780_viability', 'TOX21_ERa_LUC_VM7_Antagonist_0.1nM_E2_viability', 'TOX21_ERa_LUC_VM7_Antagonist_0.5nM_E2_viability', 'TOX21_ERb_BLA_Agonist_ch1', 'TOX21_ERb_BLA_Agonist_ch2', 'TOX21_ERb_BLA_Agonist_viability', 'TOX21_ERb_BLA_Antagonist_ch1', 'TOX21_ERb_BLA_Antagonist_ch2', 'TOX21_ERb_BLA_Antagonist_viability', } FALSE_POSITIVE_QC_PAIRS = { 'TOX21_ERa_BLA_Antagonist_ratio': 'TOX21_ERa_BLA_Antagonist_viability', 'TOX21_ERb_BLA_Agonist_ratio': 'TOX21_ERb_BLA_Agonist_viability', 'TOX21_ERb_BLA_Antagonist_ratio': 'TOX21_ERb_BLA_Antagonist_viability', 'TOX21_ERa_LUC_VM7_Agonist_10nM_ICI182780': 'TOX21_ERa_LUC_VM7_Agonist_10nM_ICI182780_viability', 'TOX21_ERa_LUC_VM7_Antagonist_0.5nM_E2': 'TOX21_ERa_LUC_VM7_Antagonist_0.5nM_E2_viability', } LUC_VM7_AGONIST_ENDPOINT = 'TOX21_ERa_LUC_VM7_Agonist' LUC_VM7_0_1NM_QC_ENDPOINT = 'TOX21_ERa_LUC_VM7_Antagonist_0.1nM_E2' LUC_VM7_0_1NM_QC_VIABILITY_ENDPOINT = 'TOX21_ERa_LUC_VM7_Antagonist_0.1nM_E2_viability' ENDPOINT_GROUPS = {endpoint: 'agonist' for endpoint in AGONIST_ENDPOINTS} ENDPOINT_GROUPS.update({endpoint: 'antagonist' for endpoint in ANTAGONIST_ENDPOINTS}) REFERENCE_BY_GROUP = { 'agonist': '17beta-Estradiol', 'antagonist': '4-Hydroxytamoxifen', } # Scoring endpoints with no reference-compound AC50 data: the reference compound # 17beta-Estradiol was not assayed for these endpoints in ToxCast, so # ScoreSingle (=100*ref_AC50/chem_AC50) can never be computed - a confirmed data gap. # These endpoints are silently skipped (no Missing reference AC50 warning) to avoid # uninformative warning spam. KNOWN_MISSING_REFERENCE_ENDPOINTS = { 'ATG_chERa_XSP1_up', 'ATG_chERa_XSP2_up', 'ATG_frER1_XSP1_up', 'ATG_frER1_XSP2_up', 'ATG_frER2_XSP1_up', 'ATG_frER2_XSP2_up', 'ATG_hERa_XSP1_up', 'ATG_hERa_XSP2_up', 'ATG_hERb_XSP1_up', 'ATG_hERb_XSP2_up', 'ATG_ERb_TRANS2_up', } DETAIL_COLUMNS = [ 'CompoundName', OUTPUT_CONCENTRATION_COLUMN, 'CSraw_agonist', 'CSraw_antagonist', 'EC_agonist', 'EC_antagonist', 'RQmodified_ag', 'RQmodified_ant', ] SUMMARY_COLUMNS = ['SUM_RQmodified_ag', 'SUM_RQmodified_ant'] def read_csv(path: str | Path) -> list[dict[str, str]]: with Path(path).open(newline='', encoding='utf-8-sig') as f: return list(csv.DictReader(f)) def write_csv(path: str | Path, rows: Iterable[dict[str, object]], fieldnames: list[str]) -> None: with Path(path).open('w', newline='', encoding='utf-8') as f: writer = csv.DictWriter(f, fieldnames=fieldnames) writer.writeheader() for row in rows: writer.writerow(row) def parse_positive_float(value: object) -> float | None: if value is None: return None text = str(value).strip() if not text: return None try: parsed = float(text) except ValueError: return None if parsed <= 0: return None return parsed def convert_ng_ml_to_ng_l(value: object) -> float | None: concentration_ng_ml = parse_positive_float(value) if concentration_ng_ml is None: return None # project convention: ng/L = ng/mL * 1000 return concentration_ng_ml * 1000 def normalize_hitcall(value: object) -> str: return str(value or '').strip().lower() def harmonic_mean(values: Iterable[float]) -> float: positives = [value for value in values if value > 0] if not positives: return 0.0 return len(positives) / sum(1 / value for value in positives) def build_reference_map(reference_rows: Iterable[dict[str, str]]) -> dict[tuple[str, str], float]: reference: dict[tuple[str, str], float] = {} scoring_endpoints = set(ENDPOINT_GROUPS) for row in reference_rows: if normalize_hitcall(row.get('HitCall')) != 'active': continue compound = str(row.get('CompoundName', '')).strip() endpoint = str(row.get('Endpoint', '')).strip() if endpoint not in scoring_endpoints: continue ac50 = parse_positive_float(row.get('AC50')) if compound and endpoint and ac50 is not None: reference.setdefault((compound, endpoint), ac50) return reference def build_endpoint_ac50_map(input_rows: Iterable[dict[str, str]]) -> dict[tuple[str, str, str], float]: endpoint_ac50: dict[tuple[str, str, str], float] = {} for row in input_rows: compound = str(row.get('CompoundName', '')).strip() endpoint = str(row.get('Endpoint', '')).strip() concentration_ng_l = convert_ng_ml_to_ng_l(row.get(INPUT_CONCENTRATION_COLUMN)) concentration = '' if concentration_ng_l is None else format_number(concentration_ng_l) ac50 = parse_positive_float(row.get('AC50')) if compound and endpoint and ac50 is not None: endpoint_ac50.setdefault((compound, concentration, endpoint), ac50) return endpoint_ac50 def is_false_positive_endpoint( compound: str, concentration: str, endpoint: str, ac50_chem: float, endpoint_ac50: dict[tuple[str, str, str], float], ) -> tuple[bool, float | None]: viability_endpoint = FALSE_POSITIVE_QC_PAIRS.get(endpoint) if viability_endpoint is None: return False, None viability_ac50 = endpoint_ac50.get((compound, concentration, viability_endpoint)) if viability_ac50 is None: return False, None return ac50_chem >= viability_ac50 / 3, viability_ac50 def should_skip_luc_vm7_agonist( compound: str, concentration: str, agonist_ac50: float, endpoint_ac50: dict[tuple[str, str, str], float], ) -> tuple[bool, str | None]: qc_ac50 = endpoint_ac50.get((compound, concentration, LUC_VM7_0_1NM_QC_ENDPOINT)) if qc_ac50 is None: return False, None viability_ac50 = endpoint_ac50.get((compound, concentration, LUC_VM7_0_1NM_QC_VIABILITY_ENDPOINT)) if viability_ac50 is not None and qc_ac50 >= viability_ac50 / 3: return False, ( 'QC endpoint false positive; agonist retained: ' f'CompoundName={compound}; {OUTPUT_CONCENTRATION_COLUMN}={concentration}; ' f'AgonistEndpoint={LUC_VM7_AGONIST_ENDPOINT}; ' f'QCEndpoint={LUC_VM7_0_1NM_QC_ENDPOINT}; QCAC50={format_number(qc_ac50)}; ' f'ViabilityEndpoint={LUC_VM7_0_1NM_QC_VIABILITY_ENDPOINT}; ' f'ViabilityAC50={format_number(viability_ac50)}; Rule=QC AC50 >= viability AC50 / 3' ) if qc_ac50 >= 10 * agonist_ac50: return True, ( 'False positive agonist endpoint skipped: ' f'CompoundName={compound}; {OUTPUT_CONCENTRATION_COLUMN}={concentration}; ' f'Endpoint={LUC_VM7_AGONIST_ENDPOINT}; AC50={format_number(agonist_ac50)}; ' f'QCEndpoint={LUC_VM7_0_1NM_QC_ENDPOINT}; QCAC50={format_number(qc_ac50)}; ' 'Rule=0.1nM_E2 AC50 >= 10 * agonist AC50' ) return False, None def format_number(value: float) -> str: if abs(value) < 1e-15: value = 0.0 return f'{value:.12g}' def max_rqmodified(rqmodified_ag: object, rqmodified_ant: object) -> float | None: values = [] for value in (rqmodified_ag, rqmodified_ant): try: if value is not None: values.append(float(value)) except (TypeError, ValueError): pass if not values: return None return max(values) def calculate_rqmodified( input_rows: Iterable[dict[str, str]], reference_rows: Iterable[dict[str, str]], pnec_e2: float = DEFAULT_PNEC_E2, pnec_4oht: float = DEFAULT_PNEC_4OHT, ) -> tuple[list[dict[str, object]], list[str], dict[str, float]]: input_rows = list(input_rows) reference = build_reference_map(reference_rows) endpoint_ac50 = build_endpoint_ac50_map(input_rows) warnings: list[str] = [] grouped_scores: dict[tuple[str, str], dict[str, list[float]]] = defaultdict( lambda: {'agonist': [], 'antagonist': []} ) seen_compound_concentrations: set[tuple[str, str]] = set() concentrations_by_compound: dict[str, set[str]] = defaultdict(set) missing_reference_seen: set[tuple[str, str, str]] = set() for row in input_rows: compound = str(row.get('CompoundName', '')).strip() concentration_ng_l = convert_ng_ml_to_ng_l(row.get(INPUT_CONCENTRATION_COLUMN)) concentration = '' if concentration_ng_l is None else format_number(concentration_ng_l) if not compound: continue seen_compound_concentrations.add((compound, concentration)) concentrations_by_compound[compound].add(concentration) for compound in sorted(concentrations_by_compound): concentrations = sorted(concentrations_by_compound[compound]) if len(concentrations) > 1: warnings.append( 'Concentration not unique: ' f'CompoundName={compound}; {OUTPUT_CONCENTRATION_COLUMN}={"|".join(concentrations)}' ) for row in input_rows: if normalize_hitcall(row.get('HitCall')) != 'active': continue compound = str(row.get('CompoundName', '')).strip() concentration_ng_l = convert_ng_ml_to_ng_l(row.get(INPUT_CONCENTRATION_COLUMN)) concentration = '' if concentration_ng_l is None else format_number(concentration_ng_l) endpoint = str(row.get('Endpoint', '')).strip() if not compound: continue if endpoint in NON_SCORING_ENDPOINTS: continue group = ENDPOINT_GROUPS.get(endpoint) if group is None: continue ac50_chem = parse_positive_float(row.get('AC50')) if ac50_chem is None: warnings.append( f'Invalid AC50 skipped: CompoundName={compound}; ' f'{INPUT_CONCENTRATION_COLUMN}={row.get(INPUT_CONCENTRATION_COLUMN)}; Endpoint={endpoint}; AC50={row.get("AC50")}' ) continue if endpoint == LUC_VM7_AGONIST_ENDPOINT: skip_agonist, qc_warning = should_skip_luc_vm7_agonist( compound, concentration, ac50_chem, endpoint_ac50, ) if qc_warning is not None: warnings.append(qc_warning) if skip_agonist: continue is_false_positive, viability_ac50 = is_false_positive_endpoint( compound, concentration, endpoint, ac50_chem, endpoint_ac50, ) if is_false_positive: warnings.append( 'False positive endpoint skipped: ' f'CompoundName={compound}; {OUTPUT_CONCENTRATION_COLUMN}={concentration}; ' f'Endpoint={endpoint}; AC50={format_number(ac50_chem)}; ' f'ViabilityEndpoint={FALSE_POSITIVE_QC_PAIRS[endpoint]}; ' f'ViabilityAC50={format_number(viability_ac50)}; Rule=endpoint AC50 >= viability AC50 / 3' ) continue reference_compound = REFERENCE_BY_GROUP[group] reference_ac50 = reference.get((reference_compound, endpoint)) if reference_ac50 is None: if endpoint in KNOWN_MISSING_REFERENCE_ENDPOINTS: continue # known no-reference-data endpoint, skip silently without warning key = (endpoint, group, reference_compound) if key not in missing_reference_seen: warnings.append( 'Missing reference AC50: ' f'Endpoint={endpoint}; Group={group}; Reference={reference_compound}' ) missing_reference_seen.add(key) continue score_single = 100 * reference_ac50 / ac50_chem if score_single > 0: grouped_scores[(compound, concentration)][group].append(score_single) for key in seen_compound_concentrations: grouped_scores.setdefault(key, {'agonist': [], 'antagonist': []}) results: list[dict[str, object]] = [] for compound, concentration in sorted(grouped_scores, key=lambda item: (item[0], item[1])): scores = grouped_scores[(compound, concentration)] csraw_agonist = harmonic_mean(scores['agonist']) csraw_antagonist = harmonic_mean(scores['antagonist']) concentration_value = parse_positive_float(concentration) ec_agonist: float | None = None ec_antagonist: float | None = None rqmodified_ag: float | None = None rqmodified_ant: float | None = None if concentration_value is not None: ec_agonist = concentration_value * (csraw_agonist / 100) ec_antagonist = concentration_value * (csraw_antagonist / 100) if pnec_e2 > 0 and csraw_agonist > 0: rqmodified_ag = ec_agonist / pnec_e2 if pnec_4oht > 0 and csraw_antagonist > 0: rqmodified_ant = ec_antagonist / pnec_4oht results.append( { 'CompoundName': compound, OUTPUT_CONCENTRATION_COLUMN: concentration, 'CSraw_agonist': csraw_agonist, 'CSraw_antagonist': csraw_antagonist, 'EC_agonist': ec_agonist if csraw_agonist > 0 else None, 'EC_antagonist': ec_antagonist if csraw_antagonist > 0 else None, 'RQmodified_ag': rqmodified_ag, 'RQmodified_ant': rqmodified_ant, } ) totals = { 'SUM_RQmodified_ag': sum(float(row['RQmodified_ag']) for row in results if row['RQmodified_ag'] is not None), 'SUM_RQmodified_ant': sum(float(row['RQmodified_ant']) for row in results if row['RQmodified_ant'] is not None), } return results, warnings, totals def format_detail_rows(rows: Iterable[dict[str, object]]) -> list[dict[str, str]]: formatted_rows: list[dict[str, str]] = [] for row in rows: formatted_rows.append( { 'CompoundName': str(row['CompoundName']), OUTPUT_CONCENTRATION_COLUMN: str(row[OUTPUT_CONCENTRATION_COLUMN]), 'CSraw_agonist': format_number(float(row['CSraw_agonist'])), 'CSraw_antagonist': format_number(float(row['CSraw_antagonist'])), 'EC_agonist': '' if row['EC_agonist'] is None else format_number(float(row['EC_agonist'])), 'EC_antagonist': '' if row['EC_antagonist'] is None else format_number(float(row['EC_antagonist'])), 'RQmodified_ag': '' if row['RQmodified_ag'] is None else format_number(float(row['RQmodified_ag'])), 'RQmodified_ant': '' if row['RQmodified_ant'] is None else format_number(float(row['RQmodified_ant'])), } ) return formatted_rows