Un ghid pas cu pas pentru analiza datelor de secvențiere a ADN-ului

Dr. Javier Quilez Oliete, un consultant experimentat în bioinformatică pe Kolabtree, oferă un ghid cuprinzător pentru analiza datelor de secvențiere a ADN-ului, inclusiv instrumentele și software-ul utilizat pentru a citi datele.

Introducere

Acidul dezoxiribonucleic (ADN) este molecula care poartă cea mai mare parte a informațiilor genetice ale unui organism. (La unele tipuri de virusuri, informația genetică este transportată de acidul ribonucleic (ARN)). Nucleotidele (reprezentate în mod convențional prin literele A, C, G sau T) sunt unitățile de bază ale moleculelor de ADN. Din punct de vedere conceptual, secvențierea ADN este procesul de citire a nucleotidelor care alcătuiesc o moleculă de ADN (de exemplu, „GCAAACCAAT” este un șir de 10 nucleotide de ADN). Tehnologiile actuale de secvențiere produc milioane de astfel de citiri de ADN într-un timp rezonabil și la un cost relativ scăzut. Ca referință, costul secvențierii unui genom uman – un genom este setul complet de molecule de ADN dintr-un organism – a coborât bariera de 100 de dolari și poate fi realizat în câteva zile. Acest lucru contrastează cu prima inițiativă de secvențiere a genomului uman, care a fost finalizată într-un deceniu și a avut un cost de aproximativ 2,7 miliarde de dolari.

Această capacitate de a secvenția ADN la un randament ridicat și la un cost redus a permis dezvoltarea unui număr tot mai mare de metode și aplicații bazate pe secvențiere. De exemplu, secvențierea unor genomuri întregi sau a regiunilor lor codificatoare de proteine (două abordări cunoscute sub numele de secvențiere a genomului întreg și, respectiv, secvențiere a exomului) la persoane bolnave și la persoane sănătoase poate indica alterări ale ADN-ului care cauzează boli. De asemenea, secvențierea ARN-ului care este transcris din ADN – o tehnică cunoscută sub numele de secvențiere a ARN-ului – este utilizată pentru a cuantifica activitatea genelor și modul în care aceasta se modifică în diferite condiții (de exemplu, netratată față de tratament). Pe de altă parte, metodele de secvențiere a capturii conformației cromozomilor detectează interacțiunile dintre moleculele de ADN apropiate și ajută astfel la determinarea distribuției spațiale a cromozomilor în interiorul celulei.

Comună cu aceste și alte aplicații ale secvențierii ADN este generarea unor seturi de date de ordinul gigabytes și care cuprind milioane de secvențe citite. Prin urmare, a da sens experimentelor de secvențiere de mare randament (HTS) necesită capacități substanțiale de analiză a datelor. Din fericire, există instrumente computaționale și statistice dedicate și fluxuri de lucru de analiză relativ standard pentru majoritatea tipurilor de date HTS. În timp ce unele dintre etapele (inițiale) de analiză sunt comune pentru majoritatea tipurilor de date de secvențiere, o analiză mai în aval va depinde de tipul de date și/sau de scopul final al analizei. Mai jos, ofer o introducere în etapele fundamentale ale analizei datelor HTS și fac referire la instrumente populare.

Câteva dintre secțiunile de mai jos se concentrează pe analiza datelor generate de tehnologiile de secvențiere cu citire scurtă (în principal Illumina), deoarece acestea au dominat în mod istoric piața HTS. Cu toate acestea, tehnologiile mai noi care generează citiri mai lungi (de exemplu, Oxford Nanopore Technologies, PacBio) câștigă rapid teren. Deoarece secvențierea cu citire lungă are unele particularități (de exemplu, rate de eroare mai mari), se dezvoltă instrumente specifice pentru analiza acestui tip de date.

Controlul de calitate (QC) al citirilor brute

