Pokrycie, mapowalność, duplikaty i bias referencji

Streszczenie

„Pokrycie 30×” brzmi jak kompletna informacja o jakości sekwencjonowania, ale jest tylko średnią. Ten sam wynik może powstać z trzydziestu niezależnych cząsteczek, setek kopii PCR jednej molekuły albo nadmiaru reads w łatwych regionach przy braku danych w trudnych. Wariant można wykryć tylko tam, gdzie próbka, biblioteka, instrument i mapowanie dostarczyły informacyjnego dowodu.

Artykuł rozdziela depth, breadth, physical i molecular coverage, duplikaty, mapowalność, mapping quality, callable region, GC bias, allele dropout, paralogi, pseudogeny oraz reference bias. Pokazuje rachunki, diagnostykę i sposób raportowania. Celem nie jest osiągnięcie jednej wysokiej liczby, lecz ustalenie, które hipotezy biologiczne dane rzeczywiście rozstrzygają.

Pokrycie ma kilka znaczeń

Najprostsze pokrycie bazowe pozycji:

depth = liczba reads obejmujących pozycję

Jeżeli 30 reads przecina bazę, depth wynosi 30×. To nie mówi jeszcze:

• czy reads są niezależne;
• czy mapują jednoznacznie;
• czy pochodzą z obu alleli;
• czy mają właściwą orientację;
• czy zawierają wysokiej jakości bazę;
• czy nie są artefaktem PCR;
• czy obejmują cały wariant.

Mean, median i rozkład

Średnie pokrycie:

mean coverage = suma baz w alignments / rozmiar analizowanego regionu

Może je zawyżyć kilka regionów o bardzo dużej głębokości. Mediana jest odporniejsza, ale też nie pokazuje dolnego ogona. Dla diagnostyki ważny bywa udział celu pokryty co najmniej 10×, 20×, 30× lub progiem wynikającym z walidacji.

Raport powinien zawierać:

• średnią;
• medianę;
• percentyle;
• histogram;
• procent celu powyżej progów;
• listę regionów poniżej progu;
• coverage po filtrach.

Breadth of coverage

Breadth odpowiada na pytanie, jaki udział celu ma co najmniej określoną liczbę reads.

Przykład:

• panel ma 1 000 000 zasad;
• 980 000 ma ≥20×;
• breadth at 20× = 98%.

Pozostałe 2% może obejmować geny nieważne albo jedyny krytyczny ekson. Globalny procent nie zastępuje listy miejsc.

Model Lander–Waterman i jego granice

W idealnym losowym shotgun sequencing liczba reads pokrywających pozycję może być modelowana rozkładem Poissona. Prawdopodobieństwo braku pokrycia przy średnim C:

P(0) = e^-C

Dla 5×:

e^-5 ≈ 0,0067

czyli około 0,67% baz byłoby niepokrytych w idealnym modelu. Klasyczny model Lander–Waterman połączył głębokość, liczbę fragmentów i oczekiwaną ciągłość mapowania genomu.¹

Rzeczywiste biblioteki naruszają założenia:

• fragmentacja nie jest idealnie losowa;
• GC wpływa na odzysk;
• PCR faworyzuje część cząsteczek;
• capture ma nierówne sondy;
• mapowalność zależy od sekwencji;
• ploidalność i CNV zmieniają oczekiwaną głębokość;
• filtry odrzucają część reads.

Dlatego większe średnie coverage nie usuwa systematycznych dziur.

Sequence coverage i physical coverage

Sequence coverage liczy bazy faktycznie odczytane.

Physical coverage liczy fragmenty, których znane końce obejmują region, także przestrzeń między paired-end reads.

Para 2 × 150 bp z insertem 500 bp daje:

• 300 bp sequence coverage;
• około 500 bp physical coverage.

Physical coverage pomaga wykrywać warianty strukturalne i łączyć loci, ale nie dostarcza bazy w nieodczytanym środku.

Spanning coverage

Dla powtórzenia lub breakpointu liczy się liczba cząsteczek, które:

• zaczynają się w unikalnej kotwicy przed elementem;
• przechodzą przez całość;
• kończą w unikalnej kotwicy za elementem.

Read pokrywający środek powtórzenia bez kotwic nie rozstrzyga jego lokalizacji. Średnie 30× może dać zero spanning reads dla elementu dłuższego niż insert.

Haplotype coverage

W heterozygotycznym locus oczekuje się dowodu dla obu alleli. Przy 30 niezależnych cząsteczkach liczba alternatywnych reads losowo się waha. Allele imbalance może jednak wynikać z:

• małej liczby;
• reference bias;
• allele dropout;
• CNV;
• mosaicism;
• strand bias;
• mapping różnicowego.

Coverage trzeba analizować per haplotyp, jeśli wniosek dotyczy cis/trans lub allele-specific expression.

Pokrycie cząsteczkowe

W PCR wiele reads może pochodzić z jednej pierwotnej cząsteczki. Sequence depth rośnie, ale liczba niezależnych obserwacji nie.

