From fb99880b80fc416392bcb760d85d7e28c09614d1 Mon Sep 17 00:00:00 2001 From: Matiq Date: Mon, 7 Sep 2026 01:05:49 +0300 Subject: [PATCH] README: threads, VCF filtering, annotation (GENCODE/OncoKB), VAF filter --- README.md | 77 ++++++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 68 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 2af70a3..7095897 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ - расстояния от **5′- и 3′-конца** read; - профиль частот (не просто количество mismatch); - BED-треки для IGV; -- CSV/TSV для сравнения образцов между собой. +- CSV для сравнения образцов между собой. Классификация `low / moderate / high` назначается **только** при сравнении с валидированными контрольными образцами, приготовленными тем же протоколом (`ffpe_compare.py`). Абсолютных порогов типа «>5% = плохо» нет — у них нет валидированного основания. @@ -33,12 +33,29 @@ G>A_i = #(G→A на позиции i) / #(G на позиции i) ```bash cd ~/Projects/ffpe_damage python3 -m venv --system-site-packages .venv -.venv/bin/pip install pysam==0.24.0 # плюс pandas, matplotlib, scipy +.venv/bin/pip install pysam==0.24.0 openpyxl # плюс pandas, matplotlib, scipy, requests ``` Требования на входные данные: -- BAM отсортирован и проиндексирован (`.bai`); -- референсная FASTA проиндексирована (`.fai`). +- BAM отсортирован и проиндексирован (`.bai` — создаётся автоматически если отсутствует); +- референсная FASTA проиндексирована (`.fai` — создаётся автоматически). + +### Референс hg38 (UCSC, `chr`-префикс) + +```bash +wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz +gunzip hg38.fa.gz && samtools faidx hg38.fa +# md5 hg38.fa.gz: 1c9dcaddfa41027f17cd8f7a82c7293b +``` + +Для офлайн-аннотации дополнительно: + +```bash +mkdir -p references +wget -O references/gencode.v44.annotation.gtf.gz \ + https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.annotation.gtf.gz +# 48M, 252k транскриптов — используется для развёртки всех транскриптов без сети +``` ## 1. Анализ образца @@ -49,14 +66,26 @@ python3 -m venv --system-site-packages .venv --outdir FFPE_sample ``` -Строгий вариант: +Строгий вариант + многопоток (auto = `ядра-1`, cap 8): ```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 + --mapq 30 --baseq 20 --min-depth 30 --min-alt-count 5 --threads 7 ``` +Фильтрация повреждённых прямо по BAM с опциональным VCF: + +```bash +# удалить повреждённые из VCF: +.venv/bin/python ffpe_damage_v2.py --bam sample.bam --reference hg38.fa --outdir FFPE_sample --vcf input.vcf +# -> FFPE_sample/input.clean.vcf + input.ffpe_flagged.vcf +# оставить но пометить FILTER=FFPE: +.venv/bin/python ffpe_damage_v2.py ... --vcf input.vcf --keep-damaged +``` + +Настройка детекции повреждённых: `--filter-end 10` (расстояние от конца read) и `--filter-strand-p 0.05`. + ### Выходные файлы | Файл | Содержимое | @@ -72,7 +101,8 @@ python3 -m venv --system-site-packages .venv | `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.csv` | контекст каждой позиции: VAF, средняя позиция в read, strand bias, BQ, раскладка R1/R2 | +| `candidate_CtoT_clean.bed`, `candidate_GtoA_clean.bed` | те же BED без FFPE-подозрительных | +| `candidate_variants.csv` | контекст каждой позиции: VAF, средняя позиция в read, strand bias, BQ, раскладка R1/R2, `is_ffpe_suspect` | | `clean_variants.csv` / `ffpe_suspect_variants.csv` | разделение на чистые и подозрительные | | `CtoT_R1.png` … `GtoA_R2_3prime.png` | графики частоты по позиции в read (5′ и 3′) | @@ -80,7 +110,7 @@ python3 -m venv --system-site-packages .venv - **Классическая FFPE-сигнатура**: пик `C>T` у 5′-конца R1 и `G>A` у 3′-конца R2 — две цепи одного и того же дезаминирования цитозина. Проверяйте профили R1/R2 и по strand. - Если C>T/G>A повышены равномерно по всей длине read (без пика у концов) — это не классическое end-associated повреждение, а другой артефакт/ошибки. -- **Отличить настоящий вариант от дезаминирования** (`candidate_variants.csv`): настоящий вариант имеет ожидаемый VAF, распределён по read (средняя позиция ~ середина), встречается на обеих цепях (strand bias не значим) и не привязан к концам. Дезаминирование концентрируется у концов read и на одной цепи. +- **Отличить настоящий вариант от дезаминирования** (`candidate_variants.csv`): настоящий вариант имеет ожидаемый VAF, распределён по read (средняя позиция ~ середина), встречается на обеих цепях (strand bias не значим) и не привязан к концам. Дезаминирование концентрируется у концов read и на одной цепи. Критерий по умолчанию: `C>T/G>A` + `mean_pos ≤10` от любого конца → `is_ffpe_suspect`. ## 2. Классификация относительно контролей @@ -112,7 +142,36 @@ low : z < -k > Зачем сравнивать только с контролями? FFPE-повреждение зависит от протокола (UDG-обработка, end repair, длина фрагментов, панель). Одно и то же значение частоты может быть нормой для одного протокола и тяжёлым повреждением для другого. Поэтому сравнивать стоит только образцы, приготовленные одинаково, а классификация всегда относительная. -## 3. Синтетические данные для проверки +## 3. Аннотация чистых вариантов → Excel + +Из `clean_variants.csv` (0-based `pos`, `VAF = alt_count/depth`) — таблица 8 колонок: + +`Ген | HGVS (c.6713C>T p.(Pro2238Leu)) | Тип варианта и эффект | VAF | PAF | ACMG значимость | AMP уровень | Уровень онкогенности` + +```bash +# Быстро офлайн, 1 лучший транскрипт на вариант (по умолчанию, MANE > appris > canonical): +.venv/bin/python annotate_clean.py \ + --clean FFPE_sample/clean_variants.csv \ + --out FFPE_sample/clean.annotated.xlsx \ + --offline --min-vaf 0.10 --min-depth 20 + +# Все транскрипты (1 вариант → N строк, ~9.5 среднее): +.venv/bin/python annotate_clean.py ... --offline --min-vaf 0.10 --all-transcripts + +# Онлайн pan-cancer (OncoKB токен ~/.config/oncokb/token, Ensembl REST для HGVS): +.venv/bin/python annotate_clean.py --clean FFPE_sample/clean_variants.csv --out FFPE_sample/clean.annotated.xlsx --min-vaf 0.10 +# токен: echo "a6061cfc-..." > ~/.config/oncokb/token && chmod 600 ~/.config/oncokb/token +# другой tumor_type: --tumor-type "LUAD" +``` + +Детали: +- **Ген/HGVS/Эффект** — офлайн: `GENCODE v44` GTF (252k транскриптов, `chr`-префикс, `MANE_Select` > `appris_principal_1` > `Ensembl_canonical` > `basic` > `protein_coding`), `HGVS = ENST...:c.25234445C>T p.(?) (chr1:g.25234445C>T)`, `Тип = transcript_type` или `missense_variant`. Онлайн: `Ensembl REST /vep/homo_sapiens/hgvs` с тем же приоритетом. +- **VAF** — из `clean_variants.csv` как `65.02%`, `PAF` — пусто (по запросу). +- **ACMG** — `ClinVar` + `InterVar` placeholder (`Uncertain significance` без локальной БД, требует ручной курации для `Pathogenic`). +- **AMP/Онкогенность** — `OncoKB` (`pan-cancer` / `All Solid Tumors` по умолчанию, `LEVEL_1→Tier I` и т.д.) + `CIViC` fallback; в `--offline` пропускается (пусто). Большие списки (>500) пропускают OncoKB с предупреждением — фильтруйте по `VAF` перед онлайн. +- **Excel** — `openpyxl`, цвет по `ACMG`/`AMP`/`Onco`, автофильтр, заморозка, ширина колонок. При `VAF>10%` из 3501 → ~318 строк (1:1), с `--all-transcripts` → ~3018. + +## 4. Синтетические данные для проверки ```bash .venv/bin/python make_test_data.py --outdir test_data/damaged --damage-rate 0.5