Analiștii dornici vor începe analiza pornind de la fișiere FASTQ; formatul FASTQ este de mult timp standardul pentru stocarea datelor de secvențiere cu citire scurtă. În esență, fișierele FASTQ conțin secvența de nucleotide și calitatea de apelare per-bază pentru milioane de citiri. Deși dimensiunea fișierului va depinde de numărul real de citiri, fișierele FASTQ sunt de obicei mari (de ordinul megabytes și gigabytes) și sunt comprimate. De menționat că majoritatea instrumentelor care utilizează fișiere FASTQ ca intrare le pot gestiona în format comprimat, astfel încât, pentru a economisi spațiu pe disc, se recomandă să nu le decomprimați. Ca o convenție, aici voi echivala un fișier FASTQ cu un eșantion de secvențiere.

FastQC este probabil cel mai popular instrument pentru a efectua QC al citirilor brute. Acesta poate fi rulat prin intermediul unei interfețe vizuale sau în mod programatic. În timp ce prima opțiune poate fi mai convenabilă pentru utilizatorii care nu se simt confortabil cu mediul liniei de comandă, cea de-a doua oferă o scalabilitate și o reproductibilitate incomparabile (gândiți-vă cât de anevoioasă și predispusă la erori poate fi rularea manuală a instrumentului pentru zeci de fișiere). Oricum ar fi, principala ieșire a FastQC este un fișier HTML care raportează statistici rezumative cheie despre calitatea generală a citirilor de secvențiere brute de la o anumită probă. Inspectarea zecilor de rapoarte FastQC unul câte unul este plictisitoare și complică comparația între eșantioane. Prin urmare, este posibil să doriți să utilizați MultiQC, care agregă rapoartele HTML de la FastQC (precum și de la alte instrumente utilizate în aval, de exemplu, tăierea adaptoarelor, alinierea) într-un singur raport.

Informațiile QC au rolul de a permite utilizatorului să aprecieze dacă eșantioanele au o calitate bună și, prin urmare, pot fi utilizate pentru etapele ulterioare sau trebuie să fie eliminate. Din nefericire, nu există un prag consensual bazat pe metricile FastQC pentru a clasifica eșantioanele ca fiind de calitate bună sau proastă. Abordarea pe care o folosesc eu este următoarea. Mă aștept ca toate eșantioanele care au trecut prin aceeași procedură (de exemplu, extracția ADN, pregătirea bibliotecii) să aibă statistici de calitate similare și o majoritate de indicatori „pass”. Dacă unele eșantioane au o calitate mai mică decât media, le voi folosi totuși în analiza din aval ținând cont de acest lucru. Pe de altă parte, dacă toate eșantioanele din experiment primesc în mod sistematic stegulețe „warning” sau „fail” în mai multe măsurători (a se vedea acest exemplu), suspectez că ceva nu a mers bine în experiment (de exemplu, calitatea proastă a ADN-ului, pregătirea bibliotecii etc.) și recomand repetarea acestuia.

Read trimming

QC of raw reads ajută la identificarea eșantioanelor problematice, dar nu îmbunătățește calitatea reală a lecturilor. Pentru a face acest lucru, trebuie să tăiem citirile pentru a elimina secvențele tehnice și capetele de calitate scăzută.

Secvențele tehnice sunt resturi din procedura experimentală (de exemplu, adaptoare de secvențiere). Dacă astfel de secvențe sunt adiacente la adevărata secvență a citirii, alinierea (a se vedea mai jos) poate cartografia citirile la o poziție greșită în genom sau poate scădea încrederea într-o anumită aliniere. Pe lângă secvențele tehnice, este posibil să dorim, de asemenea, să eliminăm secvențele de origine biologică dacă acestea sunt foarte prezente printre citiri. De exemplu, procedurile suboptime de preparare a ADN-ului pot lăsa în eșantion o proporție mare de ARN ribozomal (ARNr) transformat în ADN. Cu excepția cazului în care acest tip de acid nucleic este ținta experimentului de secvențiere, păstrarea citirilor derivate din ARNr nu va face decât să mărească sarcina de calcul a etapelor din aval și poate confunda rezultatele. De remarcat, dacă nivelurile de secvențe tehnice, ARNr sau alți contaminanți sunt foarte ridicate, care probabil vor fi fost deja evidențiate de QC, este posibil să doriți să aruncați întreaga probă de secvențiere.