Molecular coverage liczy pierwotne molekuły, zwykle z pomocą UMI.

Przykład:

• 100 reads obejmuje locus;
• po grupowaniu UMI powstaje 12 rodzin;
• surowy depth = 100×;
• molecular depth = 12×.

Granica detekcji rzadkiego wariantu zależy bardziej od 12 niż od 100.

UMI ma skończoną przestrzeń

UMI długości 8 ma teoretycznie:

4^8 = 65 536

kodów. Nie wszystkie są równie prawdopodobne, a dwie molekuły mogą dostać ten sam kod. Błąd sekwencjonowania może rozdzielić jedną rodzinę na dwie.

UMI-tools formalizuje grupowanie i deduplikację z wykorzystaniem struktury sieci kodów, zamiast traktować każdy różny UMI jako niezależny.⁸

Duplex molecular coverage

Jeśli oznaczono obie nici pierwotnego dupleksu, konsensus duplex wymaga zgodnego wariantu na komplementarnych rodzinach. To silnie ogranicza błędy uszkodzenia i PCR, ale zmniejsza liczbę cząsteczek spełniających kryterium.

Raport rzadkiego wariantu powinien podać:

• raw reads;
• rodziny UMI;
• duplex families;
• liczebność alternatywy;
• liczebność referencji;
• próg jakości rodziny.

Duplikat optyczny, PCR i biologiczny

Nie każda identyczna para współrzędnych ma tę samą przyczynę.

Optical duplicate powstaje, gdy jeden klaster jest rozpoznany jako kilka obiektów.

PCR duplicate pochodzi z amplifikacji tej samej cząsteczki.

Natural collision występuje, gdy niezależne fragmenty mają te same końce, szczególnie w panelu celowanym lub przy bardzo dużej głębokości.

Deduplikacja po współrzędnych

W losowo fragmentowanej bibliotece WGS para pozycji start/koniec jest przybliżeniem tożsamości molekuły. W ampliconie wszystkie prawidłowe molekuły mają te same końce. Usunięcie „duplikatów” po współrzędnych skasowałoby biologiczny dowód.

Strategia zależy od biblioteki:

• shotgun bez UMI — współrzędne;
• capture z UMI — UMI + współrzędne;
• amplicon bez UMI — nie można wiarygodnie rozpoznać molekuł;
• duplex — UMI obu nici;
• long reads — inne źródła powtórzeń.

Mark czy remove

Marking zachowuje reads w BAM z flagą. Caller może je ignorować lub analizować osobno. Fizyczne usunięcie utrudnia audyt.

Duplikaty są również metryką biblioteki. Wysoki udział może oznaczać:

• mało materiału wejściowego;
• nadmierny PCR;
• zbyt głębokie sekwencjonowanie;
• mały target;
• silne wzbogacenie;
• bottleneck capture.

Nie należy interpretować procentu bez przewidywanej złożoności.

Library complexity

Złożoność to liczba różnych molekuł dostępnych do sekwencjonowania. Kolejne reads coraz częściej trafiają do już obserwowanych fragmentów.

Krzywa complexity pokazuje:

• unikalne reads przy obecnej głębokości;
• oczekiwany przyrost przy dalszym sekwencjonowaniu;
• punkt nasycenia;
• różnice między bibliotekami.

Więcej sekwencjonowania może nie pomóc

Jeśli biblioteka ma mało molekuł, podwojenie liczby reads głównie zwiększy duplikaty. Rozwiązaniem może być:

• nowa ekstrakcja;
• większy input;
• mniej PCR;
• lepszy capture;
• inny protokół;
• więcej niezależnych bibliotek.

Dokupienie lane nie odtwarza utraconej różnorodności.

Mapowanie jest klasyfikacją hipotez

Aligner nie „wie”, skąd pochodzi read. Porównuje możliwe miejsca i wybiera dopasowanie według modelu mismatch, gap, seeds i powtórzeń.

Wynik zawiera:

• pozycję;
• CIGAR;
• orientation;
• edit distance;
• alignment score;
• mapping quality;
• secondary/supplementary alignments.

Base quality i mapping quality

Base quality estymuje prawdopodobieństwo błędnej litery.

Mapping quality estymuje niepewność lokalizacji alignera.

Read może mieć Q40 dla każdej bazy i MAPQ 0, jeśli pasuje identycznie do dziesięciu miejsc. Może też mapować jednoznacznie, ale zawierać słabą bazę.

MAPQ nie jest uniwersalną skalą między programami

Nominalnie MAPQ ma formę Phred:

MAPQ = -10 log10 P(złe mapowanie)

Implementacje różnią się sposobem estymacji i wartością maksymalną. Próg MAPQ 20 w dwóch alignerach nie musi wybierać identycznych reads.

Raport powinien podać aligner, wersję, parametry i dystrybucję MAPQ.

CIGAR opisuje relację, nie prawdę biologiczną

