Składanie i adnotacja genomu

Streszczenie

Złożenie genomu rekonstruuje długą sekwencję z zachodzących na siebie odczytów. Adnotacja przypisuje tej sekwencji geny, transkrypty, powtórzenia i elementy regulatorowe. Oba etapy są wnioskowaniem z niepełnych danych, dlatego kompletność i dokładność trzeba mierzyć osobno.

Nie istnieje jedna liczba oznaczająca „dobry genom”. Złożenie może być bardzo ciągłe, a jednocześnie zawierać błędne połączenie chromosomów; może mieć niemal wszystkie geny, a gubić centromery; może trafnie reprezentować większość sekwencji, lecz mieszać dwa haplotypy. Adnotacja może zawierać tysiące modeli genów, ale nadawać funkcje na podstawie słabego podobieństwa.

Ten tekst prowadzi przez pełny projekt: wymagania próbki, wybór odczytów, graf złożenia, usuwanie zanieczyszczeń i duplikatów haplotypowych, fazowanie, scaffolding, polishing, ocenę bez referencji i względem referencji, adnotację powtórzeń oraz genów, wersjonowanie i publikację. Celem nie jest lista programów, lecz umiejętność rozpoznania, jaki dowód wspiera każdą decyzję.

Najpierw definiuje się produkt złożenia

„Złożyć genom” może oznaczać bardzo różne produkty:

• jeden konsensus wystarczający do identyfikacji bakterii;
• zamknięty chromosom i plazmidy izolatu;
• draft genomu nieznanego gatunku;
• chromosomowe złożenie do badań populacyjnych;
• dwie rozdzielone sekwencje haplotypów osobnika diploidalnego;
• złożenie telomere-to-telomere;
• referencję pangenomową z wielu haplotypów;
• metagenome-assembled genomes z mieszaniny organizmów.

Każdy produkt ma inne kryteria ukończenia. Dla bakterii pojedynczy kolisty kontig może być realistycznym celem, ale trzeba jeszcze wykazać, że nie jest błędnie zamknięty. Dla dużego genomu roślinnego z poliploidią „liczba chromosomów” nie mówi, czy homeologi zostały rozdzielone. Dla człowieka primary i alternate contigs nie są tym samym co dwa kompletne haplotypy.

Przed wyborem assemblera zapisuje się:

1. oczekiwany rozmiar i ploidalność;
2. poziom heterozygotyczności;
3. udział powtórzeń;
4. organella i elementy pozachromosomowe;
5. oczekiwaną ciągłość;
6. wymaganie fazy;
7. dopuszczalne luki;
8. zastosowanie końcowe;
9. kryteria poprawności;
10. format publikowanego produktu.

Zastosowanie ustala tolerancję błędu

Genom do szybkiej identyfikacji gatunku może tolerować fragmentację, której nie można zaakceptować przy analizie ewolucji rodzin genowych. Referencja do mapowania RNA-seq potrzebuje dobrze opisanych loci i izoform. Genom do projektowania CRISPR wymaga poprawnych paralogów i sekwencji off-target. Złożenie do analizy wariantów strukturalnych musi zachowywać rzeczywistą architekturę, a nie tylko wysokie N50.

Nie istnieje więc uniwersalny próg „assembly quality”. Kryteria powinny wynikać z decyzji, którą genom ma wspierać.

Próbka ogranicza graf przed uruchomieniem programu

Assembler nie odróżni mechanicznie dwóch osobników tego samego gatunku od dwóch haplotypów jednego osobnika, jeśli metadane tego nie ujawnią. Materiał powinien pochodzić z jasno zdefiniowanej jednostki biologicznej.

Trzeba ustalić:

• czy próbka jest klonalna, inbredowa czy outbredowa;
• czy zawiera wiele tkanek lub osobników;
• czy gatunek jest diploidalny, poliploidalny albo aneuploidalny;
• czy występuje mozaicyzm;
• czy materiał zawiera symbionty, pasożyty lub mikrobiom;
• czy organella mają liczne haplotypy;
• czy DNA pochodzi z jednej izolacji;
• czy dostępne są próbki rodziców.

Pool wielu osobników zwiększa liczbę alleli w grafie. Może poprawić reprezentację zmienności populacji, ale zwykle utrudnia złożenie pojedynczego genomu.

HMW, HiFi i Hi-C wymagają zgodności tożsamości

Projekt hybrydowy często używa:

• DNA HMW do długich reads;
• osobnej ekstrakcji do krótkich reads;
• świeżej tkanki do Hi-C;
• RNA z wielu tkanek do adnotacji;
• materiału rodziców do trio binning.

Wszystkie te źródła powinny dotyczyć tego samego osobnika albo jasno opisanej relacji. Zamiana próbki między biblioteką HiFi i Hi-C może tworzyć pozornie trudny problem fazowania. RNA innego genotypu nadal pomoże przewidzieć geny, ale warianty i brakujące eksony mogą mylić mapowanie.