În secvențierea cu citire scurtă, secvența de ADN este determinată nucleotidă cu nucleotidă (din punct de vedere tehnic, o nucleotidă la fiecare ciclu de secvențiere). Cu alte cuvinte, numărul de cicluri de secvențiere determină lungimea de citire. O problemă cunoscută a metodelor de secvențiere HTS este diminuarea acurateței cu care sunt determinate nucleotidele pe măsură ce se acumulează ciclurile de secvențiere. Acest lucru se reflectă într-o scădere generală a calității de apelare per bază, în special spre sfârșitul citirii. Așa cum se întâmplă în cazul secvențelor tehnice, încercarea de a alinia lecturi care conțin capete de calitate scăzută poate duce la erori de poziționare sau la o calitate slabă a cartografierii.

Pentru a elimina secvențele tehnice/contaminante și capetele de calitate scăzută, există și sunt utilizate pe scară largă instrumente de tăiere a lecturilor precum Trimmomatic și Cutadapt. În esență, astfel de instrumente vor elimina secvențele tehnice (disponibile la nivel intern și/sau furnizate de utilizator) și vor tăia citirile pe baza calității, maximizând în același timp lungimea citirilor. Lecturile care rămân prea scurte după tăiere sunt eliminate (lecturile excesiv de scurte, de exemplu <36 nucleotide, complică etapa de aliniere, deoarece acestea se vor mapa probabil la mai multe situri din genom). Este posibil să doriți să vă uitați la procentul de lecturi care supraviețuiesc tăierii, deoarece o rată ridicată de lecturi eliminate este probabil un semn de date de proastă calitate.

În cele din urmă, de obicei, reiau FastQC pe citirile tăiate pentru a verifica dacă această etapă a fost eficientă și a îmbunătățit sistematic metricile QC.

Alignment

Cu excepții (de exemplu, asamblarea de novo), alinierea (denumită și cartografiere) este, de obicei, următoarea etapă pentru majoritatea tipurilor de date și aplicațiilor HTS. Alinierea citirilor constă în determinarea poziției în genom din care a derivat secvența citirii (exprimată de obicei ca cromozom:start-end). Prin urmare, în această etapă este necesară utilizarea unei secvențe de referință pentru alinierea/cartarea citirilor.

alegerea secvenței de referință va fi determinată de mai mulți factori. Pentru unul, specia din care provine ADN-ul secvențiat. Deși numărul de specii pentru care este disponibilă o secvență de referință de înaltă calitate este în creștere, este posibil ca acest lucru să nu fie încă valabil pentru unele organisme mai puțin studiate. În aceste cazuri, este posibil să doriți să aliniați citirile la o specie apropiată din punct de vedere evolutiv pentru care este disponibil un genom de referință. De exemplu, deoarece nu există o secvență de referință pentru genomul coiotului, îl putem folosi pe cel al câinelui, care este foarte apropiat, pentru alinierea citirilor. În mod similar, este posibil să dorim totuși să ne aliniem citirile la o specie înrudită îndeaproape pentru care există o secvență de referință de calitate superioară. De exemplu, deși a fost publicat genomul gibonului, acesta este fragmentat în mii de fragmente care nu recapitulează pe deplin organizarea acelui genom în zeci de cromozomi; în acest caz, efectuarea alinierii folosind secvența de referință umană poate fi benefică.

Un alt factor de luat în considerare este versiunea ansamblului secvenței de referință, deoarece sunt publicate noi versiuni pe măsură ce secvența este actualizată și îmbunătățită. Este important faptul că coordonatele unei anumite alinieri pot varia de la o versiune la alta. De exemplu, mai multe versiuni ale genomului uman pot fi găsite în UCSC Genome Browser. În orice specie, favorizez cu tărie migrarea la cea mai nouă versiune de asamblare, odată ce aceasta este complet lansată. Acest lucru poate cauza unele neplăceri în timpul tranziției, deoarece rezultatele deja existente vor fi relative la versiunile mai vechi, dar se amortizează pe termen lung.