Operacje CIGAR rozdzielają:

• M — pozycje alignment;
• I — insercję względem referencji;
• D — delecję;
• N — pominięcie, typowo intron;
• S — soft clipping;
• H — hard clipping.

Ten sam wariant może mieć kilka równoważnych alignments w powtórzeniu. Soft-clipped read może wskazywać breakpoint, adapter, słabą jakość albo chimerę.

Secondary i supplementary

Secondary alignment jest alternatywną lokalizacją tego samego readu.

Supplementary alignment reprezentuje inną część readu, na przykład po drugiej stronie wariantu strukturalnego.

Odrzucenie wszystkich dodatkowych alignments usuwa dowód split-read. Zliczenie ich jak niezależnych reads zawyża coverage.

Pair geometry

Paired-end reads dostarczają:

• orientację;
• odległość;
• chromosome;
• insert distribution.

Discordant pair może wskazywać:

• deletion;
• inversion;
• translocation;
• chimeryczną bibliotekę;
• mapping error.

Dowód jest statystyczny względem rozkładu normalnych par w tej bibliotece.

Mapowalność jest właściwością sekwencji i odczytu

Region jest mapowalny, jeśli read o określonej długości i błędzie ma jednoznaczne miejsce w referencji.

Mapowalność zależy od:

• długości readu;
• paired-end insert;
• profilu błędu;
• referencji;
• dopuszczalnych mismatch;
• haplotypu;
• algorytmu.

Nie jest stałą etykietą locus.

Genomic dark matter

Krótkie reads słabo oceniają:

• duplikacje segmentowe;
• centromery;
• telomery;
• tandem repeats;
• pseudogeny;
• rodziny genów;
• HLA/KIR;
• SMN1/SMN2;
• CYP2D6 i inne paralogi.

Lee i Schatz opisali „genomową ciemną materię” jako sekwencje trudne do analizy krótkimi reads ze względu na powtórzenia i brak jednoznacznego mapowania.³

Dłuższy read pomaga tylko z kotwicą

Read może być długi, ale w całości leżeć w identycznej duplikacji. Informację daje fragment sięgający unikalnej różnicy lub granicy.

Minimalna informacyjna długość:

powtórzenie + unikalna kotwica po jednej lub obu stronach

To dlatego read N50 nie zastępuje analizy konkretnego locus.

Paralog i pseudogen mogą tworzyć fałszywy wariant

Jeśli paralog ma stałą różnicę względem genu docelowego, reads źle przypisane do celu wyglądają jak heterozygotyczny wariant.

Sygnały ostrzegawcze:

• VAF około udziału paralogów;
• wiele „wariantów” odpowiadających paralogowi;
• niskie MAPQ;
• depth większy od oczekiwanego;
• imbalance strand;
• brak potwierdzenia primerem swoistym;
• zależność od alignera.

Wysoko homologiczne geny wymagają dedykowanych strategii: long-range PCR, swoistych kotwic, long reads, copy-number model i bazy paralog-specific variants. Przegląd trudnych klinicznie genów pokazał, że standardowe krótkie reads pozostawiają ważne regiony poza wiarygodną analizą.⁴

Gene conversion

Paralogi mogą wymieniać fragmenty przez gene conversion. Różnica podobna do paralogu nie zawsze jest artefaktem mapowania; może być prawdziwym zdarzeniem. Rozstrzygnięcie wymaga cząsteczki łączącej zmianę z unikalną kotwicą genu.

Reference bias

Read z allelem zgodnym z referencją ma mniej mismatch niż read alternatywny. Aligner może:

• łatwiej go zmapować;
• nadać wyższy MAPQ;
• odrzucić read alternatywny;
• przesunąć alignment indela;
• wybrać inne miejsce.

W rezultacie VAF jest zaniżony, a allele-specific expression lub methylation faworyzuje referencję. Efekt read mapping bias w analizie allelospecyficznej był wykazany już dla krótkich reads RNA-seq.⁵

Indel reference bias

Read z insercją lub delecją wymaga gap. W powtórzeniu można umieścić gap w kilku równoważnych miejscach. Normalizacja reprezentacji pomaga porównać VCF, ale nie odzyskuje reads odrzuconych przy mapowaniu.

Bias zależy od:

• długości indela;
• kontekstu;
• read length;
• local assembly;
• scoring alignera;
• jakości końców.

Strategie ograniczania

• spersonalizowana referencja;
• graph/pangenome reference;
• remapping po podmianie allelu;
• WASP-like filtering;
• local assembly;
• dłuższe reads;
• oba haplotypy;
• benchmark na wariantach heterozygotycznych.

Każda ma koszt. Odrzucenie reads z biasem może zmniejszyć moc.

Referencja jest częścią pomiaru

Ta sama próbka mapowana do dwóch wersji genomu może dać:

• inne współrzędne;
• inne paralogi;
• naprawione błędy;
• nowe contigi;
• odmienny decoy;
• inne alt loci;
• inny callable region.

