Przeglądarka genomu nie jest mapą z jedną prawdą

Ensembl i UCSC Genome Browser pozwalają spojrzeć na ten sam fragment genomu, ale nie są dwiema kopiami jednej bazy. Oba systemy układają informacje względem sekwencji referencyjnej i pokazują je jako warstwy, nazywane zwykle ścieżkami (tracks). Jedna warstwa może przedstawiać modele genów, inna warianty, sygnał z ATAC-seq, powtórzenia, konserwację ewolucyjną albo wynik własnego eksperymentu.

To rozróżnienie jest podstawowe:

  • assembly dostarcza układu odniesienia i sekwencji;
  • współrzędne wskazują odcinek w tym układzie;
  • anotacja przypisuje odcinkowi obiekty i znaczenia;
  • warstwa dowodowa pokazuje pomiar albo wynik obliczeń;
  • interfejs decyduje, które elementy są widoczne i jak są narysowane.

To, że prostokąt genu nie pojawia się w danym widoku, nie oznacza jeszcze, że genu „nie ma”. Warstwa mogła być wyłączona, filtr mógł ukryć transkrypt albo dwie anotacje mogły inaczej wyznaczyć granice eksonów. Z kolei kolorowy pik w warstwie regulacyjnej nie dowodzi sam przez się, że region jest aktywnym enhancerem w każdej komórce.

Przeglądarka genomu jest więc stołem roboczym do zestawiania danych w kontekście, a nie wyrocznią. Najwięcej daje wtedy, gdy czytelnik potrafi powiedzieć, co jest osią współrzędnych, skąd pochodzi każda ścieżka i jaki rodzaj wniosku rzeczywiście wspiera.

Schemat czytania przeglądarki genomu. Współrzędna jest poprawnie określona dopiero razem z wersją assembly, nazwą sekwencji, konwencją przedziału i nicią. Na tej osi nakładane są niezależne warstwy: modele transkryptów, warianty, regulacja, konserwacja i własne dane. LiftOver mapuje odcinki między układami, lecz nie zastępuje kontroli sekwencji i ponownej anotacji.
Schemat czytania przeglądarki genomu. Współrzędna jest poprawnie określona dopiero razem z wersją assembly, nazwą sekwencji, konwencją przedziału i nicią. Na tej osi nakładane są niezależne warstwy: modele transkryptów, warianty, regulacja, konserwacja i własne dane. LiftOver mapuje odcinki między układami, lecz nie zastępuje kontroli sekwencji i ponownej anotacji.

Cztery składniki adresu genomowego

Zapis chr7:140453136-140453145 wygląda precyzyjnie, ale bez dodatkowych informacji jest niepełny. Poprawny adres powinien zawierać co najmniej:

  1. organizm i assembly, np. człowiek GRCh38;
  2. nazwę sekwencji, np. chromosom, contig lub accession z wersją;
  3. konwencję przedziału;
  4. nić, jeśli obiekt ma orientację biologiczną.

Można to zapisać w notatkach następująco:

Homo sapiens, GRCh38, chr7:140453136–140453145, współrzędne 1-based inclusive, nić minus.

Jeśli wynik ma wejść do programu, lepszy bywa identyfikator sekwencji z wersją, np. accession RefSeq albo GenBank assembly. Nazwa chr7 jest wygodna, lecz sposób nazwania sekwencji zależy od źródła pliku. Jeden plik może używać chr1, drugi 1, a trzeci accessionu NC_000001.11. Zanim uzna się brak dopasowania za wynik biologiczny, trzeba sprawdzić słownik nazw.

Assembly jest częścią wyniku, nie ustawieniem kosmetycznym

Assembly to uporządkowana kolekcja sekwencji referencyjnych. Dwa wydania genomu człowieka, takie jak GRCh37 i GRCh38, nie różnią się tylko etykietą. Mogą mieć:

  • poprawione sekwencje;
  • zmienione długości i granice regionów;
  • inne reprezentacje obszarów złożonych;
  • dodatkowe alternatywne loci i sekwencje łat;
  • inne położenia tego samego biologicznego motywu.

W UCSC dla głównych ludzkich wydań spotyka się potoczne nazwy hg19 i hg38. W Ensembl używa się nazw zespołu referencyjnego, takich jak GRCh37 i GRCh38. Te pary są blisko powiązane, lecz do pracy automatycznej nie wystarczy bezrefleksyjnie zamienić napis. Trzeba sprawdzić zestaw sekwencji, źródło FASTA oraz sposób obsługi sekwencji dodatkowych.

Zmiana assembly może spowodować, że:

  • pozycja przesunie się o tysiące zasad;
  • odcinek zostanie rozdzielony;
  • region nie będzie miał jednoznacznego odpowiednika;
  • allel referencyjny okaże się inny;
  • gen albo transkrypt otrzyma zmienione granice;
  • odczyty zmapują się inaczej.