În plus, tipul de date de secvențiere contează și el. Lecturile generate de protocoalele DNA-seq, ChIP-seq sau Hi-C vor fi aliniate la secvența de referință a genomului. Pe de altă parte, deoarece ARN-ul transcris din ADN este prelucrat ulterior în ARNm (adică introni eliminați), multe lecturi RNA-seq nu se vor alinia la o secvență de referință a genomului. În schimb, trebuie fie să le aliniem la secvențele de referință ale transcriptomului, fie să folosim alinieri care iau în considerare divizarea (a se vedea mai jos) atunci când folosim secvența genomului ca referință. Legat de acest aspect este alegerea sursei pentru adnotarea secvenței de referință, și anume baza de date cu coordonatele genelor, transcriptelor, centromerilor etc. Eu folosesc de obicei adnotarea GENCODE, deoarece combină o adnotare cuprinzătoare a genelor și a secvențelor de transcripte.

A fost dezvoltată o listă lungă de instrumente de aliniere a secvențelor cu citire scurtă (a se vedea secțiunea Short-read sequence alignment aici). Trecerea lor în revistă depășește scopul acestui articol (detalii despre algoritmii din spatele acestor instrumente pot fi găsite aici). Din experiența mea, printre cele mai populare se numără Bowtie2, BWA, HISAT2, Minimap2, STAR și TopHat. Recomandarea mea este să vă alegeți alignerul pe baza luării în considerare a unor factori cheie, cum ar fi tipul de date HTS și aplicația, precum și acceptarea de către comunitate, calitatea documentației și numărul de utilizatori. De exemplu, este nevoie de alinieri precum STAR sau Bowtie2, care sunt conștienți de joncțiunile exon-exon atunci când cartografiază RNA-seq la genom.

Comună pentru majoritatea mapatorilor este nevoia de a indexa secvența utilizată ca referință înainte de a avea loc alinierea efectivă. Acest pas poate consuma mult timp, dar trebuie făcut doar o singură dată pentru fiecare secvență de referință. Majoritatea cartografilor vor stoca alinierile în fișiere SAM/BAM, care respectă formatul SAM/BAM (fișierele BAM sunt versiuni binare ale fișierelor SAM). Alinierea este unul dintre pașii care necesită cel mai mult timp și calcul în analiza datelor de secvențiere, iar fișierele SAM/BAM sunt grele (de ordinul gigabytes). Prin urmare, este important să vă asigurați că dispuneți de resursele necesare (a se vedea secțiunea finală de mai jos) pentru a executa alinierea într-un timp rezonabil și pentru a stoca rezultatele. În mod similar, din cauza dimensiunii și a formatului binar al fișierelor BAM, evitați să le deschideți cu editori de text; în schimb, folosiți comenzi Unix sau instrumente dedicate, cum ar fi SAMtools.

După alinieri

Am spune că nu există o etapă comună clară după aliniere, deoarece în acest punct este punctul în care fiecare tip de date HTS și fiecare aplicație pot fi diferite.

O analiză comună în aval pentru datele ADN-seq este apelarea variantelor, adică identificarea pozițiilor din genom care variază în raport cu genomul de referință și între indivizi. Un cadru de analiză popular pentru această aplicație este GATK pentru polimorfism de un singur nucleotid (SNP) sau inserții/deleții mici (indels) (figura 2). Variantele care cuprind bucăți mai mari de ADN (denumite și variante structurale) necesită metode de apelare dedicate (a se vedea acest articol pentru o comparație cuprinzătoare). Ca și în cazul alinierelor, vă sfătuiesc să selectați instrumentul potrivit luând în considerare factori cheie precum tipul de variante (SNP, indel sau variante structurale), acceptarea de către comunitate, calitatea documentației și numărul de utilizatori.