Raport „GRCh38” może być niewystarczający, jeśli pipeline używa innego zestawu decoy, HLA contigs lub mask.

Alt contigs i decoys

Alternate loci reprezentują odmienne haplotypy trudnych regionów. Decoy sequences przechwytują reads, które inaczej błędnie mapowałyby do głównych chromosomów.

Dodanie decoy może:

• zmniejszyć fałszywe warianty;
• obniżyć depth w locus;
• zmienić MAPQ;
• ujawnić multi-mapping.

Spadek liczby wariantów po lepszej referencji może oznaczać poprawę, nie utratę czułości.

Pangenom jako zbiór ścieżek

Graph reference pozwala readowi dopasować się do ścieżki zawierającej allel alternatywny. Może zmniejszać reference bias i reprezentować duże insercje nieobecne w liniowej referencji.

Nowe pytania:

• która ścieżka jest współrzędną raportową;
• jak normalizować wariant;
• jak porównać pipeline liniowy i grafowy;
• jak oceniać mapowalność w grafie;
• jak liftować wynik;
• jakie populacje reprezentuje graf.

Pangenom nie usuwa biasu automatycznie. Niewystarczająca reprezentacja może faworyzować częstsze ścieżki.

GC bias

Regiony bardzo bogate lub ubogie w GC mogą mieć mniejsze pokrycie z powodu:

• ekstrakcji;
• fragmentacji;
• PCR;
• hybrydyzacji capture;
• klastrów;
• basecallingu;
• mapowania.

Ross i współpracownicy porównali platformy i pokazali, że bias sekwencyjny ma charakter zależny od technologii oraz biblioteki.²

Krzywa GC

Dzieli się genom na okna, liczy GC oraz normalized coverage. Idealna krzywa byłaby płaska. Rzeczywista ma spadki na krańcach i może różnić się między próbkami.

Metryki:

• normalized coverage per GC bin;
• dropout;
• median absolute deviation;
• fold-80 base penalty;
• coverage uniformity.

Korekcja GC w CNV modelu nie czyni regionu wiarygodnym dla SNV, jeśli fizycznie brak reads.

Capture bias

W panelu i eksomie sonda musi hybrydyzować do celu. Wariant pod sondą może osłabić wiązanie jednego allelu. Powstaje capture allele dropout.

Inne źródła:

• gęstość sond;
• repeats;
• GC;
• fragment length;
• pseudogen;
• niedokładna granica eksonu;
• lot panelu;
• off-target hybridization.

On-target nie równa się kompletność

Biblioteka może mieć 90% on-target i nadal nie pokrywać jednego eksonu. On-target mówi o efektywności wykorzystania reads, breadth — o kompletności.

Projekt panelu powinien zawierać:

• padding wokół eksonów;
• alternatywne sondy;
• regiony pseudogenowe;
• znane polimorfizmy;
• kontrolę aktualnej adnotacji transkryptów;
• listę regionów niewalidowalnych.

PCR allele dropout

Wariant pod starterem może zmniejszyć amplifikację jednego allelu. Heterozygota wygląda wtedy jak homozygota referencyjna lub alternatywna.

Ryzyko rośnie przy:

• amplicon panels;
• Sanger;
• long-range PCR;
• niskim input;
• zdegradowanym DNA;
• dużym indelu;
• zmianie copy number.

Potwierdzenie tym samym starterem powtarza błąd. Ortogonalny test powinien użyć innych miejsc wiązania lub metody bez PCR.

Jakość materiału zmienia użyteczne coverage

FFPE

Formalina tworzy:

• fragmentację;
• crosslinks;
• deaminację C;
• uszkodzone końce;
• nierówną amplifikację.

Wysoki depth może zawierać wiele kopii krótkich uszkodzonych molekuł. UMI pomaga rozpoznać rodziny, a matched normal i filtry kontekstowe ograniczają artefakty.

W raporcie FFPE potrzebne są:

• insert size;
• molecular depth;
• deamination profile;
• duplication;
• on-target;
• time/age blocku;
• tumor purity;
• limit detekcji.

cfDNA

Cell-free DNA ma naturalny rozkład fragmentów związany z nukleosomami. Cząsteczki guza mogą stanowić mały udział.

Coverage musi być liczone po:

• collapse UMI;
• kontroli długości;
• usunięciu kontaminacji leukocytarnej;
• rozpoznaniu clonal hematopoiesis;
• duplex consensus.

Wariant w osoczu nie zawsze pochodzi z guza. Matched white blood cells pomagają rozdzielić źródło.

Ancient i environmental DNA

Krótkie uszkodzone fragmenty mają charakterystyczne zmiany na końcach. Mapowanie do odległej referencji zwiększa bias. Kontaminacja współczesnym DNA może mieć lepsze mapowanie i dominować mimo małej masy.

Authenticity wymaga profilu uszkodzeń, długości, blanków i niezależnych bibliotek.

Coverage mitochondrialne nie zachowuje się jak jądrowe

