| |
| """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', |
| } |
|
|
| |
| 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', |
| } |
|
|
| |
| |
| |
| |
| |
| 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 |
| |
| 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 |
| 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 |
|
|
