""" petgluco_core.py — PETGlucoFlux 정량 엔진 (표준 라이브러리만 사용). rodent micro-PET 18F-FDG 정량의 핵심 계산을 담는다: - 18F decay 보정 (반감기 109.77 min) - SUV (BW / lean / glucose-corrected) - Patlak graphical analysis -> net influx Ki - lumped-constant 보정 조직 MRGlu - 종단(within-animal) paired 변화 - cohort 요약 통계 (paired t-test, scipy 없이 stdlib t-분포 근사) - mechanism 분류 설계 원칙: numpy/scipy/pandas 없이 순수 python 으로 동작한다. app.py 는 이 모듈을 그대로 재사용하되, 시각화/업로드만 streamlit + matplotlib 로 감싼다. 본 도구는 연구용·참고용이며, 정량 결과는 사용자 검증이 필요하다. 원시 DICOM 재구성은 범위 밖이며, 정량 CSV ingest 를 전제로 한다. 합성(synthetic) 데이터를 기본 샘플로 제공한다. """ import csv import math # --------------------------------------------------------------------------- # 물리 상수 / 모델 기본값 # --------------------------------------------------------------------------- F18_HALF_LIFE_MIN = 109.77 # 18F 반감기 (min) F18_DECAY_LAMBDA = math.log(2) / F18_HALF_LIFE_MIN # 조직별 lumped constant (LC) 기본값 (rodent 문헌 기반 대표 가정값; 사용자 검증 필요) DEFAULT_LUMPED_CONSTANT = { "skeletal_muscle": 1.00, "myocardium": 1.00, "brain": 0.60, "BAT": 1.10, "liver": 1.00, } DEFAULT_LC_FALLBACK = 1.00 # 모델별 조직 SUV / Ki 참고 범위 (대표 가정값 — 연구 설계 sanity-check 용도) # 형식: model -> tissue -> {"suv":(lo,hi), "ki":(lo,hi)} (BW-SUV g/mL, Ki mL/min/mL) MODEL_REFERENCE = { "C57BL/6": { "skeletal_muscle": {"suv": (0.3, 1.5), "ki": (0.005, 0.030)}, "myocardium": {"suv": (3.0, 10.0), "ki": (0.05, 0.15)}, "BAT": {"suv": (1.0, 8.0), "ki": (0.02, 0.12)}, "brain": {"suv": (1.0, 3.5), "ki": (0.03, 0.09)}, "liver": {"suv": (0.3, 1.2), "ki": (0.004, 0.020)}, }, "db/db": { "skeletal_muscle": {"suv": (0.2, 1.0), "ki": (0.003, 0.020)}, "myocardium": {"suv": (2.0, 7.0), "ki": (0.03, 0.10)}, "BAT": {"suv": (0.5, 4.0), "ki": (0.01, 0.07)}, "brain": {"suv": (1.0, 3.5), "ki": (0.03, 0.09)}, "liver": {"suv": (0.3, 1.5), "ki": (0.004, 0.025)}, }, "ob/ob": { "skeletal_muscle": {"suv": (0.2, 1.0), "ki": (0.003, 0.020)}, "myocardium": {"suv": (2.0, 7.0), "ki": (0.03, 0.10)}, "BAT": {"suv": (0.5, 4.0), "ki": (0.01, 0.07)}, "brain": {"suv": (1.0, 3.5), "ki": (0.03, 0.09)}, "liver": {"suv": (0.3, 1.5), "ki": (0.004, 0.025)}, }, "DIO": { "skeletal_muscle": {"suv": (0.2, 1.2), "ki": (0.004, 0.025)}, "myocardium": {"suv": (2.5, 8.0), "ki": (0.04, 0.12)}, "BAT": {"suv": (0.5, 5.0), "ki": (0.01, 0.09)}, "brain": {"suv": (1.0, 3.5), "ki": (0.03, 0.09)}, "liver": {"suv": (0.3, 1.4), "ki": (0.004, 0.022)}, }, } # Patlak 선형 구간 기본 t* (min). 보통 10-15 min 이후가 선형. DEFAULT_TSTAR_MIN = 12.0 # --------------------------------------------------------------------------- # CSV ingest # --------------------------------------------------------------------------- def _to_float(x): try: return float(x) except (TypeError, ValueError): return None def load_csv(path): """CSV -> list[dict]. 헤더 그대로 키로 사용.""" with open(path, newline="") as f: return list(csv.DictReader(f)) def load_tac(path): rows = load_csv(path) out = [] for r in rows: out.append({ "animal_id": r["animal_id"], "group": r["group"], "timepoint": r["timepoint"], "tissue": r["tissue"], "time_min": _to_float(r["time_min"]), "activity_bqml": _to_float(r["activity_bqml"]), }) return out def load_meta(path): rows = load_csv(path) out = {} for r in rows: key = (r["animal_id"], r["timepoint"]) out[key] = { "animal_id": r["animal_id"], "group": r["group"], "timepoint": r["timepoint"], "dose_mbq": _to_float(r["dose_mbq"]), "body_weight_g": _to_float(r["body_weight_g"]), "lean_mass_g": _to_float(r["lean_mass_g"]), "glucose_mgdl": _to_float(r["glucose_mgdl"]), "scan_temp_c": _to_float(r.get("scan_temp_c")), "anesthesia": r.get("anesthesia", ""), } return out def load_input_function(path): rows = load_csv(path) out = {} for r in rows: key = (r["animal_id"], r["timepoint"]) out.setdefault(key, []).append({ "time_min": _to_float(r["time_min"]), "plasma_activity_bqml": _to_float(r["plasma_activity_bqml"]), }) for k in out: out[k].sort(key=lambda d: d["time_min"]) return out # --------------------------------------------------------------------------- # decay 보정 + QC # --------------------------------------------------------------------------- def decay_correct(activity, time_min): """raw(측정) -> 주사 시점(t=0) 기준 decay-corrected 활성도. A_corr = A_meas * exp(+lambda * t).""" return activity * math.exp(F18_DECAY_LAMBDA * time_min) def decay_correct_tac(tac_rows): """TAC 각 행에 decay_corrected_bqml 추가.""" for r in tac_rows: if r["time_min"] is None or r["activity_bqml"] is None: r["decay_corrected_bqml"] = None else: r["decay_corrected_bqml"] = decay_correct(r["activity_bqml"], r["time_min"]) return tac_rows def qc_flags(tac_rows, input_funcs, meta): """간단한 QC 플래그. - extravasation 의심: 매우 낮은 input function peak (주사 외삼출 -> 혈중 활성도 저하) - 움직임 의심: TAC 가 단조 증가가 아니라 후반 frame 에서 급격 변동 반환: list[str] 경고 메시지. """ flags = [] # input peak 기반 extravasation (군 평균 대비 50% 미만) peaks = {} for key, rows in input_funcs.items(): peaks[key] = max((r["plasma_activity_bqml"] or 0) for r in rows) if rows else 0 if peaks: mean_peak = sum(peaks.values()) / len(peaks) for key, pk in peaks.items(): if mean_peak > 0 and pk < 0.5 * mean_peak: flags.append( "QC[extravasation 의심] %s/%s: input peak %.0f Bq/mL " "(cohort 평균 %.0f 의 50%% 미만)" % (key[0], key[1], pk, mean_peak) ) # 움직임: decay-corrected TAC 후반(>=t*) 프레임 간 상대 변동 > 25% by_curve = {} for r in tac_rows: by_curve.setdefault((r["animal_id"], r["timepoint"], r["tissue"]), []).append(r) for key, rows in by_curve.items(): rows = sorted([x for x in rows if x.get("decay_corrected_bqml") is not None], key=lambda d: d["time_min"]) late = [x for x in rows if x["time_min"] >= DEFAULT_TSTAR_MIN] for i in range(1, len(late)): prev, cur = late[i - 1]["decay_corrected_bqml"], late[i]["decay_corrected_bqml"] if prev and abs(cur - prev) / prev > 0.25: flags.append( "QC[움직임/재구성 의심] %s/%s/%s: %.1f->%.1f min 사이 decay-corrected " "활성도 %.0f%% 변동" % (key[0], key[1], key[2], late[i - 1]["time_min"], late[i]["time_min"], 100 * (cur - prev) / prev) ) break return flags # --------------------------------------------------------------------------- # SUV # --------------------------------------------------------------------------- def suv(tissue_activity_bqml, dose_mbq, mass_g): """SUV = (조직 활성도 / 주사량) x 체중. 단위 정합: 조직 활성도 Bq/mL, 주사량 MBq(=1e6 Bq), 질량 g(~mL). SUV[g/mL] = activity[Bq/mL] / (dose[Bq]) * mass[g] """ if dose_mbq is None or mass_g is None or dose_mbq <= 0: return None dose_bq = dose_mbq * 1.0e6 return tissue_activity_bqml / dose_bq * mass_g def static_activity_from_tac(tac_rows, animal_id, timepoint, tissue, window_min=(40.0, 1e9)): """정적 SUV용 조직 활성도: 후반(uptake plateau) 프레임의 decay-corrected 평균.""" sel = [r for r in tac_rows if r["animal_id"] == animal_id and r["timepoint"] == timepoint and r["tissue"] == tissue and r.get("decay_corrected_bqml") is not None and window_min[0] <= r["time_min"] <= window_min[1]] if not sel: # fallback: 마지막 프레임 sel = [r for r in tac_rows if r["animal_id"] == animal_id and r["timepoint"] == timepoint and r["tissue"] == tissue and r.get("decay_corrected_bqml") is not None] if not sel: return None sel = [max(sel, key=lambda d: d["time_min"])] return sum(r["decay_corrected_bqml"] for r in sel) / len(sel) def suv_table(tac_rows, meta): """조직 x animal x timepoint 별 SUV(BW/lean/glucose-corrected) 테이블. SUVglu = SUV_bw * (glucose_mgdl / 100) (혈당 정규화; 100 mg/dL 기준).""" out = [] keys = sorted({(r["animal_id"], r["timepoint"], r["tissue"]) for r in tac_rows}) for aid, tp, tissue in keys: m = meta.get((aid, tp)) if not m: continue act = static_activity_from_tac(tac_rows, aid, tp, tissue) if act is None: continue suv_bw = suv(act, m["dose_mbq"], m["body_weight_g"]) suv_lean = suv(act, m["dose_mbq"], m["lean_mass_g"]) suv_glu = None if suv_bw is not None and m["glucose_mgdl"]: suv_glu = suv_bw * (m["glucose_mgdl"] / 100.0) out.append({ "animal_id": aid, "group": m["group"], "timepoint": tp, "tissue": tissue, "tissue_activity_bqml": act, "SUV_bw": suv_bw, "SUV_lean": suv_lean, "SUV_glu": suv_glu, }) return out # --------------------------------------------------------------------------- # Patlak graphical analysis # --------------------------------------------------------------------------- def _linreg(xs, ys): """OLS 선형회귀 -> (slope, intercept, r2).""" n = len(xs) if n < 2: return None, None, None mx = sum(xs) / n my = sum(ys) / n sxx = sum((x - mx) ** 2 for x in xs) sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) if sxx == 0: return None, None, None slope = sxy / sxx intercept = my - slope * mx ss_tot = sum((y - my) ** 2 for y in ys) ss_res = sum((y - (slope * x + intercept)) ** 2 for x, y in zip(xs, ys)) r2 = 1 - ss_res / ss_tot if ss_tot > 0 else None return slope, intercept, r2 def _interp(curve_t, curve_v, t): """선형 보간 (curve 는 time 정렬 가정).""" if t <= curve_t[0]: return curve_v[0] if t >= curve_t[-1]: return curve_v[-1] for i in range(1, len(curve_t)): if t <= curve_t[i]: t0, t1 = curve_t[i - 1], curve_t[i] v0, v1 = curve_v[i - 1], curve_v[i] if t1 == t0: return v0 return v0 + (v1 - v0) * (t - t0) / (t1 - t0) return curve_v[-1] def patlak(tissue_tac, input_func, tstar_min=None, auto_tstar=True): """Patlak graphical analysis. 입력: tissue_tac: list[(time_min, decay_corrected_activity)] (조직) input_func: list[(time_min, plasma_activity)] (decay-corrected plasma Cp) Patlak 변수: x = (integral_0^t Cp dt') / Cp(t) y = C_tissue(t) / Cp(t) 선형 구간 (t >= t*) 의 slope = Ki (net influx, mL/min/mL), intercept = 분포부피항. 반환 dict: Ki, intercept, r2, tstar, n_points, points(list of (x,y,t)). """ tac = sorted([(t, v) for t, v in tissue_tac if t is not None and v is not None]) inp = sorted([(t, v) for t, v in input_func if t is not None and v is not None]) if len(tac) < 3 or len(inp) < 2: return None it = [p[0] for p in inp] iv = [p[1] for p in inp] # Cp 누적적분 (사다리꼴) — 조직 프레임 시간에서 평가 times = [t for t, _ in tac] # 적분은 0..t. 미세 그리드로 누적. cum = {} integral = 0.0 grid = sorted(set([0.0] + times)) prev_t = 0.0 prev_cp = _interp(it, iv, 0.0) running = 0.0 int_at = {0.0: 0.0} for g in grid[1:]: cp_g = _interp(it, iv, g) running += 0.5 * (prev_cp + cp_g) * (g - prev_t) int_at[g] = running prev_t, prev_cp = g, cp_g xs_all, ys_all, ts_all = [], [], [] for t, ct in tac: cp = _interp(it, iv, t) if cp <= 0: continue xs_all.append(int_at[t] / cp) ys_all.append(ct / cp) ts_all.append(t) if len(xs_all) < 3: return None # 선형 구간 선택 chosen_tstar = tstar_min if tstar_min is not None else DEFAULT_TSTAR_MIN if auto_tstar and tstar_min is None: # 후보 t* 들 중 r2 최대(최소 3점)인 구간 선택 best = None candidate_ts = sorted(set(ts_all)) for cand in candidate_ts: idx = [i for i, tt in enumerate(ts_all) if tt >= cand] if len(idx) < 3: continue sx = [xs_all[i] for i in idx] sy = [ys_all[i] for i in idx] slope, intercept, r2 = _linreg(sx, sy) if slope is None or r2 is None: continue if best is None or r2 > best[0]: best = (r2, cand, slope, intercept, len(idx)) if best: r2, chosen_tstar, slope, intercept, npts = best return {"Ki": slope, "intercept": intercept, "r2": r2, "tstar": chosen_tstar, "n_points": npts, "points": list(zip(xs_all, ys_all, ts_all))} idx = [i for i, tt in enumerate(ts_all) if tt >= chosen_tstar] if len(idx) < 2: idx = list(range(len(ts_all))) sx = [xs_all[i] for i in idx] sy = [ys_all[i] for i in idx] slope, intercept, r2 = _linreg(sx, sy) if slope is None: return None return {"Ki": slope, "intercept": intercept, "r2": r2, "tstar": chosen_tstar, "n_points": len(idx), "points": list(zip(xs_all, ys_all, ts_all))} def mrglu(ki, glucose_mgdl, lumped_constant): """조직 포도당 대사율 MRGlu = (Ki x 혈당) / LC. 단위: Ki mL/min/mL, glucose mg/dL -> mg/(min·100mL)? 여기서는 관례적으로 MRGlu = Ki[mL/min/g] * Glu[mg/dL] / LC, glucose 를 mg/dL/100 = mg/mL 환산해 umol 계열 대신 상대 정량(연구간 비교)용 지표로 출력한다. 출력 단위: (mg glucose) / (min · 100 mL tissue) 근사.""" if ki is None or glucose_mgdl is None or not lumped_constant: return None return (ki * glucose_mgdl) / lumped_constant def kinetic_table(tac_rows, input_funcs, meta, lumped_constants=None): """animal x timepoint x tissue 별 Patlak Ki + MRGlu 테이블.""" if lumped_constants is None: lumped_constants = DEFAULT_LUMPED_CONSTANT out = [] keys = sorted({(r["animal_id"], r["timepoint"], r["tissue"]) for r in tac_rows}) for aid, tp, tissue in keys: tac = [(r["time_min"], r.get("decay_corrected_bqml")) for r in tac_rows if r["animal_id"] == aid and r["timepoint"] == tp and r["tissue"] == tissue] inp = input_funcs.get((aid, tp)) m = meta.get((aid, tp)) if not inp or not m: continue inp_pairs = [(d["time_min"], d["plasma_activity_bqml"]) for d in inp] res = patlak(tac, inp_pairs) if res is None: continue lc = lumped_constants.get(tissue, DEFAULT_LC_FALLBACK) mr = mrglu(res["Ki"], m["glucose_mgdl"], lc) out.append({ "animal_id": aid, "group": m["group"], "timepoint": tp, "tissue": tissue, "Ki": res["Ki"], "patlak_r2": res["r2"], "tstar": res["tstar"], "n_points": res["n_points"], "lumped_constant": lc, "glucose_mgdl": m["glucose_mgdl"], "MRGlu": mr, }) return out # --------------------------------------------------------------------------- # 종단 within-animal (paired baseline vs post) # --------------------------------------------------------------------------- def longitudinal_changes(metric_rows, value_key): """동일 개체 baseline vs post paired 변화. metric_rows: dict 들 (animal_id, group, timepoint, tissue, value_key 포함). 반환: list[dict] (animal_id, group, tissue, baseline, post, delta, pct_change).""" by = {} for r in metric_rows: by.setdefault((r["animal_id"], r["group"], r["tissue"]), {})[r["timepoint"]] = r.get(value_key) out = [] for (aid, group, tissue), tps in sorted(by.items()): b = tps.get("baseline") p = tps.get("post") if b is None or p is None: continue delta = p - b pct = (delta / b * 100.0) if b else None out.append({"animal_id": aid, "group": group, "tissue": tissue, "baseline": b, "post": p, "delta": delta, "pct_change": pct}) return out # --------------------------------------------------------------------------- # cohort 통계 (stdlib paired t-test) # --------------------------------------------------------------------------- def _mean(xs): xs = [x for x in xs if x is not None] return sum(xs) / len(xs) if xs else None def _std(xs): xs = [x for x in xs if x is not None] n = len(xs) if n < 2: return None mu = sum(xs) / n return math.sqrt(sum((x - mu) ** 2 for x in xs) / (n - 1)) def _student_t_sf(t, df): """two-sided p-value 근사 (Student t, stdlib only). incomplete beta 기반 정확식.""" if df <= 0: return None t = abs(t) x = df / (df + t * t) # regularized incomplete beta I_x(df/2, 1/2) p = _betai(df / 2.0, 0.5, x) return p # two-sided def _betacf(a, b, x): MAXIT = 200; EPS = 3e-12; FPMIN = 1e-300 qab = a + b; qap = a + 1.0; qam = a - 1.0 c = 1.0; d = 1.0 - qab * x / qap if abs(d) < FPMIN: d = FPMIN d = 1.0 / d; h = d for m in range(1, MAXIT + 1): m2 = 2 * m aa = m * (b - m) * x / ((qam + m2) * (a + m2)) d = 1.0 + aa * d if abs(d) < FPMIN: d = FPMIN c = 1.0 + aa / c if abs(c) < FPMIN: c = FPMIN d = 1.0 / d; h *= d * c aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)) d = 1.0 + aa * d if abs(d) < FPMIN: d = FPMIN c = 1.0 + aa / c if abs(c) < FPMIN: c = FPMIN d = 1.0 / d; de = d * c; h *= de if abs(de - 1.0) < EPS: break return h def _betai(a, b, x): if x <= 0.0: return 0.0 if x >= 1.0: return 1.0 lbeta = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) bt = math.exp(lbeta + a * math.log(x) + b * math.log(1.0 - x)) if x < (a + 1.0) / (a + b + 2.0): return bt * _betacf(a, b, x) / a return 1.0 - bt * _betacf(b, a, 1.0 - x) / b def paired_ttest(pairs): """pairs: list[(baseline, post)] -> dict(n, mean_delta, t, df, p).""" deltas = [p - b for b, p in pairs if b is not None and p is not None] n = len(deltas) if n < 2: return {"n": n, "mean_delta": _mean(deltas), "t": None, "df": None, "p": None} md = sum(deltas) / n sd = math.sqrt(sum((d - md) ** 2 for d in deltas) / (n - 1)) if sd == 0: return {"n": n, "mean_delta": md, "t": None, "df": n - 1, "p": None} se = sd / math.sqrt(n) t = md / se df = n - 1 p = _student_t_sf(t, df) return {"n": n, "mean_delta": md, "t": t, "df": df, "p": p} def cohort_summary(metric_rows, value_key): """group x tissue x timepoint 별 mean/std/n + (baseline vs post) paired t-test.""" summary = {} for r in metric_rows: k = (r["group"], r["tissue"], r["timepoint"]) summary.setdefault(k, []).append(r.get(value_key)) rows = [] for (group, tissue, tp), vals in sorted(summary.items()): rows.append({"group": group, "tissue": tissue, "timepoint": tp, "n": len([v for v in vals if v is not None]), "mean": _mean(vals), "std": _std(vals)}) # paired tests per group x tissue tests = [] by = {} for r in metric_rows: by.setdefault((r["animal_id"], r["group"], r["tissue"]), {})[r["timepoint"]] = r.get(value_key) grp = {} for (aid, group, tissue), tps in by.items(): if "baseline" in tps and "post" in tps: grp.setdefault((group, tissue), []).append((tps["baseline"], tps["post"])) for (group, tissue), pairs in sorted(grp.items()): res = paired_ttest(pairs) res.update({"group": group, "tissue": tissue}) tests.append(res) return {"descriptive": rows, "paired_tests": tests} # --------------------------------------------------------------------------- # mechanism 분류 # --------------------------------------------------------------------------- def classify_mechanism(long_changes_ki): """treatment 군의 종단 Ki(또는 SUV) 변화로 mechanism 분류. 입력: longitudinal_changes() 결과 (Ki 기준). 규칙(treatment 군 평균 pct_change): - skeletal_muscle 상승 우세 -> '근육 흡수 개선형' - BAT 상승 우세 -> 'BAT 활성형' - myocardium 변화 우세 -> '심근 변화형' """ tx = [r for r in long_changes_ki if r["group"] == "treatment"] agg = {} for r in tx: agg.setdefault(r["tissue"], []).append(r["pct_change"]) means = {t: _mean(v) for t, v in agg.items()} candidates = { "근육 흡수 개선형 (skeletal muscle uptake)": means.get("skeletal_muscle"), "BAT 활성형 (brown adipose activation)": means.get("BAT"), "심근 변화형 (myocardial shift)": means.get("myocardium"), } valid = {k: v for k, v in candidates.items() if v is not None} if not valid: return {"label": "분류 불가 (treatment paired 데이터 부족)", "scores": candidates} label = max(valid, key=lambda k: valid[k]) return {"label": label, "scores": candidates, "tissue_pct_change": means} # --------------------------------------------------------------------------- # 참고 범위 / 누락 경고 # --------------------------------------------------------------------------- def reference_check(suv_rows, model="C57BL/6"): """SUV_bw 가 모델 참고 범위를 벗어나는 항목 경고.""" warnings = [] ref = MODEL_REFERENCE.get(model) if not ref: return ["참고 모델 '%s' 미정의" % model] for r in suv_rows: tr = ref.get(r["tissue"]) if not tr or r["SUV_bw"] is None: continue lo, hi = tr["suv"] if r["SUV_bw"] < lo or r["SUV_bw"] > hi: warnings.append( "참고범위 이탈 [%s] %s/%s %s: SUV_bw=%.2f (%s 참고 %.2f-%.2f)" % (model, r["animal_id"], r["timepoint"], r["tissue"], r["SUV_bw"], model, lo, hi)) return warnings def metadata_warnings(meta, input_funcs, lumped_constants, tac_rows): """input function / lumped constant / BAT 조건 누락 경고.""" warns = [] tissues = sorted({r["tissue"] for r in tac_rows}) for tissue in tissues: if tissue not in (lumped_constants or DEFAULT_LUMPED_CONSTANT): warns.append("lumped constant 누락: '%s' -> fallback %.2f 사용" % (tissue, DEFAULT_LC_FALLBACK)) for key, m in meta.items(): if key not in input_funcs: warns.append("input function 누락: %s/%s -> Patlak Ki 산출 불가" % (key[0], key[1])) if m.get("scan_temp_c") in (None, ""): warns.append("BAT 촬영 온도(scan_temp_c) 미기록: %s/%s" % (key[0], key[1])) if not m.get("anesthesia"): warns.append("마취(anesthesia) 조건 미기록: %s/%s" % (key[0], key[1])) return warns