Dlatego assembly zapisuje się obok każdej tabeli wariantów, pliku BED, zestawu starterów i wyniku CRISPR. „Współrzędne ludzkie” nie są kompletną metadną.

Chromosom nie wyczerpuje zawartości assembly

Na szkolnym schemacie genom składa się z chromosomów 1–22, X, Y i genomu mitochondrialnego. Referencyjne assembly zawiera jednak także sekwencje potrzebne do przedstawienia regionów, których nie da się uczciwie opisać jedną linią:

  • unlocalized sequences wiadomo przypisać do chromosomu, ale nie do dokładnego miejsca;
  • unplaced sequences należą do assembly, lecz nie zostały przypisane do chromosomu;
  • alternate loci przedstawiają alternatywny haplotyp złożonego regionu;
  • fix patches korygują błąd bez natychmiastowej zmiany głównej wersji;
  • novel patches dodają alternatywną reprezentację.1

Te sekwencje nie są „śmieciowymi contigami”. Mogą zawierać ważne geny i wpływać na jednoznaczność mapowania. Jeśli pipeline wyrzuca wszystko poza chr1–chr22,X,Y, powinien robić to jawnie i z uzasadnieniem.

Przy porównywaniu wyników trzeba zanotować, czy indeks alignera obejmował sekwencje alternatywne, decoye i genom mitochondrialny. Dwie analizy nazwane „GRCh38” mogą używać różnych pakietów FASTA.

1-based inclusive i 0-based half-open

Jednym z najczęstszych źródeł błędów jest mieszanie dwóch konwencji przedziałów.

W widoku graficznym UCSC zapis regionu jest 1-based, domknięty na obu końcach. Pierwsza zasada chromosomu ma pozycję 1, a przedział chr1:1-100 zawiera sto zasad. W formacie BED początek jest natomiast 0-based, a koniec jest wyłączony: chr1 0 100 opisuje ten sam odcinek.2

W konwencji half-open:

długość = end − start

Przykłady:

Obiekt Zapis 1-based inclusive BED 0-based half-open
pierwsza zasada chromosomu 1–1 0 1
sto pierwszych zasad 1–100 0 100
jedna zasada na pozycji 101 101–101 100 101
dziesięć zasad 101–110 101–110 100 110

Przy przejściu z zapisu 1-based inclusive do BED odejmuje się 1 od początku, ale pozostawia koniec. Nie należy odejmować 1 od obu wartości.

Ta sama strona może używać obu konwencji: jednej w pasku pozycji, drugiej w tabeli albo pobranym pliku. Dokumentacja UCSC wprost odróżnia współrzędne wyświetlane od wewnętrznej konwencji baz danych i BED.3 Ensembl w widoku regionu i wielu endpointach operuje na pozycjach 1-based inclusive, lecz każdy format eksportu trzeba sprawdzić osobno.

Praktyczna zasada brzmi: nigdy nie kopiuj liczb bez zapisania formatu źródłowego.

Nić nie zmienia porządku współrzędnych

Gen na nici minus jest transkrybowany w kierunku przeciwnym do rosnących współrzędnych genomowych. Nie zapisuje się jednak jego przedziału jako 900–100. Początek przedziału w pliku nadal jest liczbowo mniejszy od końca, a orientację przechowuje osobne pole + albo .

Ma to kilka konsekwencji:

  • promotor genu na nici plus leży zwykle przed niższą współrzędną początku transkrypcji;
  • promotor genu na nici minus leży w stronę wyższych współrzędnych;
  • „upstream” i „downstream” są pojęciami biologicznymi, a nie synonimami lewo/prawo na ekranie;
  • sekwencję genu na nici minus trzeba odwrócić i skomplementować, aby czytać ją 5′→3′;
  • numeracja pozycji w transkrypcie nie jest tym samym co numeracja genomowa.

To szczególnie ważne przy projektowaniu starterów, RNA przewodnikowego, sond i interpretacji wariantu względem kodonu.

Gen jest zbiorem modeli, a nie jednym prostokątem

Warstwa genowa może narysować kilka transkryptów w jednym locus. Każdy model ma własny zestaw eksonów, granice UTR i często własną sekwencję kodującą.

Warto rozróżniać:

  • gene — locus grupujące powiązane transkrypty;
  • transcript — konkretny model dojrzałego RNA;
  • exon — fragment zachowany w danym transkrypcie;
  • CDS — część transkryptu tłumaczona na białko;
  • UTR — część eksonowa, ale nietłumaczona;
  • intron — odcinek między eksonami w modelu pre-mRNA;
  • protein product — sekwencję aminokwasową wynikającą z wybranego CDS.

Ekson nie jest synonimem fragmentu kodującego. Pierwszy i ostatni ekson często zawierają UTR, a transkrypty niekodujące mogą być w całości eksonowe bez CDS.