mtDNA ma wiele kopii per komórka. Depth może być tysiące razy większy niż genomu jądrowego.

Problemy:

• NUMTs;
• heteroplasmy;
• tissue-specific copy number;
• bottleneck;
• strand bias;
• circular reference boundary;
• PCR amplification.

Read podobny do mtDNA może pochodzić z NUMT. Długi fragment lub paired-end kotwica w jądrze pomaga rozstrzygnąć.

Heteroplasmy

VAF mitochondrialny nie odpowiada diploidalnemu genotypowi. Heteroplasmy 5% wymaga:

• molecular depth;
• tła błędu;
• niezależnych nici;
• właściwej tkanki;
• kontroli NUMT;
• powtarzalności.

Wysoki surowy depth bez rodzin nie wystarcza.

Coverage RNA-seq jest zależne od ekspresji

Brak reads w eksonie może znaczyć:

• gen nie jest eksprymowany;
• izoforma go pomija;
• biblioteka ma bias 3′;
• RNA jest zdegradowane;
• mapowanie splice zawiodło;
• ekson jest homologiczny.

RNA-seq nie jest pełnym testem wariantów genomowych. Może wspierać:

• ekspresję allelu;
• splice;
• fuzję;
• RNA editing;
• nonsense-mediated decay.

Allele-specific expression

Heterozygotyczny wariant DNA pozwala zliczyć allele w RNA. Reference bias, mapping i mała liczba reads mogą udawać imbalance.

Projekt wymaga:

• genotypu DNA;
• phasing, gdy analizuje się wiele loci;
• korekty mapowania;
• count per donor;
• modelu overdispersion;
• kontroli imprintingu i X.

Gene coverage nie jest równomierne

Poly(A) selection, degradacja i fragmentacja tworzą bias 3′ lub 5′. Junction reads nie pokrywają intronu, ale potwierdzają splicing.

Metryki:

• gene body coverage;
• junction saturation;
• exonic/intronic/intergenic fraction;
• strandedness;
• duplication;
• rRNA;
• insert size.

Long-read coverage ma własne pułapki

Długi read zwiększa kontekst, ale:

• ma wyższy błąd określonego typu;
• length distribution jest szeroki;
• HMW extraction selekcjonuje;
• ultradługie molekuły są nieliczne;
• aligner może splitować;
• polishing może mieszać allele.

Read depth i molecule depth są bliskie, ale nie identyczne

Biblioteka natywna bez PCR ogranicza klasyczne duplikaty. Jednak:

• re-sequencing tej samej cząsteczki;
• consensus reads;
• duplex;
• read splitting;
• chimeras,

wymagają określenia jednostki dowodu.

Coverage długości

Dla SV raportuje się:

• reads spanning;
• reads split;
• breakpoint-support;
• reference-support;
• phase-support;
• length distribution around locus;
• both-strand evidence.

Dziesięć reads kończących się wewnątrz powtórzenia nie zastępuje dwóch pełnych spanning molecules.

Prawdopodobieństwo obserwacji allelu

Przy n niezależnych molekuł i prawdziwej frakcji p liczba alternatywnych obserwacji ma w prostym modelu rozkład dwumianowy.

Prawdopodobieństwo zobaczenia zera:

P(X=0) = (1-p)^n

Dla heterozygoty p=0,5 i n=10:

P(0) = 0,5^10 ≈ 0,00098

Dla mozaiki 1% i 100 molekuł:

P(0) = 0,99^100 ≈ 0,366

Czyli nawet idealny assay ma około 36,6% szans nie zobaczyć żadnej cząsteczki alternatywnej.

Ile molekuł dla 95% szans co najmniej jednej obserwacji

Szukamy n, dla którego:

1 - (1-p)^n ≥ 0,95

Dla p=0,01:

n ≥ ln(0,05) / ln(0,99) ≈ 299

To tylko wykrycie co najmniej jednej cząsteczki. Jedna obserwacja nie odróżnia wariantu od błędu. Wiarygodna detekcja wymaga większej liczby, niskiego tła i kryterium obu nici.

Overdispersion

Rzeczywiste dane nie są idealnie dwumianowe z powodu:

• allele bias;
• PCR;
• batch;
• różnej wydajności molekuł;
• lokalnych artefaktów.

Walidacja empiryczna jest ważniejsza niż czysty rachunek.

Strand i orientation bias

Prawdziwy wariant powinien mieć dowód zgodny z geometrią biblioteki. Jeśli alternatywa występuje tylko:

• w forward reads;
• w R1;
• na końcu reads;
• na jednej tile;
• w jednym kierunku insertu,

może być artefaktem.

Nie każdy wariant ma idealną równowagę przy małej liczbie. Test wymaga modelu i liczebności, nie prostego progu 50:50.

Read position bias

Błędy rosną na końcach reads. Indel może być widoczny głównie przy soft-clipping. Rozkład pozycji alternatywy względem referencji pomaga odróżnić sygnał.

