Format mówi, na którym etapie analizy jesteś

FASTA, FASTQ, BAM i VCF nie są czterema sposobami zapisu tego samego. Reprezentują inne obiekty: sekwencję referencyjną, surowe odczyty z jakością, odczyty uliniowione do referencji i wywołane warianty. Rozpoznanie formatu pozwala od razu ustalić, co już zrobiono z danymi i jakich informacji po drodze ubyło — a także jakim narzędziem wolno go otworzyć. Ten artykuł przechodzi przez formaty w kolejności potoku, pokazując przy każdym realne programy, którymi się go czyta, waliduje i przekształca, oraz pułapki, które cicho psują wyniki.

Cztery formaty to cztery etapy tego samego potoku. Surowe odczyty (FASTQ) mapuje się do referencji (FASTA), tworząc uliniowienia (BAM), z których wywołuje się warianty (VCF). Na każdym etapie coś dochodzi, a coś ubywa.
Cztery formaty to cztery etapy tego samego potoku. Surowe odczyty (FASTQ) mapuje się do referencji (FASTA), tworząc uliniowienia (BAM), z których wywołuje się warianty (VCF). Na każdym etapie coś dochodzi, a coś ubywa.

Większość tych formatów obsługuje jeden wspólny ekosystem narzędzi: htslib i zbudowane na nim samtools (SAM/BAM/CRAM), bcftools (VCF/BCF) oraz tabix/bgzip (kompresja i indeksowanie). Do przedziałów służy bedtools, do FASTA/FASTQ — seqkit i seqtk, do wizualizacji — IGV (Integrative Genomics Viewer) i przeglądarki UCSC/Ensembl. Znajomość tych kilku programów pokrywa większość codziennej pracy z plikami.

FASTA: identyfikator i sekwencja

Rekord FASTA zaczyna się wierszem > z identyfikatorem, po którym następuje sekwencja nukleotydowa lub białkowa. Jeden plik może przechowywać pojedynczy gen, cały zestaw chromosomów albo proteom. Nie zawiera jakości zasad ani informacji o pochodzeniu odczytu. To, czy A oznacza nukleotyd, czy aminokwas, wynika z kontekstu; równie ważna jest wersja referencji — współrzędna chr1:1000 ma sens dopiero wraz z nazwą zestawu (np. GRCh38) i wersją sekwencji. Do jednoznacznej identyfikacji treści przydaje się suma kontrolna (np. md5), a nie nazwa pliku.

W praktyce genomy referencyjne indeksuje się poleceniem samtools faidx (tworzy .fai), co pozwala błyskawicznie wyciąć dowolny region bez wczytywania całego chromosomu. Do manipulacji (filtrowanie, przeliczanie długości, wyciąganie podzbiorów, konwersje) służą seqkit i seqtk.

FASTQ: sekwencja plus jakość pomiaru

Każdy rekord FASTQ ma zwykle cztery wiersze: identyfikator rozpoczynający się @, sekwencję, separator + oraz ciąg znaków kodujących jakość poszczególnych zasad. Wartość Phred jest logarytmiczną oceną prawdopodobieństwa błędu, nie gwarancją prawdziwości zasady. Warto zapamiętać kilka progów, bo wracają w każdym raporcie: Q20 ≈ 1% błędu, Q30 ≈ 0,1%, Q40 ≈ 0,01% — każde 10 punktów Phred to dziesięciokrotnie mniejsze ryzyko. Q30 bywa progiem „dobrej" jakości, ale właściwy próg zależy od zastosowania.

Historia formatu pozostawiła różne kodowania jakości: współczesny standard to Phred+33 (Sanger / Illumina 1.8+), ale w starych danych spotyka się Phred+64 (dawne Illumina/Solexa) — pomylenie ich zawyża lub zaniża jakość, dlatego przy archiwalnych plikach trzeba sprawdzić konwencję.2 Dane paired-end występują w dwóch powiązanych plikach (_R1/_R2), których kolejność i identyfikatory par muszą się zgadzać.