Jeżeli wariant jest „w eksonie”, nadal trzeba zapytać: w którym transkrypcie? W jednym może leżeć w CDS, w drugim w 5′ UTR, w trzecim w intronie, a jeszcze inny transkrypt może go nie obejmować.

„Canonical” jest wyborem operacyjnym

Interfejsy próbują ograniczyć przeciążenie informacją i wyróżniają transkrypt reprezentatywny. Etykieta canonical nie oznacza jednak:

  • jedynego biologicznie prawdziwego transkryptu;
  • najczęstszego transkryptu w każdej tkance;
  • izoformy odpowiedzialnej za każdą chorobę;
  • zawsze najdłuższego białka;
  • transkryptu najlepszego dla wszystkich analiz.

Kryteria mogą obejmować konserwację, ekspresję, zgodność z innymi zasobami, kompletność CDS i dowody eksperymentalne. Kryteria oraz wyniki zmieniają się wraz z wydaniem anotacji.

Do raportu klinicznego lub funkcjonalnego lepiej podać dokładny identyfikator transkryptu z wersją. Jeżeli użyto MANE Select, RefSeq Select albo innego zestawu, należy zapisać również nazwę tego standardu i wydanie. Sam symbol genu nie wystarcza.

Stabilny identyfikator nie oznacza niezmiennego modelu

Ensembl nadaje stabilne identyfikatory genom, transkryptom, białkom i eksonom. Do identyfikatora można dołączyć numer wersji po kropce. Reguły wersjonowania są zależne od poziomu obiektu:

  • wersja genu rośnie, gdy zmienia się zestaw przypisanych transkryptów;
  • wersja transkryptu rośnie przy zmianie splicingu, położenia lub sekwencji cDNA;
  • wersja białka rośnie, gdy zmienia się sekwencja peptydu;
  • wersja eksonu rośnie przy zmianie jego sekwencji genomowej.4

Zmiana powiązanej nazwy, opisu albo zewnętrznego odsyłacza nie musi zwiększyć wersji. Dlatego identyfikator z wersją przypina model sekwencyjny lepiej niż sam rdzeń, ale nie zamraża wszystkich metadanych.

Symbol genu, np. TP53, jest wygodny dla człowieka, lecz może być niejednoznaczny między gatunkami, zmieniać się historycznie albo występować w nazwach pseudogenów. W pliku wynikowym warto przechowywać równolegle:

  • gatunek;
  • symbol;
  • stabilny identyfikator genu;
  • stabilny identyfikator transkryptu z wersją;
  • identyfikator białka, jeśli analizowane jest białko;
  • release anotacji;
  • assembly.

Dlaczego warstwy genowe się różnią

W UCSC można zestawić np. modele NCBI RefSeq, GENCODE i inne warstwy genowe. Ensembl prezentuje własny build anotacji oraz powiązania z zewnętrznymi zasobami. Różnice mogą wynikać z:

  • innego zestawu danych wejściowych;
  • odmiennych algorytmów budowania modeli;
  • innego progu dowodowego;
  • kuracji ręcznej;
  • daty wydania;
  • sposobu reprezentacji pseudogenów i transkryptów niekodujących;
  • traktowania częściowych lub problematycznych modeli.

Rozbieżność nie musi być błędem. Może ujawniać realną niepewność. Najlepsze pytanie nie brzmi „która przeglądarka ma rację?”, lecz:

Który model jest właściwy dla mojego pytania i jakie dowody go wspierają?

Do projektowania konstruktu potrzebna jest konkretna sekwencja transkryptu. Do zliczania RNA-seq potrzebny jest zamrożony zestaw anotacji. Do interpretacji klinicznej potrzebne są standardy właściwe dla raportowania wariantów. Każde zadanie może uzasadniać inny wybór.

Warstwa jest wynikiem procesu

Kliknięcie nazwy ścieżki powinno prowadzić do strony opisowej. Przed interpretacją trzeba znaleźć:

  • właściciela i źródło danych;
  • organizm oraz assembly;
  • datę albo release;
  • metodę eksperymentalną;
  • pipeline obliczeniowy;
  • typ obiektu;
  • znaczenie wartości i kolorów;
  • filtry jakości;
  • ograniczenia;
  • sposób pobrania danych.

Ten nawyk chroni przed pomyleniem sygnału z interpretacją. Ścieżka Conservation może łączyć uliniowienia wielu gatunków i wynik konkretnego modelu. Ścieżka wariantów może zawierać zgłoszenia o bardzo różnej jakości. Warstwa regulatorowa może być modelem powstałym z wielu eksperymentów, a nie bezpośrednim odczytem pojedynczej próbki.

Sposób rysowania nie zmienia danych