Copy number zmienia oczekiwane coverage

Deletion obniża depth, duplication podnosi. Jednak coverage zależy też od GC, capture i mapowalności. CNV caller porównuje:

• biny w próbce;
• kohortę referencyjną;
• allelic balance;
• segmentation;
• czasem split/discordant reads.

Normalizacja kohortowa

Próbki referencyjne powinny mieć:

• ten sam panel;
• podobny protokół;
• podobny input;
• ten sam instrument lub walidowaną zgodność;
• podobną jakość;
• brak dużej zmiany biologicznej w większości celu.

Batch może udawać CNV. Z kolei prawdziwy częsty CNV w kohorcie może zostać znormalizowany.

Mosaicism wymaga molekuł

W diploidalnym locus heterozygotyczny wariant germinalny ma oczekiwany VAF około 50%, z losowym i systematycznym odchyleniem. Mozaika 1% oznacza około 1 alternatywny allel na 100 odczytanych kopii, ale granica zależy od liczby molekuł i tła błędu.

Przy 1000 reads z 20 rodzin UMI nie ma 1000 niezależnych szans. Przy 1000 niezależnych molekuł można dopiero sensownie modelować 1%, jeśli metoda ma tło poniżej tego poziomu.

Limit of blank, detection i quantification

• LoB — sygnał obserwowany w blankach bez wariantu;
• LoD — najmniejszy poziom wiarygodnie odróżniany od blanku;
• LoQ — poziom, który można oszacować z określoną precyzją.

LoD nie powinien być deklarowany z jednego dilution. Wymaga wielu próbek, dni, operatorów, lotów i kontekstów.

Kontaminacja zmienia VAF i coverage

Domieszka innej próbki:

• rozcieńcza warianty;
• tworzy niskie VAF obcych alleli;
• zmienia heterozygotyczność;
• może udawać mozaikę;
• zaburza CNV.

Wskaźniki:

• fingerprint;
• allele balance;
• płeć chromosomalna;
• mitochondrial haplotype;
• contamination estimate;
• index hopping;
• blank.

Tumor purity

VAF w guzie zależy od:

• purity;
• local copy number;
• clonality;
• ploidy;
• normal contamination;
• loss of heterozygosity.

Somatyczny heterozygotyczny wariant klonalny nie musi mieć VAF 50%. Interpretacja bez CN i purity jest uproszczona.

Lejek filtrów trzeba zachować

Kolejne etapy odrzucają reads:

1. demultipleksowanie;
2. adapter/quality;
3. mapping;
4. MAPQ;
5. duplicate;
6. target;
7. base quality;
8. overlap clipping;
9. caller filters.

Raport powinien pokazać liczbę przed i po. Inaczej nie wiadomo, czy niski depth powstał w laboratorium czy analizie.

Overlap clipping

Paired-end reads krótszego insertu mogą nakładać się. Obie obserwują tę samą bazę tej samej cząsteczki. Zliczenie obu jako niezależnych zawyża depth i może podwoić błąd.

Caller lub preprocessing powinien rozpoznać overlap.

Callable region jest zależny od klasy wariantu

Pozycja callable dla SNV nie musi być callable dla:

• indela 20 bp;
• repeat expansion;
• CNV;
• inversion;
• mosaic variant;
• phasing;
• methylation.

Definicja powinna zawierać:

• klasę;
• zakres wielkości;
• minimalny depth;
• molecular depth;
• MAPQ;
• base quality;
• strand/orientation;
• mapowalność;
• kontekst;
• walidowany pipeline.

No-call jest uczciwym wynikiem

Brak wariantu w VCF może oznaczać:

• homozygotę referencyjną;
• niewystarczające dane;
• region wyłączony;
• filtr;
• brak modelu dla klasy wariantu.

Raport powinien rozróżniać negative i not assessed.

Trudne geny i benchmark

Benchmark na łatwych regionach zawyża wrażenie kompletności. Genome in a Bottle tworzy benchmark sets dla trudnych genów i rozszerzonych regionów, pozwalając ocenić pipeline tam, gdzie mapowanie oraz reprezentacja są problemem.⁶

Truth set też ma granice

Benchmark może:

• wykluczać część genomu;
• obejmować tylko określone klasy;
• zależeć od referencji;
• nie reprezentować mozaiki;
• mieć niepewne SV;
• pochodzić z dobrze scharakteryzowanej linii.

Precision i recall powinny być raportowane wraz z benchmark region.

Artefakty rosną z głębokością

Wysokie coverage ujawnia systematyczne błędy, których nie widać przy małej głębokości. Li opisał artefakty wariantów w bardzo głęboko sekwencjonowanych próbkach, związane między innymi z mapowaniem, paralogami i błędami technicznymi.⁷

Więcej reads może zwiększyć statystyczną „pewność” błędnego modelu. Jakość danych nie rośnie liniowo z depth.

Diagnostyka według wzoru

Niski depth w tej samej pozycji we wszystkich próbkach