Fingerprint genetyczny na wspólnym zestawie markerów jest tańszy niż późniejsze wyjaśnianie sprzecznego grafu.

Ilość danych liczy się w informacyjnych cząsteczkach

Średnie pokrycie można oszacować:

C = liczba odczytanych baz / haploidalny rozmiar genomu

Jeżeli powstaje 90 Gb danych dla genomu 3 Gb, nominalne pokrycie wynosi 30×. Ten rachunek zakłada jednak:

• właściwy rozmiar genomu;
• brak dużej kontaminacji;
• równomierne losowanie;
• brak duplikatów;
• poprawne przypisanie reads;
• reprezentację obu haplotypów.

Po filtrach może pozostać znacznie mniej. Długie reads o tej samej liczbie baz mogą różnić się zdolnością do przekraczania powtórzeń. Projekt powinien raportować rozkład długości i pokrycie fizyczne, nie tylko sumę gigabaz.

Rozmiar i złożoność genomu można oszacować z k-merów

Przed assembly krótkie, dokładne reads rozbija się na wszystkie słowa długości k. Histogram liczebności k-merów pokazuje, ile słów występuje raz, kilka razy lub bardzo często.

W prostym genomie diploidalnym mogą pojawić się:

• lewy ogon k-merów błędnych o małej liczebności;
• pik heterozygotyczny przy około połowie głównego pokrycia;
• pik homozygotyczny;
• piki powtórzeń o wielokrotnościach pokrycia;
• bardzo częste k-mery organelli i satelitów.

Z modelu można oszacować:

• rozmiar haploidalny;
• heterozygotyczność;
• udział sekwencji powtarzalnej;
• duplikację;
• poziom błędu;
• możliwą kontaminację.

To estymacja zależna od k, jakości i modelu. Genomy poliploidalne, silnie powtarzalne lub mieszaniny mogą mieć wiele nakładających się pików.

K-mer jest później niezależną jednostką kontroli

K-mery z dokładnych reads można porównać ze złożeniem bez mapowania do referencji. Każdy k-mer obserwowany w danych powinien być:

• obecny w odpowiedniej liczbie kopii;
• przypisany do właściwego haplotypu;
• niepowielony artefaktycznie;
• wolny od k-merów powstałych wyłącznie przez błąd assembly.

Taki bilans wykrywa utratę allelu, haplotig duplication i błędy zasad, które mogą umknąć przy samym alignmencie.

Dane wejściowe muszą przejść własne QC

Przed assembly ocenia się osobno każdy typ danych.

Długie odczyty

• yield po filtrach;
• N50 i pełny rozkład długości;
• jakość per read;
• obecność adapterów;
• chimery;
• kontaminację;
• pokrycie obu haplotypów;
• udział reads organellowych;
• sygnał degradacji próbki.

Krótkie odczyty

• Q per cycle;
• adapter read-through;
• duplikaty;
• bias GC;
• insert size;
• kontaminację;
• zgodność par;
• histogram k-merów.

Hi-C

• odsetek ważnych par;
• duplikaty;
• cis/trans ratio;
• rozkład odległości kontaktów cis;
• zgodność miejsca restrykcyjnego lub chemii;
• biblioteki chimeryczne;
• kontakt między haplotypami i chromosomami.

RNA do adnotacji

• tkanki, warunki i etapy rozwoju;
• integralność RNA;
• strandedness;
• długość i typ biblioteki;
• pokrycie końców transkryptów;
• rRNA i kontaminację;
• zgodność osobnika.

Assembler uruchomiony na niezweryfikowanych danych może zakończyć się poprawnie technicznie i wiernie zrekonstruować zanieczyszczenie.

Od odczytu do grafu

Assembler nie układa plików jak puzzli według obrazu pudełka. Buduje graf podobieństw. W grafie nakładania wierzchołkami są odczyty, a krawędzie oznaczają zgodne końce. W grafie de Bruijna odczyty rozbija się na k-mery, a ścieżka przechodzi przez ich kolejne nakładania.

Sekwencja unikalna daje prostą ścieżkę. Powtórzenie dłuższe niż informacja w odczycie tworzy rozgałęzienie. Heterozygotyczność może wyglądać jak dwie alternatywne ścieżki, a błąd jak rzadka odnoga.

Graf de Bruijna

W grafie de Bruijna k-mery są redukowane do relacji prefiks–sufiks. Zwiększenie k może rozdzielić krótkie powtórzenia, lecz wymaga większego pokrycia i sprawia, że błąd rozbija więcej k-merów. Mniejsze k lepiej łączy dane o małym pokryciu, ale skleja więcej podobnych regionów.

Krótkie reads nie są przechowywane w grafie jako pełne historie cząsteczek. Para końców i insert size dodają informację o połączeniu, ale długie powtórzenie nadal może pozostać nierozstrzygnięte.

Overlap–layout–consensus

Dla długich reads można szukać nakładania całych cząsteczek. Etapy obejmują:

1. wykrycie kandydatów overlaps;
2. odrzucenie słabych i powtarzalnych dopasowań;
3. zbudowanie layoutu;
4. wybranie ścieżek;
5. utworzenie konsensusu.