UCSC pozwala wyświetlać ścieżkę w trybach takich jak hide, dense, squish, pack i full. Zmieniają one zagęszczenie oraz etykiety, a nie zawartość źródłowej tabeli.

To ważne, ponieważ widok dense może scalić wiele elementów w jeden pasek. Zrzut ekranu wygląda wtedy jak jeden obiekt, choć tabela zawiera setki rekordów. Przy dużym oddaleniu przeglądarka może też agregować sygnał.

Podobnie kolor nie ma uniwersalnego znaczenia. Może kodować:

  • nić;
  • typ transkryptu;
  • wynik jakości;
  • poziom sygnału;
  • klasę funkcjonalną;
  • ustawienie autora własnej ścieżki.

Legendę czyta się dla konkretnej warstwy. Czerwony nie zawsze oznacza patogeniczność, a ciemniejszy prostokąt nie zawsze silniejszy dowód.

Warstwy regulatorowe: hipotezy z kontekstem komórkowym

Anotacje regulacyjne integrują dane takie jak ATAC-seq, DNase-seq, ChIP-seq czynników transkrypcyjnych, modyfikacje histonów i dostępne informacje o aktywności w tkankach. Ensembl opisuje regulatory build jako wynik analizy pierwotnych danych, z której powstają cechy regulacyjne i ich anotacja.5

Z takiej warstwy można wnioskować, że region:

  • ma dostępny chromatynowo sygnał w określonym materiale;
  • wiąże badane białko w danym doświadczeniu;
  • nosi modyfikację skorelowaną z klasą elementów;
  • został sklasyfikowany przez pipeline jako promotor, enhancer albo inny element;
  • wykazuje aktywność w wybranych typach komórek.

Nie można bez dodatkowego eksperymentu stwierdzić, że:

  • element jest aktywny we wszystkich tkankach;
  • kontroluje najbliższy gen;
  • pojedynczy pik jest przyczyną ekspresji;
  • wariant w piku zmienia fenotyp;
  • etykieta computational enhancer jest bezpośrednim dowodem mechanizmu.

Odległość wzdłuż chromosomu pomaga, ale regulacja działa w trzech wymiarach. Enhancer może oddziaływać z dalszym promotorem, omijając bliższy gen. Wnioski wzmacniają dane kontaktów chromatyny, eQTL, perturbacje CRISPR, reportery, ekspresja i zgodność kontekstu komórkowego.

Wariant jest relacją do referencji

Wariant powinien być opisany co najmniej przez:

  • assembly;
  • sekwencję/chromosom;
  • pozycję lub przedział;
  • allel referencyjny;
  • allel alternatywny;
  • sposób normalizacji;
  • źródło i wersję anotacji.

Sama pozycja nie wystarcza. Po zmianie assembly można otrzymać tę samą liczbę w innym kontekście, a przy insercjach i delecjach kilka równoważnych zapisów zależnych od wyrównania.

Warstwy dbSNP, ClinVar, gnomAD czy katalogów nowotworowych odpowiadają na inne pytania. Obecność wariantu w dbSNP nie oznacza, że jest częsty ani łagodny. Obecność w ClinVar nie oznacza jednomyślnej klasyfikacji. Częstość populacyjna nie dowodzi braku wpływu w każdym modelu dziedziczenia.

Przeglądarka dobrze nadaje się do sprawdzenia, czy wariant:

  • leży w CDS, UTR, intronie czy poza transkryptem;
  • przecina miejsce splicingowe;
  • wypada w powtórzeniu albo regionie słabo mapowalnym;
  • sąsiaduje z dowodami regulacyjnymi;
  • jest konserwowany;
  • ma różne konsekwencje w kilku transkryptach.

Nie zastępuje jednak wyspecjalizowanej anotacji wariantów i kryteriów interpretacyjnych.

Powtórzenia i mapowalność są częścią interpretacji

Region może wyglądać interesująco biologicznie, a jednocześnie być technicznie trudny. Ścieżki powtórzeń, segmental duplications i mapowalności pomagają ocenić, czy:

  • krótkie odczyty można przypisać jednoznacznie;
  • starter PCR ma unikalne miejsce wiązania;
  • sonda może hybrydyzować w kilku loci;
  • pik ChIP-seq albo ATAC-seq jest podatny na artefakt mapowania;
  • wariant pochodzi z właściwego paralogu;
  • delecja lub duplikacja leży w regionie sprzyjającym rekombinacji.

Brak odczytów w regionie powtórzonym nie zawsze oznacza delecję. Nadmiar odczytów nie zawsze oznacza amplifikację. Najpierw trzeba znać długość odczytów, reguły alignera, jakość mapowania i zastosowane filtry.

Konserwacja nie jest testem funkcji

Ścieżki comparative genomics zestawiają uliniowienia między gatunkami i wyniki mierzące zachowanie sekwencji. Silna konserwacja może wskazywać selekcję i pomagać priorytetyzować elementy funkcjonalne. Nie jest jednak dowodem, że region pełni dokładnie tę samą rolę w badanym typie komórki.