Podejrzenia:

• projekt sond;
• GC;
• repeat;
• mapowalność;
• adnotacja;
• stały dropout.

Niski depth tylko w jednej próbce

Podejrzenia:

• input;
• degradacja;
• biblioteka;
• CNV;
• wariant pod sondą;
• zamiana próbki.

Wysoki depth i niskie MAPQ

Podejrzenia:

• collapsed paralogs;
• multi-mapping;
• pseudogene;
• repeat.

Wysoki depth i wiele duplikatów

Podejrzenia:

• mała złożoność;
• nadmierny PCR;
• zbyt głębokie sekwencjonowanie;
• mały target.

VAF zaniżony dla heterozygot

Podejrzenia:

• reference bias;
• capture dropout;
• PCR dropout;
• CNV;
• mosaicism;
• kontaminacja.

Studium 1: heterozygota w trudnym eksonie

Panel raportuje 180×, ale wariant ma VAF 12%.

Krok po kroku:

1. Globalne 180× nie jest coverage locus.
2. Locus ma 95 reads po jakości.
3. Tylko 28 ma MAPQ ≥20.
4. Po UMI zostaje 14 rodzin.
5. Większość alternatywnych reads ma secondary alignment do pseudogenu.
6. Region zawiera wariant pod sondą capture.
7. CNV model nie wskazuje delecji.

Możliwe hipotezy:

• reference/capture bias;
• pseudogene contamination;
• mozaika;
• prawdziwy allele imbalance.

Potwierdzenie powinno użyć long-range PCR z unikalnych kotwic albo long reads, a nie tego samego ampliconu.

Studium 2: mozaika 1% w FFPE

Surowy depth wynosi 5000×, a 42 reads zawierają alternatywę: VAF 0,84%.

Po analizie:

• 5000 reads tworzy 320 rodzin UMI;
• alternatywa jest w 5 rodzinach;
• wszystkie mają tę samą orientację;
• zmiana C>T leży na końcach;
• blank FFPE ma podobny sygnał;
• nie ma duplex families.

Wysokie raw depth nie pokonuje artefaktu deaminacji. Wynik jest poniżej zwalidowanego LoD lub no-call, a nie potwierdzoną mozaiką.

Studium 3: brak wariantu w ekspansji

WGS ma mean 35×, ale żaden read nie obejmuje całego powtórzenia.

Dane:

• unikalne kotwice leżą 2 kb od siebie;
• insert short-read wynosi 450 bp;
• reads wewnątrz repeat mają MAPQ 0;
• lokalne depth jest wysokie;
• caller nie raportuje wariantu.

Brak w VCF nie jest wynikiem ujemnym. Assay nie ma informacyjnej długości. Potrzebny jest celowany long-read, repeat-primed PCR lub inna zwalidowana metoda.

Studium 4: apparent deletion w panelu

Jeden ekson ma połowę depth względem kohorty.

Możliwe przyczyny:

• heterozygotyczna delecja;
• słaba biblioteka;
• wariant pod sondą;
• GC;
• batch;
• różnica transkryptu;
• błędna normalizacja.

Dowody CNV:

• spójny spadek wielu sond;
• allele balance sąsiednich SNP;
• split reads;
• powtórzenie w drugiej bibliotece;
• qPCR/MLPA.

Pojedynczy słaby bait nie wystarcza.

Studium 5: allele-specific expression

DNA pokazuje heterozygotę 50:50. RNA ma 90:10.

Przed wnioskiem o regulacji:

1. sprawdza się genotype;
2. usuwa mapping bias;
3. analizuje informative reads;
4. kontroluje copy number;
5. rozważa nonsense-mediated decay;
6. sprawdza imprinting/X;
7. modeluje dawców;
8. fazuje z wariantem regulatorowym.

Jedna osoba i jeden SNP tworzą hipotezę, nie uniwersalny mechanizm.

Decyzja: czy zwiększyć głębokość

Tak, jeśli

• biblioteka ma wysoką złożoność;
• locus jest mapowalny;
• molecular depth rośnie;
• saturation jest niska;
• brak wynika z losowania;
• assay ma odpowiednią klasę wariantu.

Nie, jeśli

• wszystkie nowe reads są duplikatami;
• region ma MAPQ 0;
• allele dropout powstał przed biblioteką;
• wariant jest dłuższy niż informacyjna cząsteczka;
• tło przekracza oczekiwany VAF;
• target nie jest objęty;
• referencja nie reprezentuje sekwencji.

Alternatywą jest zmiana biblioteki, długości, enrichment, referencji lub metody.

Decyzja: czy region jest callable

1. Czy target istnieje w referencji i panelu?
2. Czy reads mają unikalne kotwice?
3. Czy jest wystarczający molecular depth?
4. Czy oba allele są widoczne?
5. Czy klasa i wielkość wariantu były walidowane?
6. Czy strand, position i orientation są ocenialne?
7. Czy CN i ploidalność są znane?
8. Czy tło błędu jest poniżej limitu?
9. Czy benchmark obejmuje ten kontekst?
10. Czy negatywny wynik ma odpowiednią czułość?