Koszt porównania każdej pary byłby ogromny, więc programy używają minimizerów i indeksów. Powtórzenia, chimery i nierówne pokrycie nadal tworzą fałszywe krawędzie.

String graph i phased assembly graph

Graf może usuwać reads całkowicie zawarte w innych i zachowywać istotne relacje nakładania. Warianty heterozygotyczne tworzą bubbles. Program fazujący próbuje rozpoznać, które warianty współwystępują na tych samych cząsteczkach, i prowadzić osobne ścieżki haplotypów.

Bubble nie zawsze jest allelem. Może być:

• błędem sekwencjonowania;
• powtórzeniem;
• duplikacją segmentową;
• kontaminacją;
• wariantem somatycznym;
• prawdziwą różnicą haplotypów.

Decyzja wymaga pokrycia, jakości, długości kotwic i spójności z innymi markerami.

Wybór realnego assemblera zależy od typu odczytów. Dla dokładnych długich odczytów PacBio HiFi standardem jest hifiasm (złożenie z rozdzieleniem haplotypów).⁴ Dla odczytów nanopore używa się Flye lub Canu (podejście nakładania, OLC). Krótkie, dokładne odczyty i małe genomy (bakterie) składa SPAdes (graf de Bruijna). Do złożeń kompletnych, „telomer-do-telomeru", łączy się HiFi z ultradługimi odczytami nanopore w narzędziach takich jak Verkko i hifiasm (UL). Nazwa programu jest jednak mniej ważna niż dopasowanie do danych i zapisane parametry.

Kontig, scaffold i chromosom

Kontig jest ciągłą sekwencją bez przerw. Scaffold porządkuje i orientuje kontigi, czasem zostawiając nieznane fragmenty. Dane Hi-C o kontaktach chromatyny pomagają przypisać scaffoldy do chromosomów — służą do tego narzędzia SALSA2, YaHS oraz Juicer/3D-DNA, a mapy optyczne Bionano dają niezależną informację o strukturze. Błędne połączenie może jednak wyglądać wiarygodnie, dlatego scaffolding weryfikuje się mapami kontaktów.

Złożenie telomere-to-telomere prowadzi od jednego końca chromosomu do drugiego bez luk. Pierwsza kompletna ludzka sekwencja T2T dotyczyła szczególnej, niemal homozygotycznej linii CHM13, co zmniejszało problem dwóch haplotypów.¹ Genom zwykłej diploidalnej osoby wymaga rozdzielenia obu wersji.

Luka jest informacją o nierozstrzygnięciu

Scaffold może zawierać znaki N, które reprezentują sekwencję o nieznanej treści lub długości. Ich obecność nie jest wyłącznie estetycznym brakiem. Mówi, gdzie dane pozwoliły uporządkować kontigi, lecz nie odtworzyć łącznika.

Raport powinien podać:

• liczbę luk;
• łączną liczbę N;
• szacowaną lub nieznaną długość;
• typ dowodu łączącego;
• pozycję względem centromerów, telomerów i genów;
• czy luka występuje w obu haplotypach.

Zastąpienie nieznanej sekwencji dowolnym fragmentem tylko po to, aby „zamknąć” chromosom, jest gorsze niż jawna luka.

Scaffolding dodaje relacje dalekiego zasięgu

Kontig jest wynikiem lokalnej ciągłości sekwencji. Scaffolding używa informacji, że dwa fragmenty są blisko na tej samej cząsteczce lub w tym samym chromosomie.

Źródła mogą obejmować:

• ultradługie reads;
• linked reads i barcodes cząsteczkowe;
• mapy optyczne;
• Hi-C;
• mapy genetyczne;
• synteny z bliskim gatunkiem;
• ręczną analizę centromerów i telomerów.

Każde ma inny zasięg i rozdzielczość. Hi-C daje relacje chromosomowe, ale nie sekwencję luki. Mapa optyczna mierzy wzór znaczników na długich molekułach, lecz nie rozstrzyga pojedynczych zasad. Bliski genom referencyjny może pomóc w orientacji, ale może też wymusić strukturę, której badany osobnik nie ma.

Hi-C nie jest prostą linijką

Komórki są utrwalane, chromatyna cięta i ligowana tak, by zachować informację o kontaktach przestrzennych. Pary reads wskazują fragmenty, które były blisko w jądrze. Kontakty wewnątrz chromosomu są zwykle częstsze niż między chromosomami, a częstość maleje z odległością genomową.

Na mapie prawidłowego chromosomu oczekuje się silnej przekątnej. Potencjalna błędna orientacja, translokacja lub misjoin może tworzyć przerwę albo nietypowy blok kontaktów. Jednak biologia jądra także tworzy domeny, przedziały i kontakty dalekiego zasięgu.

Kuracja Hi-C wymaga więc rozróżnienia:

• błędu assembly;
• rzeczywistego wariantu strukturalnego;
• haplotypu;
• biologicznej organizacji chromatyny;
• artefaktu biblioteki.