Brak głębokiej konserwacji także nie wyklucza funkcji:

  • element może być ewolucyjnie młody;
  • szybka ewolucja może zamazać podobieństwo;
  • uliniowienie może być niepewne;
  • funkcja może zależeć od ogólnej właściwości, nie dokładnych zasad;
  • region może być charakterystyczny dla jednej linii ewolucyjnej.

Warto oddzielić konserwację sekwencji, konserwację położenia, konserwację ekspresji i konserwację funkcji.

Ortolog nie znaczy „ten sam gen”

Ensembl Compara buduje drzewa genów i relacje ortologii oraz paralogii. Dla części par ocenia pewność z użyciem m.in. zgodności otoczenia genowego, uliniowień genomowych, identyczności sekwencji i zgodności drzewa.6

Pojęcia:

  • ortologi rozdzieliły się wskutek specjacji;
  • paralogi powstały wskutek duplikacji;
  • relacja one-to-one oznacza brak dodatkowych zachowanych duplikacji w rozważanych liniach;
  • relacje one-to-many i many-to-many wymagają ostrożniejszego przenoszenia funkcji.

Wysoki procent identyczności nie wystarcza do rozpoznania ortologu. Po duplikacji dwa paralogi mogą być bardzo podobne, a prawdziwy ortolog w innym gatunku może ewoluować szybciej. Przy planowaniu eksperymentu międzygatunkowego trzeba sprawdzić drzewo, synteny, domeny, ekspresję i literaturę funkcjonalną.

Ensembl: od obiektu biologicznego do sąsiedztwa

Ensembl jest szczególnie wygodny, gdy punkt startowy stanowi gen, transkrypt, wariant albo relacja porównawcza. Typowy tok pracy:

  1. wyszukaj symbol, stabilny identyfikator albo region;
  2. potwierdź gatunek i assembly;
  3. wejdź na stronę genu;
  4. przejrzyj listę transkryptów i ich biotypy;
  5. wybierz konkretny transkrypt;
  6. sprawdź eksony, CDS, UTR, domeny białka i warianty;
  7. przejdź do widoku regionu, aby zobaczyć sąsiedztwo;
  8. włącz dane regulacyjne i comparative genomics;
  9. wyeksportuj sekwencję albo tabelę;
  10. zapisz stabilne identyfikatory, wersje i release.

Widok genu i widok regionu odpowiadają na różne pytania. Strona genu porządkuje obiekty anotacji. Region in detail układa ścieżki na osi genomu. Do pełnej analizy zwykle potrzebne są oba.

UCSC: od regionu do zestawu warstw

UCSC Genome Browser często zaczyna się od pozycji. Jego siłą jest szybkie zestawianie wielu ścieżek, własnych danych i dużych kolekcji udostępnionych przez track hubs.

Praktyczny tok:

  1. wybierz organizm i assembly;
  2. wpisz region, nazwę genu albo identyfikator;
  3. ustal właściwą warstwę genową;
  4. ustaw czytelny tryb widoczności;
  5. otwórz opis każdej istotnej ścieżki;
  6. sprawdź źródło tabeli i znaczenie pól;
  7. dodaj warianty, regulację, repeats i conservation;
  8. załaduj własny BED, bigBed, bedGraph lub bigWig;
  9. użyj Table Browser albo Data Integrator do eksportu;
  10. zapisz sesję i osobny manifest parametrów.

UCSC udostępnia dokumentację custom tracks, track hubs, sesji, narzędzi command-line i usług programowych.7 To sprawia, że przeglądarka może być punktem łączącym oglądanie z analizą, o ile nie traktuje się zrzutu ekranu jako jedynego zapisu.

Custom track i track hub

Custom track jest dobry dla niewielkiego zestawu danych przesyłanego bezpośrednio lub wskazanego URL-em. Może zawierać np. przedziały BED, sygnał bedGraph albo odnośnik do indeksowanego pliku binarnego.

Track hub organizuje wiele zdalnych ścieżek i ich konfigurację. Nadaje się dla projektu, konsorcjum albo laboratorium, które chce publikować duże pliki bigBed i bigWig bez kopiowania ich na serwery UCSC.

Formaty:

  • BED — przedziały i opcjonalne atrybuty;
  • bedGraph — wartość numeryczna dla przedziału;
  • bigBed — indeksowana binarna postać danych przedziałowych;
  • bigWig — indeksowany sygnał ciągły;
  • BAM/CRAM — uliniowienia;
  • VCF/tabix — warianty z indeksem.

Duże dane powinny pozostać na serwerze obsługującym żądania zakresowe. Przeglądarka pobiera wtedy tylko fragment potrzebny dla bieżącego widoku.

