ffpe_damage

Методологически корректная оценка дезаминирования в FFPE-образцах (формалин-фиксированная парафинизированная ткань) по BAM.

Инструмент не изобретает универсальный «FFPE damage score». Он выдаёт измеряемые показатели:

• все 12 типов mismatch;
• нормализованную частоту C>T и G>A на каждой позиции read;
• отдельно R1/R2;
• отдельно forward/reverse alignment;
• расстояния от 5′- и 3′-конца read;
• профиль частот (не просто количество mismatch);
• BED-треки для IGV;
• CSV для сравнения образцов между собой.

Классификация 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).

Установка

cd ~/Projects/ffpe_damage
python3 -m venv --system-site-packages .venv
.venv/bin/pip install pysam==0.24.0 openpyxl  # плюс pandas, matplotlib, scipy, requests

Требования на входные данные:

• BAM отсортирован и проиндексирован (.bai — создаётся автоматически если отсутствует);
• референсная FASTA проиндексирована (.fai — создаётся автоматически).

Референс hg38 (UCSC, chr-префикс)

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

Для офлайн-аннотации дополнительно:

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. Анализ образца

.venv/bin/python ffpe_damage_v2.py \
    --bam sample.bam \
    --reference hg38.fa \
    --outdir FFPE_sample

Строгий вариант + многопоток (auto = ядра-1, cap 8):

.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 --threads 7

Фильтрация повреждённых прямо по BAM с опциональным VCF:

# удалить повреждённые из 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.

Выходные файлы

Файл	Содержимое
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_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′)

Как читать результат

• Классическая 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 и на одной цепи. Критерий по умолчанию: C>T/G>A + mean_pos ≤10 от любого конца → is_ffpe_suspect.

2. Классификация относительно контролей

.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. Аннотация чистых вариантов → Excel

Из clean_variants.csv (0-based pos, VAF = alt_count/depth) — таблица 8 колонок:

Ген | HGVS (c.6713C>T p.(Pro2238Leu)) | Тип варианта и эффект | VAF | PAF | ACMG значимость | AMP уровень | Уровень онкогенности

# Быстро офлайн, 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. Синтетические данные для проверки

.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.
• Сравнивать нужно только образцы одного протокола подготовки библиотеки.