Referencyjne scaffolding może ukryć nowość

Ułożenie kontigów według genomu bliskiego gatunku jest kuszące, bo szybko tworzy pseudomolekuły. Ale inwersja lub translokacja może zostać „naprawiona” do stanu referencji. Taki produkt nadaje się czasem jako robocza mapa, lecz nie jest niezależnym dowodem struktury badanego genomu.

Jeśli użyto reference-guided scaffolding, raport powinien to zaznaczyć i oddzielić połączenia wspierane własnymi danymi.

Haplotypy mogą być rozdzielone na kilka sposobów

Diploidalny genom zawiera pary homologów. W regionie heterozygotycznym assembler może:

• zwinąć allele do jednego konsensusu;
• wybrać jeden allel;
• utworzyć primary contig i alternate haplotig;
• rozdzielić dwie pełne ścieżki;
• przełączać się między haplotypami.

Ostatni przypadek tworzy haplotyp mozaikowy, którego nie było w komórce.

Trio binning

Dokładne krótkie reads rodziców dostarczają k-merów specyficznych dla linii matczynej i ojcowskiej. Długie reads potomka są przypisywane do koszyków, a następnie składane oddzielnie.

Zalety:

• bezpośrednie pochodzenie rodzicielskie;
• rozdzielenie długich bloków;
• mniejsze mieszanie alleli.

Ograniczenia:

• potrzeba danych rodziców;
• regiony bez informacyjnych markerów;
• de novo variants;
• błędy lub kontaminacja próbek;
• pokrewieństwo i niska różnorodność;
• brak prostej analogii w poliploidach.

Fazowanie read-backed

Długi read obejmujący wiele heterozygotycznych markerów wiąże je w cis. Nakładające się reads propagują fazę wzdłuż chromosomu. Przerwa o małej heterozygotyczności lub długie identyczne powtórzenie może zakończyć blok.

Dane Hi-C mogą połączyć fazę na większych odległościach, ale pojedyncza para kontaktowa ma inną pewność niż jedna ciągła cząsteczka.

Switch error i Hamming error

Switch error występuje, gdy od pewnego miejsca allele zostają zamienione między haplotypami. Hamming error mierzy udział markerów przypisanych do niewłaściwej ścieżki. Długi blok fazy nie gwarantuje małej liczby przełączeń.

Fazę ocenia się za pomocą:

• danych rodziców;
• markerów populacyjnych;
• niezależnych długich cząsteczek;
• Hi-C;
• porównania k-merów haplotypowych;
• znanych regionów chromosomu X lub innych kontroli.

Haplotig duplication i collapse są przeciwnymi błędami

Jeśli assembler zachowa dwa allele jako dwa niezależne locus w haploidalnej reprezentacji, rozmiar rośnie, pokrycie każdego fragmentu spada, a geny „jednokopijne” wyglądają na zduplikowane. To haplotig duplication.

Jeśli podobne paralogi albo homeologi zostaną połączone w jeden fragment, powstaje collapse. Pokrycie takiego regionu może być nienaturalnie wysokie, a warianty wyglądają jak heterozygotyczność.

Rozpoznanie wykorzystuje razem:

• depth;
• podobieństwo contigów;
• k-mery;
• BUSCO duplication;
• synteny;
• strukturę grafu;
• fazę;
• oczekiwaną ploidalność.

Automatyczne purge może usunąć prawdziwą duplikację. Przed skasowaniem sekwencji trzeba sprawdzić, czy jest allelem, paralogiem, organellum czy kontaminantem.

Poliploidia wymaga jawnego modelu

W autotetraploidzie mogą istnieć cztery podobne homologiczne kopie. W allopoliploidzie subgenomy pochodzą od różnych przodków i mogą zawierać homeologi. Model diploidalny może:

• scalić część kopii;
• sztucznie podzielić inne;
• błędnie usuwać „duplikaty”;
• tworzyć chimery między subgenomami.

Ploidalność nie jest parametrem kosmetycznym. Wyznacza liczbę oczekiwanych ścieżek i interpretację pokrycia.

Kontaminacja może wyglądać jak interesująca sekwencja

Kontigi bakterii, grzyba, pasożyta, wektora, człowieka lub organizmu pokarmowego mogą znaleźć się w assembly. Ich identyfikacja nie powinna opierać się na jednym najlepszym trafieniu BLAST.

Łączy się:

• pokrycie;
• GC;
• taksonomiczne podobieństwo;
• rozkład k-merów;
• obecność intronów;
• skład genów;
• połączenia w długich reads;
• Hi-C;
• obecność w niezależnych ekstrakcjach.

Sekwencja o bakteryjnym podobieństwie może być prawdziwym transferem horyzontalnym albo kontaminacją. Dowodem integracji są reads obejmujące granicę z unikalną sekwencją gospodarza, powtarzalność między próbkami i właściwe pokrycie.

Organella i symbionty wymagają osobnej decyzji