Trzeba pamiętać o prywatności. URL zawierający token, dane kliniczne albo niezanonimizowane identyfikatory może trafić do historii, logów serwera i zapisanej sesji. Danych wrażliwych nie publikuje się przez przypadkowy publiczny hosting.

Table Browser i Data Integrator

Widok graficzny pomaga myśleć, ale wynik do dalszej analizy powinien być tabelą. UCSC udostępnia kilka dróg:

  • Table Browser — zapytanie do pojedynczej tabeli z filtrami i przecięciami;
  • Data Integrator — łączenie danych z kilku ścieżek dla wspólnych regionów;
  • REST API — programowy dostęp do sekwencji i anotacji;
  • publiczne pliki i narzędzia wiersza poleceń;
  • bezpośredni dostęp do części tabel bazowych.8

Przed pobraniem trzeba znać tabelę stojącą za ścieżką. Nazwa widoczna w interfejsie może być grupą kilku tabel albo composite track. Eksport domyślny nie zawsze odpowiada temu, co użytkownik wyobraził sobie na podstawie rysunku.

W wynikach warto zachować:

  • pełne zapytanie;
  • assembly;
  • nazwę tabeli;
  • filtr;
  • zakres;
  • format;
  • datę;
  • liczbę rekordów;
  • checksumę pobranego pliku.

BioMart, REST i API Ensembl

Ensembl daje kilka poziomów dostępu:

  • eksport z widoku obiektu lub regionu;
  • BioMart do budowania tabel genów, transkryptów i atrybutów;
  • REST API do zapytań po identyfikatorach, sekwencji, regionach i porównaniach;
  • programistyczne API i pliki FTP dla większych analiz.

Endpoint overlap potrafi zwrócić geny, transkrypty, eksony, CDS, warianty lub cechy regulacyjne przecinające region.9 To wygodne, lecz wymaga jawnego określenia gatunku, regionu, typu cechy i wersji usługi.

Przykładowy logiczny kontrakt zapytania:

gatunek: homo_sapiens
assembly: GRCh38
region: 7:140453136-140453145
feature: transcript
strand: dowolna
release/data: zapisane w manifeście

Odpowiedź JSON nie jest samowyjaśniająca. Trzeba zachować kod, parametry, nagłówki lub wersję dokumentacji i surową odpowiedź. Zapytanie wykonane za rok może zwrócić nowszą anotację.

Przeglądarka do odkrywania, API do odtwarzania

Najefektywniejszy model pracy ma dwa etapy:

Oglądanie

  • znajdź interesujący locus;
  • włącz potrzebne warstwy;
  • zauważ konflikty i nietypowe cechy;
  • sformułuj pytanie, które da się zapisać jako regułę.

Analiza odtwarzalna

  • pobierz źródłowe rekordy;
  • zapisz wersje;
  • zastosuj jawne filtry;
  • policz wynik skryptem;
  • wygeneruj tabelę i rycinę;
  • zachowaj manifest.

Ręczne kliknięcie dziesięciu genów jest dobrą eksploracją. Ręczne kliknięcie dziesięciu tysięcy genów nie jest pipeline’em.

LiftOver nie jest zmianą etykiety

UCSC liftOver i analogiczne narzędzia Ensembl mapują współrzędne między assembly na podstawie przygotowanych relacji sekwencyjnych. To nie jest proste dodanie stałej do pozycji.

Odcinek może:

  • zmapować się jednoznacznie;
  • zmapować się na nici przeciwnej;
  • zmienić długość;
  • rozdzielić się na kilka fragmentów;
  • mieć kilka możliwych miejsc;
  • nie mieć odpowiednika.

Przyczyny obejmują poprawki assembly, luki, rearrangements, alternatywne loci i regiony powtórzone.

Dla pojedynczego wariantu po liftOver trzeba sprawdzić:

  1. czy nowa pozycja rzeczywiście istnieje;
  2. czy allel referencyjny zgadza się z nowym FASTA;
  3. czy orientacja nie wymaga komplementacji alleli;
  4. czy wariant należy ponownie znormalizować;
  5. czy anotacja transkryptowa została policzona od nowa.

Przeniesiona pozycja nie przenosi automatycznie starej konsekwencji biologicznej. Jeśli granica eksonu lub model transkryptu zmienił się w nowym wydaniu, konsekwencję trzeba ponownie anotować.

Dla regionów regulacyjnych i CNV wynik liftOver jest jeszcze bardziej delikatny. Częściowe mapowanie dużego przedziału nie oznacza, że cały element biologiczny ma odpowiednik.

Kontrola liftOver krok po kroku

