Python port měření intenzity slave kanálu z balíku cmeAnalysis (DanuserLab, MATLAB), vyvinutý pro měření dynaminu na drahách klatrinových jamek.
Pro zadané souřadnice [snímek, y, x] a TIRF video vrátí intenzitu dynaminu
spočtenou stejným postupem jako cmeAnalysis: fit skvrnky tvaru PSF s pevnou
šířkou, výsledkem je amplituda nad lokálním pozadím, její nejistota
a p-hodnota. Čistě numpy/scipy — žádná kompilovaná rozšíření, žádný MATLAB;
běží na Windows, Linuxu i obou generacích Maců. Numerická shoda s originální
zkompilovanou binárkou je změřená (viz Validace).
Obsah: Instalace · Rychlý start · Formát vstupů · Pipeline krok za krokem · Přehled skriptů · Jak port vznikal · Dokumentace měření · Poznámky k našim datasetům · Co v repozitáři není · Testy · Licence
git clone https://github.com/gulierus/CMEpython.git
cd CMEpython
pip install -e .Jádro (cmepython/) potřebuje jen numpy, scipy, tifffile.
Analytické skripty (scripts/compare_*, scripts/report_*, kohorty) navíc
pandas a matplotlib:
pip install pandas matplotlibPDF verze protokolů vyžaduje LaTeX (latexmk + xelatex); bez něj skripty
doběhnou s --skip-pdf nebo vypíší varování a nechají markdown a .tex.
Testy: pip install -e ".[dev]" a pytest.
Jeden bod:
from cmepython import dynamin_intensity
r = dynamin_intensity(video, frame, y, x,
sigma_slave=1.4356, # šířka PSF měřicího kanálu [px]
sigma_master=1.6424) # šířka PSF kanálu s detekcemi [px]
r["A"] # amplituda gaussovky nad lokálním pozadím — hledaná intenzita
r["c"] # fitované lokální pozadí
r["A_pstd"] # nejistota amplitudy
r["pval_Ar"] # p-hodnota testu proti šumu
r["hval_Ar"] # bool: signál je významnýVíc bodů najednou — čte snímky lazy z disku a volitelně paralelizuje:
from cmepython import measure_movie
import numpy as np
coords = np.array([[frame, y, x], ...]) # (N, 3)
res = measure_movie("film.tif", coords,
sigma_slave=1.4356, sigma_master=1.6424,
slave_channel=1, master_channel=0,
workers=0) # 0 = všechna jádra, None/1 = sériově
res["A"] # pole intenzit, jeden prvek na řádek coordsPod
spawn(macOS, Windows) musí být paralelní volání ve skriptu sif __name__ == "__main__":. Z interaktivní session spadne na sériový běh s varováním, ne s chybou.
Filmy: vícestránkový TIFF. Osy se čtou z metadat; podporované rozložení
TCYX, ZCYX, CYX, ZYX i YX (cmepython.measure.movie_layout). Kanály se
adresují indexem (slave_channel = měřený kanál, master_channel = kanál,
ve kterém vznikly detekce).
Trajektorie: CSV s minimálně sloupci particle (id dráhy), frame
(index snímku, 0-based), x, y (pixely; x = sloupec, y = řádek,
stejná konvence jako v obrazových polích). Volitelné sloupce (cls = shape
index, cluster, …) projdou do výstupu beze změny.
Párování film ↔ CSV: skripty párují soubory podle čísla filmu v názvu
(regulární výraz _(\d+)[-_]), např. ..._16_LRR.tif ↔
..._16-lr-trajectories.csv.
Celý postup od syrových filmů k výsledkům. Každý skript má --help
se všemi volbami.
Šířka σ se odhaduje z dat, ne z optiky — stejně jako výchozí chování
cmeAnalysis (getGaussianPSFsigmaFromData.m): detekce s pevnou σ, volný fit
xyasc, filtrace neúspěšných fitů, směs gaussovek 1–3 komponent podle BIC,
vybere se komponenta s nejvyšším vrcholem. Vzorkuje ~40 snímků napříč
všemi filmy podmínky (jedna σ na kanál na podmínku, jako cmeAnalysis):
python3 scripts/calibrate_dataset.py "lr registered" --out psf_calibration_lr.jsonVýstupní JSON obsahuje σ pro každý kanál včetně diagnostiky (BIC, komponenty,
počty spotů). cmeAnalysis hodnoty pod 1,1 px zvedá na 1,1 (runDetection.m:90);
port to dělá také a hlásí to.
python3 scripts/measure_trajectories.py \
--traj-dir "lr registered/trajectories" --movie-dir "lr registered" \
--out-dir "lr registered/measured" \
--sigma-slave 1.4356 --sigma-master 1.6424 \
--slave-channel 1 --master-channel 0 --also-master --workers 0Ke sloupcům vstupního CSV přibudou (prefix dnm_ = slave/dynamin,
clc_ = master/clathrin při --also-master):
| sloupec | význam |
|---|---|
dnm_A, dnm_c |
amplituda nad lokálním pozadím a pozadí [jednotky kamery] |
dnm_A_pstd |
nejistota amplitudy (šířená z fitu, stupně volnosti jako MEX) |
dnm_sigma_r |
směrodatná odchylka reziduí fitu |
dnm_pval, dnm_signif |
p-hodnota a verdikt testu signálu proti šumu (α = 0,05) |
dnm_x, dnm_y |
doladěná poloha (volný fit přijat při posunu < 3σ_master a vyšší A) |
dnm_valid |
fit uvnitř obrazu (False jen u okna přes okraj) |
track_len |
délka dráhy ve snímcích |
Prázdné místo není chyba: vrací A ≈ 0 a p ≈ 1. Měření a rozhodnutí o pozitivitě jsou oddělené kroky.
Před měřením proběhnou vstupní kontroly (validate=True): detekce
SIM pruhů ve spektru (na jednotlivém raw TIRF-SIM snímku fit tiše lže;
řešení je zprůměrovat devítice) a kontrola registrace kanálů FFT křížovou
korelací (posun nad 1 px varuje). Nálezy jsou varování, ne chyby.
Port Aguetovy statistické klasifikace (runSlaveChannelClassification.m).
Dva testy na dráhu, oba binomicky korigované na její délku: significant_master
(počet významných detekcí proti náhodě dané p_detection filmu)
a significant_slave (počet snímků s amplitudou nad 95. percentilem pozadí
filmu, t-test; výchozí režim „s" navazujících nástrojů cmeAnalysis).
python3 scripts/classify_trajectories.py \
--measured "lr registered/measured" --movies "lr registered" \
--masks "lr registered/masks" \
--sigma-slave 1.4356 --slave-channel 1 --master-channel 0
# citlivost testu: --alpha 0.01 jiny vystup: --out-suffix "-classified-a0.01.csv"Výstup na dráhu: particle, track_len, n_detected, thr_master, significant_master, n_above_bg, thr_slave, significant_slave, max_A, max_master_A, max_si. Pozadí filmu se počítá z míst uvnitř buněčné masky
dál než 4σ od všech detekcí; masky lze dodat (--masks), jinak se odhadnou
z maximální projekce.
Výchozí cesty skriptu míří na první dataset; pro jiný vždy zadejte
--measured/--movies/--masks/--sigma-slave/--slave-channel/--master-channel.
Port getIntensityCohorts.m + statistik z plotIntensityCohorts.m: dráhy
se rozdělí podle životnosti (meze 10/20/40/60/80/100/120 s), převzorbují
kubicky na střední délku kohorty včetně 5 bufferových snímků před vznikem
a po zániku (měřeno gap cestou portu interpTrack), a průměrují — nejdřív
v každém filmu, pak přes filmy (SEM přes filmy, jako cmeAnalysis):
python3 scripts/cohort_analysis.py \
--measured "lr registered/measured" --movies "lr registered" \
--framerate 2 --sigma-slave 1.4356 --sigma-master 1.6424 \
--slave-channel 1 --master-channel 0 --si-col cls --out out/lr_cohortsVýstup: cohort_curves.csv a dva grafy (kohorty; produktivní vs. abortivní
podle max SI > 0,7).
Tři skripty na sebe navazují:
python3 scripts/compare_classifications.py # confusion matrix + JSON
python3 scripts/compute_boxmean.py # box-mean 5×5 readout na týchž drahách
python3 scripts/report_classification_comparison.py # kompletní protokol (md + tex + pdf)Protokol obsahuje mozaiku confusion matic po délkových pásmech, kompletní
metrikovou tabulku s 95% intervaly bootstrapem po filmech, length-adjusted
OR (Mantel–Haenszel), within-band AUC, sweep přes α, Hodgesovy–Lehmannovy
posuny a van Elterenův test. Sweep přes α vyžaduje předpočítané klasifikace
(--alpha A --out-suffix=-classified-aA.csv, viz krok 3).
scripts/build_cme_corpus.py převede naměřené dráhy do formátu korpusu
Shape2Fate (dvě varianty: surová amplituda a amplituda dělená biexponenciálním
fitem průměrů snímků buňky; konvence intenzita = 1 + amplituda, práh
thresholds_q90.json). Trénink samotný běží kódem spolupracujícího projektu
(větev release/dynamin-confusion-v1 repozitáře Shape2Fate_Fake2Emulate)
na výpočetním clusteru; srovnávací protokol generuje
scripts/report_detector_comparison.py.
matlab -batch export_fits # matlab/export_fits.m (jen jednou)
python3 scripts/validate_against_matlab.pyVyžaduje MATLAB s toolboxy Image Processing a Statistics; referenční výstup
MEX binárky pro 874 oken je ale commitnutý (matlab/mex_reference.mat),
takže samotné srovnání běží i bez MATLABu.
Jádro pipeline (kroky 1–7):
| skript | účel |
|---|---|
calibrate_dataset.py |
odhad σ PSF z dat, celý dataset najednou |
measure_trajectories.py |
měření intenzity na drahách → *-dynamin.csv |
classify_trajectories.py |
Aguetova klasifikace drah → *-classified.csv |
cohort_analysis.py |
lifetime kohorty, křivky + grafy |
plot_cohorts_kamenik_style.py |
kohorty ve stylu plotIntensityCohorts (dynamin+ / dynamin−) |
compare_classifications.py |
confusion matrix SI vs. cmeAnalysis |
compute_boxmean.py |
box-mean 5×5 readout na pozicích drah |
report_classification_comparison.py |
protokol úlohy B (md/tex/pdf) |
build_cme_corpus.py |
korpusy pro trénink detektoru (2 varianty) |
build_clathrin_corpus.py |
klatrinový korpus jako pozitivní kontrola metody |
report_detector_comparison.py |
protokol úlohy A: srovnání detektorů |
validate_against_matlab.py |
numerické srovnání s MEX binárkou |
dataset_findings.py |
reprodukce datasetových nálezů (bleaching, FP rate) |
bench_measure.py |
benchmark měření |
Klatrinový experiment (jednokanálový dataset *_AVG.tif):
| skript | účel |
|---|---|
sample_clathrin_avg.py |
krok 1: klatrinová intenzita (PSF fit + box-mean) na pozicích drah |
build_clavg_corpus.py |
krok 2: korpusy klatrin → dynaminový štítek |
clavg_threshold_analysis.py |
prahová analýza klatrinové intenzity proti dynaminovému štítku |
report_clavg_experiment.py |
protokol klatrinového experimentu |
fig_logreg_matrices.py |
mozaika confusion matic logistické regrese (--lang cz|en) |
fig_summary_en.py |
souhrnný sloupcový graf (anglicky) |
Experiment dynamin → klatrin (štítek z klatrinu místo SI):
| skript | účel |
|---|---|
sample_clc_his_corpus.py |
klatrinová intenzita na pozicích korpusu spolupracujícího projektu |
dyn_to_clc_experiment.py |
logistická regrese s klatrinovým štítkem (práh podle prevalence a GMM) |
fig_dyn_to_clc.py |
obrázky experimentu (anglicky) |
report_logreg_interpretace.py |
interpretace referenční analýzy logistickou regresí |
Skripty obou experimentů čtou kód a korpus spolupracujícího projektu
z external/Shape2Fate_Fake2Emulate/ (větev release/dynamin-confusion-v1,
v repozitáři není) a potřebují navíc scikit-learn.
Nese průběh dynaminu informaci o osudu jamky?
| skript | účel |
|---|---|
shape_phase1.py |
tvar křivky vs. délka vs. úroveň (tři modely) |
shape_phase2.py |
model z množství a významnosti dynaminu |
shape_phase3_cnn.py |
kontrola 1D konvoluční sítí nad surovými stopami (vyžaduje torch) |
report_shape_story.py |
protokol všech tří testů |
Patra 1–2 potřebují scikit-learn, patro 3 torch.
Klíčová funkce cmeAnalysis fitGaussian2D existuje jen jako zkompilovaná
MEX binárka bez zdrojáku (Levenberg–Marquardt z GSL). Port ji proto
rekonstruuje podle chování a shodu měří přímo proti binárce na 874
oknech z reálných dat: amplitudy, pozadí a rezidua sedí na medián relativního
rozdílu 10⁻⁸ až 10⁻¹⁶, lokalizace na ~5×10⁻⁷ px (max ~10⁻³; poloha je nejhůř
podmíněný parametr). Akceptační test volného fitu propouští stejná okna
(567 z 874). Vyšší vrstvy (klasifikace, kohorty, kalibrace) jsou přepsané
řádek po řádku podle otevřených .m zdrojů; komentáře v kódu odkazují
soubor.m:řádek.
Poznatky, které při portování stály nejvíc času (a v kódu jsou zohledněné):
- **
padarrayXT('symmetric')je numpy'reflect'**, ne numpy'symmetric'` — MATLAB zrcadlí bez duplikace okrajového pixelu. Záměna dá chyby ~10³ v pásu u okraje. - MEX počítá stupně volnosti
npx − n_free − 1, o jeden méně než běžná konvence. Bez převzetí jsouA_pstda všechny p-hodnoty o ~0,1 % vyšší. sum(pval<Alpha)/sum(cellmask)na 2D maticích je v MATLABu maticové dělení (mrdivide), ne podíl počtů. Port reprodukuje skutečně vykonávaný výpočetp_detection.- σ se na slave kanálu nikdy nefituje — módy
Ac/xyAcmají šířku pevnou; volný fitxyascslouží jen kalibraci. σ < 1,1 px se zvedá na 1,1. - Poloha není striktně fixní: po fitu s pevnou polohou se zkusí volný
a přijme se při posunu < 3σ_master a vyšší amplitudě (
runDetection.m:180-190). - Jedna vědomá odchylka: CCS oblasti se z pozadí vylučují disky kolem pozic
z trajektorií místo detekčních masek
dmasks.tif, které bez běhu celé MATLAB pipeline neexistují.
Mapa modulů na originál:
| modul | originál v cmeAnalysis |
|---|---|
slave_intensity.dynamin_intensity |
runDetection.m:180-190 + fitGaussians2D.m |
slave_intensity.dynamin_intensity_gap |
runTrackProcessing.m:879-926 (interpTrack) |
psf_calibration.estimate_psf_sigma |
getGaussianPSFsigmaFromData.m |
classification.background_stats / classify_track |
runSlaveChannelClassification.m |
cohorts.cohort_curves |
getIntensityCohorts.m + plotIntensityCohorts.m |
background.filter_gaussian_fit_2d, mask_from_first_mode |
filterGaussianFit2D.m, maskovací větev |
edf_scaling.scale_edfs |
scaleEDFs.m |
validation.* |
vlastní vstupní kontroly (SIM pattern, registrace) |
Upstream si naklonuj zvlášť, v repozitáři není:
git clone https://github.com/DanuserLab/cmeAnalysis.gitSrozumitelný popis, co program s daty dělá — psaný i pro čtenáře bez technického zázemí (vyhledávání blobu, fitování a co je výsledná intenzita, selhání fitu, normalizace, globální statistiky):
- docs/how-dynamin-is-measured.md — anglicky, markdown
- docs/jak-merime-dynamin.html — česky, stránka s diagramy (otevři v prohlížeči)
Hodnoty σ pro oba datasety projektu jsou commitnuté
(psf_calibration.json, psf_calibration_lr.json), aby čísla ve skriptech
nebyla magické konstanty.
Dataset 1, reconstructed registered/ (TIRF-SIM, 1024², 39,55 nm/px):
ch0 = clathrin (SIM rekonstrukce), ch1 = dynamin (SIM rekonstrukce),
ch2 = zprůměrovaný raw dynamin TIRF, 2× upsamplovaný, registrovaný — na něm
se měří (σ 2,6424 px). SIM rekonstrukce mají FWHM pod difrakčním limitem
a PSF model na ně nepatří. Interval mezi snímky nedohledán (odhad 1,5–3 s).
Dataset 2, lr registered/ (512², 79,1 nm/px, 2 s/snímek, NA 1,5):
ch0 = clathrin (σ 1,6424 px), ch1 = dynamin (σ 1,4356 px); oba kanály
registrované, dynamin je průměr 9 SIM expozic (simulovaný TIRF).
Ověřené vlastnosti měření na těchto datech (reprodukce
scripts/dataset_findings.py):
- Rozptyl mezi filmy je převážně biologický — dynaminové amplitudy
kolísají 17×, pozadí ~2,3×, clathrin ~1,5–2×. EDF škálování
(
scale_edfs) by smazalo měřený efekt; používat vědomě. - Photobleaching: ~12–17 % na 100 snímků pro dráhy ≥ 15 snímků; původní vyšší odhady nafukovala levá cenzura. Port korekci nezavádí (cmeAnalysis ji také nemá); vhodnější je čas vzniku dráhy jako kovariáta.
- Test významnosti je konzervativní: na prázdném pozadí ~0,1–0,3 % falešně pozitivních při nominálních 5 %. U 2× upsamplovaného kanálu není šum sousedních pixelů nezávislý; korekce na efektivní počet pixelů mění < 1 % rozhodnutí.
Repozitář obsahuje jen kód a dokumentaci. Lokálně (mimo git) žijí:
| co | kde lokálně |
|---|---|
| filmy a masky | reconstructed registered/, lr registered/ |
| trajektorie a všechna naměřená CSV | lr registered/trajectories/, lr registered/measured/ |
| protokoly s výsledky (úlohy A i B) | lr registered/comparison_report/, CME_for_Helios/report/ |
| kohortové výstupy | lr registered/cohorts/, out/ |
klatrinový dataset (*_AVG.tif), nasamplované intenzity a protokoly experimentů |
Clathrin Analysis/ |
| balík pro výpočetní cluster (korpusy + cizí kód) | CME_for_Helios/ |
| klon spolupracujícího projektu Shape2Fate | external/ |
| klon cmeAnalysis a publikace | cmeAnalysis/, Aguet13.pdf, mmc1.pdf |
| velké validační reference | matlab/filter_ref.mat aj. (přegenerují se) |
python3 -m pytest tests/ -qTesty ověřují chování odvozené ze zdrojáku cmeAnalysis a vnitřní konzistenci
(dávka == bod po bodu). Numerickou shodu s MATLABem ověřuje zvlášť
scripts/validate_against_matlab.py (viz krok 7 pipeline).
GPL-3.0-or-later. cmeAnalysis je GPL-3.0, odvozený port proto také.