Mitochondrium i chloroplast mają zwykle większą liczbę kopii, więc tworzą wysokie pokrycie i mogą być składane oddzielnie. W genomie jądrowym istnieją fragmenty pochodzenia mitochondrialnego lub plastydowego. Automatyczne usunięcie wszystkiego podobnego do organellum skasuje NUMTs lub NUPTs.

Symbiont może być biologicznie ważnym członkiem próbki, nie „brudem”. Produkt powinien jawnie klasyfikować sekwencje:

• nuclear assembly;
• organellar assembly;
• symbiont;
• unplaced;
• contaminant;
• uncertain.

Polerowanie i kontrola jakości

Po wstępnym złożeniu odczyty mapuje się z powrotem (minimap2/BWA), korygując błędy — to polerowanie (polishing). Realne narzędzia dobiera się do danych: Pilon (krótkie, dokładne odczyty), Racon i Medaka (nanopore), NextPolish. Krótkie dokładne odczyty dobrze poprawiają pojedyncze zasady, a długie utrzymują strukturę. Zbyt agresywne polerowanie może jednak zamienić prawdziwy wariant jednego haplotypu w sekwencję drugiego.

N50 to popularna miara: wartość, dla której połowa złożenia leży w fragmentach nie krótszych niż ona. Metryki ciągłości (N50 i inne) liczy QUAST. Brzmi solidnie, ale długi, błędnie połączony scaffold daje świetny N50 — długość nie jest tym samym co poprawność. Dlatego kompletność ocenia się osobno, najczęściej narzędziem BUSCO, które sprawdza, ile z zestawu genów konserwowanych i zwykle jednokopijnych złożenie zawiera jako pełne, zduplikowane, fragmentaryczne lub brakujące.² Niezależną, wolną od referencji ocenę daje Merqury — z rozkładu k-merów odczytów szacuje jakość zasad (QV), kompletność i poprawność fazowania. Zduplikowane geny mogą zresztą sygnalizować, że dwa haplotypy zostały omyłkowo złożone jako osobne sekwencje. Potrzebne są więc zgodność z odczytami, liczba brakujących genów konserwowanych, k-mery, kompletność chromosomów i wykrywanie duplikacji. Jakość próbki często ogranicza wynik bardziej niż algorytm.

Polishing nie może nadpisywać biologii

Każda runda powinna mieć określone dane i cel:

• korekta losowych indels długiego odczytu;
• poprawa homopolimerów;
• utworzenie konsensusu HiFi;
• korekta małych błędów krótkimi reads;
• osobne traktowanie organelli;
• zachowanie fazy.

Po każdej rundzie porównuje się QV, k-mery, liczbę zmian, BUSCO i fazę. Jeśli wynik przestaje się poprawiać albo rośnie liczba switch errors, kolejna iteracja szkodzi.

Mapowanie reads do własnego assembly nie jest całkowicie niezależną kontrolą. Algorytm preferuje istniejącą sekwencję, a powtórzenia mogą mapować się niejednoznacznie. K-mery i dane ortogonalne zmniejszają tę zależność.

Jakość złożenia ma kilka osi

Ciągłość

Raportuje się:

• liczbę contigów i scaffoldów;
• contig N50 i scaffold N50;
• L50;
• największy fragment;
• liczbę i długość luk;
• udział sekwencji przypisanej do chromosomów;
• telomery i centromery.

NG50 używa oczekiwanego rozmiaru genomu zamiast sumy assembly, dzięki czemu zbyt małe złożenie nie poprawia sobie wyniku samym pominięciem trudnej sekwencji.

Poprawność zasad

QV jest logarytmiczną estymatą błędu. QV 40 odpowiada około jednemu błędowi na 10 tys. zasad, QV 50 — na 100 tys., a QV 60 — na milion. Jedna wartość nie pokazuje, czy błędy skupiają się w homopolimerach, tandem repeats albo konkretnym haplotypie.

Warto raportować:

• globalny QV;
• zgodność k-merów;
• indels i substitutions osobno;
• trudne konteksty;
• jakość organelli;
• regiony nieocenialne.

Kompletność

BUSCO pyta o zestaw genów oczekiwanych w danej linii ewolucyjnej. Wynik zależy od wybranego lineage dataset. Wysoka kompletność genowa nie dowodzi kompletności satelitów, centromerów ani intergenicznych regulatorów.

Kompletność k-merowa odpowiada na inne pytanie: jaki udział słów z reads znalazł się w assembly. Można mieć pełne BUSCO i nadal utracić jeden haplotyp w wielu regionach.

Poprawność struktury

Misassembly obejmuje:

• błędne połączenie odległych loci;
• inwersję;
• translokację;
• collapse;
• sztuczną duplikację;
• brak kopii tandemowej;
• switch haplotypu.

Kontrola używa długich reads obejmujących breakpoint, map optycznych, Hi-C, synteny, depth i grafu. Synteny sama nie rozstrzyga, bo prawdziwy genom może różnić się od krewniaka.

Faza

Raportuje się:

• phase block N50;
• switch error;
• Hamming error;
• udział markerów fazowanych;
• kompletność każdego haplotypu;
• k-mery rodzicielskie w złej ścieżce;
• regiony collapsed.

Jedna primary assembly plus zestaw alternate contigs wymaga innego opisu niż dwa symetryczne haplotypy.

Zanieczyszczenia i identyfikowalność

Należy wykazać, że:

• główne sekwencje należą do oczekiwanego taksonu;
• obce kontigi sklasyfikowano;
• organella są oznaczone;
• wektory i adaptery sprawdzono;
• wszystkie typy danych mają zgodny fingerprint;
• wersja próbki, biblioteki i pipeline jest znana.

Dwa haplotypy i fazowanie

Człowiek jest diploidalny: ma dwie kopie każdego chromosomu, nieco różne. Dobre współczesne złożenie nie „uśrednia” ich w jedną sekwencję, lecz fazuje — rozdziela na dwa haplotypy. Pomaga w tym sekwencjonowanie rodziców (tzw. trio binning), które pozwala przypisać każdy fragment do linii matczynej lub ojcowskiej. Bez fazowania warianty leżące blisko siebie mogą zostać błędnie przypisane, co ma znaczenie np. w diagnostyce.

Adnotacja powtórzeń i genów

Najpierw identyfikuje się transpozony, satelity i inne powtórzenia, aby nie pomylić ich z wieloma genami — de novo wykrywa je RepeatModeler, a maskuje RepeatMasker (dla roślin popularne jest EDTA). Dopiero potem przewiduje się geny: z sekwencji kodującej, sygnałów splicingu, danych RNA-seq, pełnych transkryptów i homologii do znanych białek. Realne predyktory to AUGUSTUS i BRAKER (łączący dowody z RNA i białek), MAKER, a coraz częściej oparty na uczeniu głębokim Helixer; dla genomów bakteryjnych standardem jest Prokka lub Bakta. Funkcję nadanym genom przypisuje się przez InterProScan (domeny) i eggNOG-mapper (ortologia → nazwy, terminy GO, szlaki).

Model genu obejmuje eksony, introny, alternatywne transkrypty i regiony niekodujące. Automatyczny predyktor może połączyć dwa geny albo pominąć specyficzny tkankowo ekson. Adnotacja jest wersjonowanym zbiorem hipotez, nie napisem ukrytym w DNA — dlatego podaje się jej wersję i dowody, na których się opiera.

Powtórzenia tworzą bibliotekę rodzin

Analiza powinna odróżnić:

• LINE i SINE;
• LTR retrotransposons;
• transpozony DNA;
• tandem repeats;
• satelity;
• proste powtórzenia;
• segmental duplications;
• pseudogeny i fragmenty genów.

Metoda de novo wykrywa rodziny częste w badanym genomie. Biblioteka kuratorowana dostarcza znanych klasyfikacji. Połączenie obu zmniejsza utratę sekwencji specyficznej dla gatunku.

Soft masking zmienia zapis zasad na małe litery, zachowując sekwencję. Hard masking zastępuje ją znakami N. Wybór wpływa na predyktory i powinien być raportowany.

Evidence-driven i ab initio odpowiadają na inne braki

Predykcja ab initio szuka wzoru start/stop, splicingu, długości eksonów i składu kodonów. Może znaleźć gen bez danych ekspresji, ale też stworzyć fałszywy model.

Dowody z RNA pokazują splice junctions i transkrypty obecne w zbadanych tkankach. Brak RNA nie dowodzi braku genu: ekspresja może być specyficzna dla etapu, stresu, płci lub rzadkiej komórki.

Homologia białkowa pomaga zachować geny konserwowane. Może jednak narzucić strukturę paralogu albo gatunku o innej organizacji eksonów.

Najlepszy model integruje niezależne dowody i zapisuje ich pochodzenie.

Krótkie i pełnej długości RNA uzupełniają się

Krótkie RNA-seq daje dużą głębokość i precyzyjne junctions, ale składa pełną izoformę pośrednio. Long-read cDNA lub direct RNA może połączyć odległe eksony na jednej cząsteczce, lecz ma inne błędy, mniejszy yield i bias końców.

Adnotacja transkryptów powinna kontrolować:

• kompletność 5′ i 3′;
• artefakty reverse transcription;
• internal priming;
• template switching;
• degradację;
• chimerę;
• liczbę niezależnych cząsteczek;
• zgodność splice sites.

Jeden locus może mieć wiele produktów

Trzeba rozdzielić:

• locus;
• gene model;
• transcript;
• CDS;
• UTR;
• białko;
• pseudogen.

Najdłuższy transkrypt nie zawsze jest biologicznie główny. Canonical transcript jest regułą wybraną dla konkretnego zastosowania, nie uniwersalnym faktem.

Adnotacja funkcji jest wnioskowaniem

Funkcję można przypisać na podstawie:

• eksperymentalnie opisanego ortologu;
• filogenezy rodziny;
• domeny białkowej;
• zachowanych reszt aktywnego miejsca;
• lokalizacji;
• ekspresji;
• sąsiedztwa genów;
• przewidywanej struktury;
• podobieństwa sekwencji.

„Najlepszy hit” nie musi być ortologiem. Po duplikacji paralogi mogą zmienić funkcję. Nazwa powinna odzwierciedlać siłę dowodu; słaby model lepiej pozostawić jako „putative” niż nadać mu precyzyjną rolę.

Kontrola adnotacji

Ocenia się:

• zgodność modeli z RNA i białkami;
• kompletność BUSCO proteomu;
• długości genów, intronów i białek;
• liczbę genów na scaffold;
• nadmiar bardzo krótkich modeli;
• geny w powtórzeniach;
• fuzje i rozcięcia modeli;
• porównanie z bliskimi gatunkami;
• udział białek z domeną lub ortologiem;
• ręczną kurację kluczowych loci.

Nienaturalnie duża liczba genów często oznacza niepełne maskowanie transpozonów albo haplotig duplication, nie biologiczną innowację.

Genom referencyjny i pangenom

Jedna liniowa referencja upraszcza mapowanie, lecz reprezentuje ograniczony zestaw alleli — a odczyty pasujące do niej „lepiej” tworzą uprzedzenie wobec sekwencji rzadszych. Pangenom używa grafu wielu haplotypów, zmniejszając to uprzedzenie; buduje się go i mapuje na niego narzędziami minigraph / minigraph-Cactus, pggb oraz vg (variation graph). W 2023 roku konsorcjum Human Pangenome Reference Consortium opublikowało pierwszy szkic ludzkiego pangenomu złożony z 47 diploidalnych, fazowanych złożeń różnorodnych osób, dokładając ponad sto milionów par zasad sekwencji nieobecnej w dotychczasowej pojedynczej referencji.³ Takie podejście lepiej oddaje zmienność populacji, ale zwiększa złożoność narzędzi i raportowania.

Graf pangenomu nie jest jednym „supergenomem”

W grafie węzły reprezentują odcinki sekwencji, krawędzie dopuszczalne przejścia, a haplotypy — ścieżki. Projekt wymaga decyzji, które złożenia włączyć, jak kontrolować ich jakość, kiedy fragmenty są allelami, a kiedy paralogami, jak reprezentować SV oraz jak utrzymać współrzędne. Niskiej jakości assembly może wnieść fałszywą gałąź. Więcej genomów nie automatycznie oznacza lepszą referencję.

Przenoszenie adnotacji między wersjami działa tylko tam, gdzie relacja sekwencji jest jednoznaczna. Nowa insercja, inwersja, naprawiony misassembly lub rozdzielenie paralogów może sprawić, że współrzędna nie ma odpowiednika, ma kilka odpowiedników albo zmienia orientację. Wariant powinien być identyfikowany także allelami i kontekstem, nie samym numerem pozycji.

Publikacja assembly wymaga więcej niż FASTA

Minimalny pakiet obejmuje:

• FASTA contigów i haplotypów;
• scaffolded FASTA oraz AGP;
• graf assembly, jeśli jest dostępny;
• GFF3/GTF, CDS i białka;
• bibliotekę repeats;
• raport QC;
• wersje narzędzi i parametry;
• provenance próbek oraz bibliotek;
• accession danych surowych;
• checksum, licencję i known issues.

Nazwy sekwencji powinny być stabilne. Zmiana sekwencji wymaga nowej wersji; cicha podmiana pliku uniemożliwia odtworzenie analiz. Dla unplaced i unlocalized scaffolds trzeba zachować jawny status. Nazwa haplotypu nie powinna sugerować pochodzenia rodzicielskiego, jeśli go nie wykazano.

Odtwarzalność wymaga zachowania decyzji

Sama lista programów nie odtwarza assembly. Potrzebne są:

• dokładne komendy i parametry;
• wersje programów, bibliotek i modeli;
• obrazy kontenerów albo opis środowiska;
• kolejność etapów;
• identyfikatory danych wejściowych;
• checksum;
• ziarna losowe, jeśli algorytm ich używa;
• reguły ręcznej kuracji;
• lista usuniętych i przywróconych contigów;
• kryteria zatrzymania polishing;
• raporty przed i po każdej większej zmianie.

Ręczna korekta mapy Hi-C może być uzasadniona, ale powinna pozostawić ślad: miejsce cięcia, połączenia, orientację, dowód i osobę zatwierdzającą. W przeciwnym razie finalny chromosom nie wynika z opublikowanego pipeline.

Assembly i adnotacja mają oddzielne cykle życia

Sekwencja może pozostać stabilna, gdy pojawiają się nowe dane RNA i zmieniają modele genów. Może też zostać poprawiona sekwencja przy zachowaniu większości adnotacji. Dlatego wersja assembly i annotation release powinny być oddzielne.

Przykładowy zapis metod powinien umożliwić odpowiedź:

• na jakiej wersji sekwencji leży gen;
• z którego release pochodzi model transkryptu;
• jaki FASTA był używany do mapowania;
• czy współrzędne zostały przeniesione;
• które modele nie przeszły lift-over;
• czy wyniki wariantów zostały przeliczone.

To szczególnie ważne, gdy poprawiony region rozdziela jeden dawny locus na dwa paralogi albo zamyka lukę w obrębie genu.

Cztery wzorcowe projekty

Bakteryjny izolat z plazmidami

Cel to zamknięty chromosom i pełne plazmidy. Sprawdza się czystość izolatu, używa długich reads do przejścia przez powtórzenia, dokładnych reads do korekty, a pokrycia do rozdzielenia chromosomu i plazmidów. Reads muszą obejmować miejsce domknięcia. Kolisty contig nie jest sam dowodem prawidłowego koła — assembler może błędnie połączyć powtórzone końce.

Diploidalny genom zwierzęcia

K-mery szacują rozmiar i heterozygotyczność, HiFi tworzy dokładne kontigi, ultradługie reads przekraczają powtórzenia, a trio lub Hi-C pomaga fazować. K-mery rodzicielskie mierzą przełączenia, depth i BUSCO wskazują collapse oraz haplotigs, a RNA wielu tkanek wspiera adnotację. Każdy haplotyp wymaga osobnego raportu kompletności.

Poliploidalna roślina

Projekt dodaje cytogenetyczne potwierdzenie ploidalności, markery subgenomowe, mapę genetyczną, analizę homeologów i szczególną ostrożność przy purge. Wysokie BUSCO duplication może być zgodne z biologią poliploidii. Automatyczna interpretacja „zduplikowane = błąd” byłaby fałszywa.

Złożenie metagenomiczne

W mieszaninie pokrycie zależy od obfitości, szczepy tworzą bubbles, a elementy mobilne łączą taksony. Po assembly wykonuje się binning według pokrycia, składu i współwystępowania, dereplikację, ocenę completeness i contamination oraz taksonomię. MAG może być kompozytem szczepów i nie jest równoważny genomowi izolatu.

Checklista pełnego projektu

Projekt i assembly

• produkt końcowy zdefiniowany;
• ploidalność i rozmiar oszacowane;
• zgodność wszystkich próbek sprawdzona;
• głębokość i długość odpowiadają repeats;
• assembler pasuje do danych;
• wersja, parametry i graf zachowane;
• źródło każdego scaffoldu znane;
• phase model jawny;
• purge kontrolowany;
• polishing oceniany po rundach.

Kontrola jakości

• contig i scaffold continuity;
• base QV;
• k-mer completeness;
• BUSCO;
• poprawność struktury;
• phase blocks i switch error;
• collapse i haplotig duplication;
• contamination;
• organella;
• zgodność ze spodziewanym rozmiarem;
• osobna ocena trudnych regionów.

Adnotacja i release

• biblioteka repeats de novo oraz kuratorowana;
• maskowanie opisane;
• RNA z odpowiednich tkanek;
• białka bliskich taksonów;
• dowody zapisane per model;
• funkcje z poziomem pewności;
• kluczowe loci ręcznie sprawdzone;
• wszystkie artefakty wersjonowane;
• dane surowe dostępne;
• known issues opublikowane.

Najczęstsze błędy

1. „Najwyższe N50 wygrywa.” Misjoin może je zwiększyć.
2. „BUSCO 99% oznacza kompletny genom.” Dotyczy wybranego zestawu genów.
3. „Więcej polishing zawsze pomaga.” Może mieszać haplotypy.
4. „Reference-guided scaffold jest prawdą.” Może skopiować strukturę referencji.
5. „Duplikat BUSCO trzeba usunąć.” Może być paralogiem lub poliploidią.
6. „Kontig bakterii to kontaminacja.” Może być symbiontem albo integracją.
7. „Dwa FASTA to dwa dobre haplotypy.” Trzeba zmierzyć switch error.
8. „Adnotacja jest odczytana z DNA.” Jest wersjonowanym modelem.
9. „Brak RNA oznacza brak genu.” Ekspresja zależy od kontekstu.
10. „Najlepszy hit nadaje funkcję.” Może być paralogiem.
11. „Pangenom automatycznie usuwa bias.” Jakość wejścia nadal decyduje.
12. „FASTA wystarczy.” Bez provenance produkt nie jest odtwarzalny.

Co warto umieć po przeczytaniu

Czytelnik powinien potrafić zdefiniować produkt assembly, interpretować histogram k-merów, rozróżnić graf de Bruijna, OLC i phased graph, odróżnić contig, scaffold i haplotyp, wyjaśnić Hi-C, rozpoznać collapse i haplotig duplication, zaplanować fazowanie, ocenić ciągłość, QV, kompletność, strukturę oraz fazę osobno, zbudować adnotację z wielu dowodów i opublikować wersjonowany produkt.

Następny tekst, Transkryptomika single-cell i spatial, pokazuje, jak adnotowany genom staje się mapą do pomiaru aktywności komórek.