Bezpieczny przebieg:

  1. zachowaj oryginalny plik i jego assembly;
  2. upewnij się, że format przedziału jest poprawny;
  3. zastosuj właściwy plik chain albo usługę mapowania;
  4. zachowaj osobno rekordy zmapowane i niezmapowane;
  5. policz odsetek oraz rozkład zmian długości;
  6. skontroluj wielokrotne i rozdzielone mapowania;
  7. dla wariantów porównaj allele z nową referencją;
  8. ponownie uruchom anotację;
  9. zapisz wersję narzędzia i chain;
  10. przejrzyj wzrokowo reprezentatywne oraz krytyczne przypadki.

Jeżeli tysiąc regionów wejściowych dało 997 regionów wyjściowych, brak trzech jest wynikiem, nie drobnym komunikatem do usunięcia.

Przykład: od symbolu genu do wariantu

Załóżmy, że publikacja podaje symbol genu i zmianę aminokwasową.

Krok 1: identyfikacja

W Ensembl wyszukaj gen, potwierdź gatunek i sprawdź wszystkie transkrypty. Zanotuj stabilny identyfikator genu oraz transkrypt, względem którego opisano białko.

Krok 2: model

Sprawdź wersję transkryptu, CDS, orientację i pozycję kodonu. Porównaj z RefSeq albo wskazanym standardem raportowania.

Krok 3: region

Przejdź do widoku genomowego. Zobacz, czy pozycja leży w unikalnym regionie, w powtórzeniu albo blisko granicy eksonu.

Krok 4: warianty

Włącz właściwe warstwy wariantowe. Otwórz rekord, zamiast interpretować sam kolor. Sprawdź allele, assembly, częstości i zgłoszenia.

Krok 5: regulacja i konserwacja

Jeśli wariant niekodujący, sprawdź dowody właściwe dla tkanki. Konserwację potraktuj jako dodatkową informację, nie werdykt.

Krok 6: odtwarzalność

Wyeksportuj sekwencję i rekordy. Zapisz release, ustawienia warstw, region, datę i identyfikatory. Zrzut ekranu dołącz jako ilustrację, nie jako jedyne źródło danych.

Przykład: własne piki ATAC-seq

Plik z pikami ATAC-seq powinien mieć jawnie przypisane assembly. Przed załadowaniem:

  • sprawdź, czy nazwy chromosomów odpowiadają przeglądarce;
  • potwierdź konwencję BED;
  • usuń lub wyjaśnij niedozwolone sekwencje;
  • zachowaj wynik kontroli jakości;
  • ustal, czy piki pochodzą z replik połączonych, consensus czy jednego eksperymentu.

Po wyświetleniu nie wystarczy znaleźć pik nad promotorem. Trzeba zapytać:

  • czy sygnał występuje we wszystkich replikach;
  • czy region jest mapowalny;
  • czy pik przecina TSS konkretnego transkryptu;
  • czy komórki referencyjnej warstwy regulacyjnej są porównywalne;
  • czy zmiana dostępności wiąże się ze zmianą ekspresji;
  • czy wynik jest specyficzny względem kontroli.

Przeglądarka podpowiada relacje przestrzenne. Statystykę i kontrolę projektu wykonuje się poza nią.

Najczęstsze pułapki

Objaw Możliwa przyczyna Kontrola
region przesunięty o jedną zasadę pomieszanie 1-based i 0-based przetestuj jednoelementowy przedział
brak danych na wszystkich chromosomach chr1 kontra 1 porównaj nagłówki FASTA i nazwy sekwencji
wariant ma inny allel referencyjny złe assembly lub nić odczytaj bazę z dokładnego FASTA
gen ma inne granice inny release lub zestaw anotacji porównaj identyfikatory i wersje transkryptów
pik zniknął po przybliżeniu agregacja albo ustawienia widoku pobierz rekordy źródłowe
setki wyników dla jednego genu wiele transkryptów i aliasów określ poziom obiektu i deduplikację
liftOver zgubił region brak jednoznacznego mapowania zachowaj plik unmapped i sprawdź chain
pozycja jest poprawna, konsekwencja inna zmieniona anotacja ponownie anotuj na docelowym release
„enhancer” nie działa w doświadczeniu inny kontekst komórkowy lub modelowa etykieta sprawdź źródła i wykonaj test funkcjonalny
gen ma kilka „najważniejszych” izoform różne kryteria canonical/MANE/ekspresji wybierz transkrypt właściwy dla pytania

Minimalny test pliku BED

Przed analizą dużego pliku warto przygotować trzy kontrolne rekordy:

chr1    0      1      first_base
chr1    99     100    base_100
chr1    100    110    ten_bases

W przeglądarce powinny odpowiadać pozycjom 1, 100 i 101–110. Jeśli nie odpowiadają, problem trzeba rozwiązać przed załadowaniem miliona rekordów.

Dodatkowe kontrole:

  • każda linia ma start < end;
  • współrzędne nie przekraczają długości sekwencji;
  • wartości start nie są ujemne;
  • nazwy sekwencji należą do assembly;
  • strand, jeśli występuje, ma dozwoloną wartość;
  • pola liczbowe nie zawierają separatorów tysięcy;
  • plik nie miesza dwóch assemblies.