Jedna odpowiedź „nie” może zmienić wynik na limited albo not assessed.

Pokrycie w eksperymentach edycji genomu

Po CRISPR trzeba ocenić:

• on-target indels;
• duże delecje;
• insercję template/wektora;
• rearanżacje;
• mozaikę klonu;
• off-target;
• copy number;
• allele dropout.

Krótki amplicon może preferencyjnie amplifikować nieedytowany allel, jeśli duża delecja usuwa starter. Wynik „100% wild type” może być fałszywy.

Silniejszy projekt łączy:

• startery poza cięciem;
• long-range assay;
• ddPCR/CNV;
• long reads;
• karyotyp;
• niezależne klony.

Pokrycie w genetyce sądowej

W mieszaninie śladów allele mają różne udziały. PCR dropout, stutter i mała liczba kopii wpływają na profile.

NGS depth nie powinien być bezpośrednio utożsamiany z ilością DNA osoby. Potrzebne są:

• locus balance;
• heterozygote balance;
• stutter model;
• noise threshold;
• replicate;
• negative controls;
• probabilistic genotyping.

W mitochondrialnym DNA dochodzi heteroplasmy i NUMTs. Chain of custody obejmuje także pliki, sample sheet i wersję pipeline.

Jak raportować jakość pokrycia

Minimalny raport:

• target i wersja adnotacji;
• raw reads/bases;
• przypisanie do próbek;
• mapped i properly paired;
• duplicate rate;
• unique/molecular reads;
• mean i median depth;
• breadth przy progach;
• uniformity;
• GC bias;
• on-target;
• MAPQ;
• callable region per klasa;
• lista dropout;
• contamination;
• ograniczenia paralogów;
• wersje pipeline.

Licznik powinien mieć mianownik

Każda liczba jakości wymaga definicji:

• procent czego;
• przed czy po filtrze;
• wszystkie reads czy unikalne;
• target projektowany czy callable;
• bazy czy molekuły;
• próbka czy kohorta;
• autosomy czy cały genom.

„98% pokrycia” bez progu i mianownika jest nieinterpretowalne. Może znaczyć 98% targetu z ≥1 readem albo 98% callable targetu z ≥30 niezależnymi molekułami — to zupełnie inne zapewnienie.

Rozkład per próbka, nie tylko batch

Zbiorczy boxplot może ukryć pojedynczą nieudaną próbkę. Dla każdej należy zachować metryki, a kohortowo monitorować:

• medianę;
• zakres;
• outliers;
• trend po dacie;
• loty;
• instrument;
• operatora.

Proces stabilny ma przewidywalny rozkład. Przesunięcie całej serii może pojawić się, zanim próbki przekroczą formalny próg odrzucenia.

Raport per gen

Dla genu klinicznego lub kluczowego eksperymentalnie trzeba podać:

• pokrycie eksonów i splice;
• regiony homologiczne;
• znane pseudogeny;
• CNV capability;
• warianty poza zakresem;
• alternatywną metodę uzupełniającą.

Najczęstsze błędy

1. „30× oznacza 30 niezależnych molekuł.” PCR może zwielokrotnić jedną.
2. „Średnia 30× oznacza każdą pozycję 30×.” Rozkład jest nierówny.
3. „Q40 baza jest dobrze zmapowana.” Base Q i MAPQ to inne osie.
4. „Dłuższy read zawsze jest unikalny.” Potrzebuje kotwicy.
5. „Duplikat po współrzędnych jest PCR.” Może być natural collision.
6. „On-target 90% oznacza kompletny panel.” Nie mówi o dropout.
7. „50% metylacji oznacza jeden allel.” Może być mieszaniną.
8. „Więcej reads usuwa bias.” Systematyczny bias może się utrwalić.
9. „Brak w VCF to wynik ujemny.” Region może być no-call.
10. „Jedno callable mask pasuje do wszystkich wariantów.” Zależy od klasy.
11. „Potwierdzenie PCR innym sekwenatorem jest ortogonalne.” Starter może powtórzyć dropout.
12. „Benchmark obejmuje cały genom.” Truth set ma granice.

Co warto umieć po przeczytaniu

Czytelnik powinien potrafić:

• rozdzielić depth, breadth, physical, spanning i molecular coverage;
• policzyć prosty model Poissona i wyjaśnić jego założenia;
• interpretować complexity i duplikaty;
• odróżnić base quality od MAPQ;
• wyjaśnić mapowalność;
• rozpoznać paralog, pseudogen i reference bias;
• ocenić GC, capture i PCR dropout;
• połączyć coverage z CNV i mozaiką;
• zdefiniować callable region per klasa;
• odróżnić negative od not assessed;
• zaplanować benchmark trudnych loci;
• przygotować raport per cel i gen.

Literatura i źródła