Przed analizą wykonuje się kontrolę jakości: FastQC (rozkład jakości, zawartość adapterów, duplikacja, skład zasad), zbiorczo MultiQC, a przycinanie adapterów i słabych końców — fastp, Trimmomatic lub cutadapt. FASTQ nie mówi jeszcze, gdzie odczyt pasuje w genomie — to zadanie mapowania.

SAM, BAM i CRAM: odczyt względem referencji

Uliniowienia zapisuje SAM (tekstowo) i BAM (binarnie, skompresowane). Powstają z FASTQ przez mapowanie programami takimi jak BWA-MEM i minimap2 (długie odczyty), Bowtie2, a dla RNA — STAR i HISAT2. Nagłówek (@) opisuje m.in. referencję i użyte programy, a każdy rekord ma 11 obowiązkowych pól: nazwę odczytu, FLAG, chromosom i pozycję, MAPQ, CIGAR, informacje o parze, sekwencję i jakości, plus pola opcjonalne.

Dwa pola bywają źródłem nieporozumień. FLAG to pole bitowe — pojedyncza liczba kodująca zestaw cech naraz (sparowany, odwrócony, niezmapowany, wtórny, duplikat itd.); rozkłada się ją na bity, nie czyta „wprost". CIGAR opisuje, jak odczyt pasuje do referencji, literami operacji: M dopasowanie/niedopasowanie, I insercja, D delecja, S przycięcie miękkie, N przeskok (istotny dla RNA na granicach egzonów). Jakość mapowania (MAPQ) to nie to samo co jakość zasad z FASTQ — dotyczy niepewności przypisania odczytu do miejsca; region powtarzalny może mieć świetne zasady, a niejednoznaczną pozycję (niskie MAPQ).

Operacje wykonuje samtools: sort (po współrzędnych), index (tworzy .bai/.csi do szybkiego pobierania regionu), view (filtrowanie po FLAG/MAPQ/regionie), flagstat i markdup. Do usuwania duplikatów PCR używa się też Picard MarkDuplicates. CRAM to oszczędna alternatywa BAM, kompresująca względem referencji (.crai do indeksu) — potrafi znacząco zmniejszyć pliki, ale bez właściwej referencji jest trudniejszy do odtworzenia, więc referencję trzeba przechowywać razem z danymi.

VCF i BCF: hipotezy o wariantach

VCF ma nagłówek z metadanymi oraz osiem stałych kolumn (CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO), a przy genotypach — kolumnę FORMAT i po jednej kolumnie na próbkę. Powstaje z BAM przez wywoływanie wariantów: GATK (HaplotypeCaller), bcftools call, DeepVariant lub FreeBayes. Pola QUAL, FILTER i INFO interpretuje się zgodnie z pipeline'em, który je wytworzył — PASS oznacza tylko przejście zdefiniowanych filtrów, nie istotność kliniczną ani przyczynowość biologiczną.

Wariant jest opisany względem konkretnej referencji, a insercje i delecje można zapisać na kilka równoważnych sposobów, dlatego przed porównaniem zbiorów wykonuje się normalizację (bcftools norm — lewostronne wyrównanie i rozbicie wariantów wieloallelicznych) oraz sprawdza zgodność allelu referencyjnego i zestawu genomu. Odmianą jest gVCF (GATK) — zapisuje też pozycje bez wariantu, co umożliwia łączne genotypowanie wielu próbek. Binarnym odpowiednikiem VCF jest BCF (szybszy w obróbce). Narzędzia: bcftools i vcftools; specyfikacje SAM/BAM/CRAM oraz VCF/BCF utrzymuje grupa hts-specs w ramach GA4GH.1

Współrzędne: 0-based kontra 1-based