Co zachować z sesji przeglądarki

Zapisana sesja UCSC lub odsyłacz Ensembl pomaga odtworzyć widok, ale zewnętrzne track hubs i bieżące bazy mogą się zmieniać. Solidny pakiet wynikowy obejmuje:

  • organizm;
  • accession i nazwę assembly;
  • datę dostępu;
  • release Ensembl/UCSC albo wersje źródłowych tracków;
  • dokładny region i konwencję;
  • listę włączonych ścieżek;
  • tryby widoczności i ważne filtry;
  • identyfikatory genów oraz transkryptów z wersjami;
  • URL lub nazwę tabeli źródłowej;
  • pliki eksportu z checksumami;
  • skrypt pobierający lub przetwarzający dane;
  • zrzut ekranu z legendą;
  • notatkę o liftOver i pliku chain, jeśli wykonano mapowanie.

Najlepiej dodać krótki README, w którym jedna osoba spoza projektu może znaleźć odpowiedź: „co dokładnie oglądam i jak mogę pobrać to jeszcze raz?”.

Jak cytować przeglądarkę i jej dane

Nie ma jednego uniwersalnego cytowania dla całego widoku. Należy rozdzielić:

  • publikację lub stronę głównego zasobu;
  • assembly;
  • anotację genową;
  • konkretną ścieżkę danych;
  • własny pipeline;
  • datę dostępu, jeśli zawartość jest aktualizowana.

Jeżeli wniosek opiera się na warstwie konsorcjum, trzeba cytować źródło tej warstwy, nie tylko UCSC jako interfejs. Jeśli Ensembl zintegrował zewnętrzny katalog, opis ścieżki wskaże źródło pierwotne.

W metodach można zapisać:

Regiony zdefiniowano na assembly GRCh38 i przecięto z tabelą X w wersji Y. Modele transkryptów pochodziły z Ensembl release Z. Współrzędne przechowywano jako BED 0-based half-open; do prezentacji konwertowano je na 1-based inclusive.

To jest bardziej informacyjne niż „geny sprawdzono w Genome Browser”.

Kiedy użyć którego narzędzia

Ensembl jest naturalnym pierwszym wyborem, gdy:

  • zaczynasz od genu lub stabilnego identyfikatora;
  • porównujesz transkrypty;
  • śledzisz relacje gen–transkrypt–białko;
  • potrzebujesz danych Compara, ortologii i drzew genów;
  • chcesz użyć BioMart lub REST dla obiektów anotacji.

UCSC jest naturalnym pierwszym wyborem, gdy:

  • zaczynasz od regionu;
  • chcesz szybko nałożyć wiele niezależnych tracków;
  • masz własne BED, bigWig, BAM albo VCF;
  • korzystasz z track hubs;
  • potrzebujesz Table Browser, Data Integrator albo narzędzi do formatów.

W praktyce najlepszy workflow często używa obu. Ważniejsze od lojalności wobec interfejsu jest zachowanie tej samej osi odniesienia oraz sprawdzenie pochodzenia warstw.

Dziesięć pytań przed wyciągnięciem wniosku

  1. Jakie assembly jest wybrane?
  2. Czy nazwa sekwencji oznacza chromosom główny, patch czy alternate locus?
  3. Czy współrzędne są 1-based inclusive, czy 0-based half-open?
  4. Czy nić ma znaczenie dla pytania?
  5. Jaki model genu i wersja transkryptu są używane?
  6. Z jakiego źródła i wydania pochodzi ścieżka?
  7. Czy oglądam pomiar, wynik pipeline’u, czy interpretacyjną klasę?
  8. Czy tryb widoku agreguje albo ukrywa rekordy?
  9. Czy region jest powtórzony lub słabo mapowalny?
  10. Czy wynik da się odtworzyć z zapisanych plików i parametrów?

Jeśli na te pytania można odpowiedzieć, przeglądarka staje się narzędziem analitycznym. Jeśli nie, pozostaje atrakcyjnym obrazkiem.

Najkrótsza reguła pracy

Współrzędna ma sens dopiero jako:

gatunek + assembly + sekwencja + konwencja przedziału + nić

Anotacja ma sens dopiero jako:

źródło + release + identyfikator obiektu + wersja

A obraz z przeglądarki ma sens dopiero jako:

zestaw warstw + ich pochodzenie + ustawienia + pliki źródłowe

Ensembl i UCSC nie usuwają niepewności z biologii. Robią coś bardziej użytecznego: pozwalają ułożyć różne rodzaje dowodów na wspólnej osi i sprawdzić, gdzie są zgodne, gdzie się różnią oraz czego jeszcze nie wiadomo.

Przypisy