From 5fe6298926e1afe7e317e10420b2a6a69cc511ae Mon Sep 17 00:00:00 2001 From: Matiq Date: Fri, 14 Aug 2026 20:27:31 +0300 Subject: [PATCH] README: methodology, usage, outputs, classification, validation --- README.md | 138 +++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 137 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 3d7e184..d1f4a13 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,139 @@ # ffpe_damage -FFPE deamination QC: normalized C>T/G>A profiling from BAM + control-based classification \ No newline at end of file +Методологически корректная оценка дезаминирования в FFPE-образцах (формалин-фиксированная парафинизированная ткань) по BAM. + +Инструмент **не** изобретает универсальный «FFPE damage score». Он выдаёт *измеряемые* показатели: + +- все 12 типов mismatch; +- нормализованную частоту `C>T` и `G>A` на каждой позиции read; +- отдельно **R1/R2**; +- отдельно **forward/reverse** alignment; +- расстояния от **5′- и 3′-конца** read; +- профиль частот (не просто количество mismatch); +- BED-треки для IGV; +- CSV/TSV для сравнения образцов между собой. + +Классификация `low / moderate / high` назначается **только** при сравнении с валидированными контрольными образцами, приготовленными тем же протоколом (`ffpe_compare.py`). Абсолютных порогов типа «>5% = плохо» нет — у них нет валидированного основания. + +## Что считаем + +Для каждой позиции *i* в read: + +``` +C>T_i = #(C→T на позиции i) / #(C на позиции i) +G>A_i = #(G→A на позиции i) / #(G на позиции i) +``` + +Знаменатель — **все** высококачественные наблюдения соответствующего reference-основания на этой позиции, **включая совпадения** (matches). Это принципиально: только так частота сопоставима между позициями и образцами. + +Фильтрация: убираются unmapped/secondary/supplementary, `MAPQ < --mapq`, основание с `BQ < --baseq`, дупликаты (кроме `--include-duplicates`). Read-позиция отсчитывается как от 5′-конца (`position_5prime`), так и от 3′-конца (`position_3prime`). + +## Установка + +```bash +cd ~/Projects/ffpe_damage +python3 -m venv --system-site-packages .venv +.venv/bin/pip install pysam==0.24.0 # плюс pandas, matplotlib, scipy +``` + +Требования на входные данные: +- BAM отсортирован и проиндексирован (`.bai`); +- референсная FASTA проиндексирована (`.fai`). + +## 1. Анализ образца + +```bash +.venv/bin/python ffpe_damage_v2.py \ + --bam sample.bam \ + --reference hg38.fa \ + --outdir FFPE_sample +``` + +Строгий вариант: + +```bash +.venv/bin/python ffpe_damage_v2.py \ + --bam sample.bam --reference hg38.fa --outdir FFPE_sample \ + --mapq 30 --baseq 20 --min-depth 30 --min-alt-count 5 +``` + +### Выходные файлы + +| Файл | Содержимое | +|---|---| +| `sample_summary.csv` | число прочитанных/использованных ридов, usable bases, глобальная частота C>T/G>A, пороги фильтрации | +| `substitution_summary.csv` | все 12 типов замен: count, знаменатель, частота | +| `substitution_all12.png` | бар-чарт частот всех 12 типов (C>T и G>A выделены) | +| `normalized_damage_by_read_position.csv` | полный профиль: R1/R2 × позиция × reference-base × замена | +| `normalized_damage_by_read_position_3prime.csv` | то же, по расстоянию от 3′-конца | +| `CtoT_GtoA_normalized_profile.csv` | ключевая таблица «позиция в read vs C>T% / G>A%» | +| `CtoT_GtoA_normalized_profile_3prime.csv` | то же по 3′-концу | +| `strand_damage_profile.csv` | C>T/G>A по цепи выравнивания (`+`/`-`) | +| `strand_damage_profile_3prime.csv` | то же по 3′-концу | +| `end_enrichment.csv` | C>T/G>A % в первых/последних 1,3,5,10,20 основаниях, по R1/R2 | +| `candidate_CtoT.bed`, `candidate_GtoA.bed` | позиции для IGV; имя = `ref>alt;depth;alt;VAF` | +| `candidate_variants.tsv` | контекст каждой позиции: VAF, средняя позиция в read, strand bias, BQ, раскладка R1/R2 | +| `CtoT_R1.png` … `GtoA_R2_3prime.png` | графики частоты по позиции в read (5′ и 3′) | + +### Как читать результат + +- **Классическая FFPE-сигнатура**: пик `C>T` у 5′-конца R1 и `G>A` у 3′-конца R2 — две цепи одного и того же дезаминирования цитозина. Проверяйте профили R1/R2 и по strand. +- Если C>T/G>A повышены равномерно по всей длине read (без пика у концов) — это не классическое end-associated повреждение, а другой артефакт/ошибки. +- **Отличить настоящий вариант от дезаминирования** (`candidate_variants.tsv`): настоящий вариант имеет ожидаемый VAF, распределён по read (средняя позиция ~ середина), встречается на обеих цепях (strand bias не значим) и не привязан к концам. Дезаминирование концентрируется у концов read и на одной цепи. + +## 2. Классификация относительно контролей + +```bash +.venv/bin/python ffpe_compare.py \ + --controls sampleA_out sampleB_out sampleC_out \ + --query query_sample_out \ + --out compare_out \ + --window 5 --k 2 +``` + +`--controls` принимает несколько папок с выходами шага 1 (в них лежит `end_enrichment.csv`) либо один файл-манифест с путями по одному в строке. + +Метрики (окно `W` от конца read): + +- `C_to_T_R1_5prime_W` — C>T в первых W основаниях R1 (forward); +- `G_to_A_R2_3prime_W` — G>A в последних W основаниях R2 (reverse); +- `combined` — среднее двух. + +Классификация против среднего M и SD S контролей, `z = (query − M)/S`: + +``` +high : z > k +moderate : -k ≤ z ≤ k +low : z < -k +``` + +Выход: `classification.csv`, `comparison_table.csv`, `control_vs_query.png`. + +> Зачем сравнивать только с контролями? FFPE-повреждение зависит от протокола (UDG-обработка, end repair, длина фрагментов, панель). Одно и то же значение частоты может быть нормой для одного протокола и тяжёлым повреждением для другого. Поэтому сравнивать стоит только образцы, приготовленные одинаково, а классификация всегда относительная. + +## 3. Синтетические данные для проверки + +```bash +.venv/bin/python make_test_data.py --outdir test_data/damaged --damage-rate 0.5 +.venv/bin/python make_test_data.py --outdir test_data/clean --damage-rate 0.0 + +.venv/bin/python ffpe_damage_v2.py \ + --bam test_data/damaged/sample.sorted.bam \ + --reference test_data/damaged/synthetic_reference.fa \ + --outdir test_data/damaged/out +``` + +Встроенный ground truth: +- damaged: C>T у 5′-конца forward-ридов (R1), G>A у 3′-конца reverse-ридов (R2), плюс один «настоящий» гетерозиготный SNP без read-position/strand bias — его инструмент **не** должен классифицировать как дезаминирование; +- clean: только фоновая ошибка секвенирования. + +Ожидаемые результаты на damaged: +- R1, 5′-профиль: C>T ~48% на позициях 1–5, ~0 дальше; +- R2, 3′-профиль: G>A ~50% в последних 5 основаниях, ~0 дальше; +- SNP на `chr1:2504`: VAF ~50%, середина read, обе цепи, без strand bias. + +## Оговорки + +- Это QC/исследовательский инструмент, а не валидированный клинический assay. +- C>T/G>A в обычном alignment нельзя автоматически интерпретировать как химическое дезаминирование при поиске соматических вариантов: нужны read-position/strand/base-quality, а в идеале — matched normal и сопоставление с variant calls. +- Сравнивать нужно только образцы одного протокола подготовки библиотеки.