BED jest 0-indeksowany i półotwarty [start, koniec), a GFF/GTF/VCF są 1-indeksowane i domknięte [start, koniec]. Ten sam trzyzasadowy odcinek ma w każdym formacie inne liczby — pomieszanie ich daje przesunięcie o jedną zasadę.
BED jest 0-indeksowany i półotwarty [start, koniec), a GFF/GTF/VCF są 1-indeksowane i domknięte [start, koniec]. Ten sam trzyzasadowy odcinek ma w każdym formacie inne liczby — pomieszanie ich daje przesunięcie o jedną zasadę.

Jedną z najcichszych pułapek jest konwencja liczenia pozycji. BED jest 0-indeksowany i półotwarty: pierwsza zasada ma współrzędną 0, a przedział [start, koniec) nie wlicza końca (długość = koniec − start). GFF/GTF i VCF są 1-indeksowane i domknięte: pierwsza zasada to 1, oba końce wliczone (długość = koniec − start + 1). Ten sam odcinek ma więc w BED i w VCF inne liczby. Narzędzia znają konwencję swojego formatu, ale ręczne mieszanie współrzędnych między formatami (albo między wersjami genomu) daje przesunięcie o jedną zasadę — wynik wygląda poprawnie, a jest błędny.

Adnotacje i cechy genomu: BED, GFF/GTF i ścieżki pokrycia

Obok formatów sekwencyjnych stale pojawiają się pliki opisujące położenie cech. BED to prosty format przedziałów (chromosom, początek, koniec, opcjonalne kolumny) do regionów zainteresowania, pików czy masek. GFF3/GTF opisują geny, egzony i transkrypty z atrybutami (GTF to starszy wariant używany m.in. przez narzędzia RNA-seq; GFF3 jest ogólniejszy). Operacje na przedziałach (przecięcia, pokrycia, sąsiedztwa) wykonuje bedtools. Ciągłe ścieżki pokrycia/sygnału dla przeglądarek zapisuje się binarnie jako BigWig i bigBed (UCSC). Wszystkie one są sensowne wyłącznie względem konkretnej wersji referencji i spójnego nazewnictwa chromosomów (chr1 vs 1).

Poza „wielką czwórką": dane jednokomórkowe, surowy sygnał i inne

Nowoczesne analizy generują formaty spoza klasycznego potoku:

  • Dane jednokomórkowe: macierze zliczeń zapisuje się jako Matrix Market (.mtx) (rzadkie macierze z Cell Ranger), a przetworzone obiekty jako AnnData (.h5ad) — oparty na HDF5 kontener używany przez scanpy — oraz Loom. HDF5 to ogólny, binarny format do dużych, wielowymiarowych danych.
  • Surowy sygnał sekwenatorów: nanopore zapisuje sygnał prądowy jako FAST5 i nowszy POD5 (basecalling zamienia go na FASTQ). Illumina generuje surowe BCL, konwertowane do FASTQ przez bcl2fastq / BCL Convert.
  • Proteomika i inne dziedziny: widma mas w mzML, a np. drzewa filogenetyczne w Newick. Każda dziedzina ma własny standard i własne walidatory.

Nagłówek jest częścią danych

Najczęstszy błąd początkującego to otworzyć plik jak tabelę i pominąć metadane. Nagłówek może zawierać wersję referencji, kolejność kontigów, definicje filtrów i historię narzędzi; jego usunięcie sprawia, że rekordy pozostają czytelne, lecz przestają być jednoznaczne. Minimalna kontrola każdego pliku obejmuje:

  • format i możliwość odczytu przez walidator (np. samtools quickcheck, bcftools view, gffread);
  • liczbę rekordów oraz oczekiwane próbki;
  • wersję referencji i nazewnictwo chromosomów;
  • sumę kontrolną po przesłaniu;
  • zgodność nagłówka z treścią;
  • pochodzenie i etap pipeline'u.

Format nie zastępuje biologicznych metadanych: FASTQ może być technicznie poprawny, a nadal nie wiadomo, z jakiej tkanki i od którego dawcy pochodzi. Gdzie szukać obu typów informacji, wyjaśnia artykuł Bazy danych dla biologa, a jak zautomatyzować te kontrole w powtarzalnym potoku — Workflow i kontrola jakości danych.