Probabil cea mai frecventă aplicație a RNA-seq este cuantificarea expresiei genice. Din punct de vedere istoric, citirile trebuiau să fie aliniate la secvența de referință și apoi numărul de citiri aliniate la o anumită genă sau transcripție era utilizat ca un proxy pentru a cuantifica nivelurile de expresie ale acesteia. Această abordare de aliniere+cuantificare este realizată de instrumente precum Cufflinks, RSEM sau featureCounts. Cu toate acestea, abordarea scuh a fost depășită din ce în ce mai mult de metode mai noi, implementate în software precum Kallisto și Salmon. Din punct de vedere conceptual, cu astfel de instrumente, nu este necesară alinierea secvenței complete a unei lecturi la secvența de referință. În schimb, trebuie doar să se alinieze suficiente nucleotide pentru a fi siguri că o citire provine dintr-un anumit transcript. Simplificând, abordarea aliniere+cuantificare este redusă la o singură etapă. Această abordare este cunoscută sub numele de pseudo-mapping și crește considerabil viteza de cuantificare a expresiei genice. Pe de altă parte, rețineți că pseudo-mapping-ul nu va fi potrivit pentru aplicațiile în care este necesară alinierea completă (de exemplu, apelarea variantelor din datele RNA-seq).

Un alt exemplu al diferențelor dintre etapele de analiză în aval și instrumentele necesare între aplicațiile bazate pe secvențiere este ChIP-seq. Citirile generate cu o astfel de tehnică vor fi utilizate pentru apelarea vârfurilor, care constă în detectarea regiunilor din genom cu un exces semnificativ de citiri care indică locul în care este legată proteina țintă. Există mai multe dispozitive de apelare a vârfurilor, iar această publicație le analizează. Ca un ultim exemplu, voi menționa datele Hi-C, în care alinierile sunt utilizate ca date de intrare pentru instrumente care determină matricele de interacțiune și, din acestea, caracteristicile 3D ale genomului. Comentariile cu privire la toate testele bazate pe secvențiere depășesc sfera de aplicare a acestui articol (pentru o listă relativ completă, a se vedea acest articol).

Înainte de a începe…

Partea rămasă a acestui articol abordează aspecte care pot să nu fie strict considerate ca etape în analiza datelor HTS și care sunt în mare parte ignorate. În schimb, eu susțin că este capital să vă gândiți la întrebările puse în Tabelul 1 înainte de a începe analiza datelor HTS (sau a oricărui tip de date, de fapt) și am scris despre aceste subiecte aici și aici.

Tabelul 1

Gândiți-vă	Acțiune propusă
Dispuneți de toate informațiile despre eșantionul dumneavoastră necesare pentru analiză?	Colectați sistematic metadatele experimentelor
Vă veți putea identifica fără echivoc eșantionul dumneavoastră?	Stabiliți un sistem pentru a atribui fiecărui eșantion un identificator unic
Unde vor fi datele și rezultatele?	Organizarea structurată și ierarhică a datelor
Vă veți putea procesa mai multe eșantioane fără probleme?	Scalabilitatea, paralelizarea, configurarea automată și modularitatea codului
Vă veți putea, dumneavoastră sau oricine altcineva, să reproduceți rezultatele?	Documentați-vă codul și procedurile!

Cum s-a menționat mai sus, datele brute HTS și unele dintre fișierele generate în timpul analizei acestora sunt de ordinul gigabytes, astfel încât nu este excepțional ca un proiect care include zeci de probe să necesite terabytes de stocare. În plus, unii pași în analiza datelor HTS sunt intensivi din punct de vedere computațional (de exemplu, alinierea). Cu toate acestea, infrastructura de stocare și de calcul necesară pentru analiza datelor HTS este un aspect important și este adesea trecut cu vederea sau nu este discutat. Ca exemplu, în cadrul unei analize recente, am analizat zeci de lucrări publicate care efectuau analiza de asociere la nivel de fenomene (PheWAS). PheWAS-urile moderne analizează 100-1.000 de variante genetice și fenotipuri, ceea ce duce la o stocare de date și o putere de calcul importante. Și totuși, practic niciunul dintre articolele pe care le-am analizat nu a comentat infrastructura necesară pentru analiza PheWAS. Deloc surprinzător, recomandarea mea este să planificați din timp cerințele de stocare și de calcul cu care vă veți confrunta și să le împărtășiți cu comunitatea.

Aveți nevoie de ajutor cu analiza datelor de secvențiere ADN? Intrați în contact cu bioinformaticieni independenți și experți în genomică pe Kolabtree.