Informație

Simularea ratei de substituție a mutațiilor neutre


Încerc să simulez computerizat o populație bazată pe modelul Wright-Fisher. Aș dori să ajung la rezultatul clasic al teoriei neutre a evoluției moleculare conform căreia rata de substituție neutră este egală cu rata mutației. Cu toate acestea, rezultatele simulărilor mele nu arată niciodată acest lucru. Voi explica pe scurt procesul meu de simulare și cum calculez ratele de mutație și substituție. Vă rog să-mi spuneți dacă există ceva în neregulă în procedura mea.

Am o populație haploidă nestructurată care crește în pași de timp discreți. Încep cu o populație de N indivizi într-un mediu cu capacitate de încărcare K. Folosesc un model de mutație alelă K. La fiecare generație, fiecare individ va reproduce n urmași și va muri imediat după aceea. O fracțiune de descendenți va fi aleasă aleatoriu pentru a avea o singură mutație în fiecare generație (u este rata mutației pe genă pe generație). La fiecare generație, înregistrez câte mutații au apărut și câte mutații au ajuns la fixare în populație (frecvența a atins mai mult de 99% din indivizi). Astfel, rata mutației și rata de substituție pe generație sunt calculate ca numărul de mutații și numărul de substituții împărțit la dimensiunea populației, respectiv. Cu această configurare, chiar și cu sute de repetări, rata de fixare nu va fi niciodată la fel de mare ca rata de mutație u.

Unde crezi că este problema? M-am luptat mult cu asta și apreciez orice ajutor.

Editați: Codul meu este scris în Python. Acesta este pseudo-codul meu:

Parametri: dimensiunea populației, dimensiunea genomului, rata_mutației, numărul_generației, pragul_fixare. "Dimensiunea_populației: aici este 1500 Dimensiunea genomului: numărul de loci care pot fi mutați. Fiecare locus poate fi mutat la oricare dintre ACTG. Aici 10 ^ 4 Rata_mutației: numărul de mutații pe genom pe generație. Aici 0.05 Număr_generație: Numărul de generații pe care populația le crește. Aici 1000 Fixare_prag: Fracțiunea indivizilor care au o alelă în populație, astfel încât alela să fie considerată fixă. Aici 0.99 "Se creează populația inițială. Toți indivizii sunt aceiași fără nicio mutație Se creează un vector care să dețină numărul de mutații și un vector care să dețină numărul de fixații. Pentru fiecare generație: stabilesc câți indivizi din fiecare generație ar trebui să fie mutați. Aici 8. Aleg aleatoriu (uniform cu înlocuitor) indivizii care urmează să fie mutați Pentru fiecare individ din populație: Va produce o descendență. Puii moștenesc genomul părintelui său, adică mutațiile părintelui. Dacă individul se numără printre cei care vor fi mutați, descendenții vor primi o mutație care este aleasă la întâmplare. Mutația poate fi o mutație din spate Individul va muri Înregistrez numărul de mutații care au avut loc (va fi întotdeauna 8 în acest caz) Înregistrez orice mutație care are o frecvență în populație mai mare de 0,99 de indivizi. Acesta este numărul de fixări

Rulez acest cod de cel puțin 100 de ori. Mă aștept ca numărul mediu de fixări să fie egal cu numărul mediu de mutații din fiecare generație sau să convergă peste generații. Deoarece teoria neutră a evoluției de către Kimura prezice că pentru siturile neutre din genom, rata de substituție este egală cu rata mutației.


Acest lucru a devenit prea mare pentru un comentariu și cred că s-ar putea califica drept răspuns.

Poți să-mi dai o sursă pentru ceea ce citești? Conform acestui fapt, care se potrivește cu definițiile pe care le utilizați, K = (Nu) / N = u unde K este rata de substituție, u este rata de mutație în populație, Nu este rata de mutație individuală și N este dimensiunea populației .

M-am întors la lucrarea originală a lui Kimura în Nature în 1968 și el derivă că u (p), rata de substituție în populație, este (aproape exact) aceeași cu rata de mutație la un individ (p). Dar acest lucru rezultă dintr-o lucrare diferită (tip prolific) în care arată că probabilitatea de fixare pentru o anumită alelă este de aproximativ 1 / N, unde N este dimensiunea populației. Rata dvs. de mutație individuală în modelul dvs. nu este u, ci Nu.

Confundați rata mutației / substituției în populație, care este utilizată pentru ceasuri moleculare și alte lucruri, și rata mutației individuale. Dacă rata mutației moleculare este, de exemplu, 1/20 (așa cum presupuneți în modelul dvs.), atunci vă așteptați la o rată de fixare pe mutație de 1, în esență. Dacă rata de substituție va fi aceeași cu rata mutației moleculare, fiecare mutație trebuie să ajungă la fixare. Dacă același locus mută de două ori în aceeași generație (destul de posibil cu suficienți indivizi), atunci ambii trebuie cumva să ajungă la fixare, chiar dacă primul care ajunge la fixare îl va distruge în mod necesar pe celălalt.

Sper că acest lucru vă ajută să rezolvați de ce Nu va fi niciodată egal cu K și vă va ajuta să obțineți rezultate ale modelului dvs.

(bara laterală: Nu nu este rata mutației moleculare brute, ci rata mutației individuale, care este de fapt mai mare deoarece, în timp ce crește până la maturitate, gametul mediu uman este implicat în aproximativ 50 de divizii celulare, ceea ce duce la câteva alte mutații per persoană.)


Nu înțeleg pe deplin cum funcționează modelul dvs.

Nu înțeleg pe deplin modul în care îți modelezi populația. Este un $ oop $ (programare orientată pe obiecte) standard în care simulați fiecare individ? Sau este o simulare în care folosiți deja un model matematic? Puteți în cele din urmă să copiați și să inserați codul. probabil că nu este foarte lung, nu? Ce limbaj ați folosit? Ai un singur locus? Simulezi mai multe poziții în cadrul acestui locus? Permiteți mutații din spate?

Se pare că permiteți modificarea dimensiunii populației în timp și nu înțeleg prea bine cum parametrul dvs. $ K $ limitează însă rata de creștere. Ar trebui să scrieți câteva ecuații despre asta. În orice caz, modificarea dimensiunii populației va modifica probabilitatea fixării mutației.

Iată probabil câteva modificări importante de făcut pe codul dvs.

Este puțin greu să comentezi simulările tale fără să vezi unele ecuații sau să vezi codul, dar cred că vei avea nevoie de două modificări importante în codul tău.

  1. fixați dimensiunea populației $ N $.
  2. Pentru a simula reproducerea, alegeți aleatoriu (distribuție uniformă) indivizi (cu înlocuire) din populație. Indivizii aleși formează noua generație. De exemplu:

În R:

nou_pop = eșantion (vechi_pop, N, înlocuiește = T)

sau direct

pop = eșantion (pop, N, înlocuiește = T)

În Piton:

din import aleatoriu randrange new_pop = [] pentru i în xrange (N): new_pop.append (old_pop [randrange (N)])

Când obțineți rezultatele, nu uitați că vorbiți despre un proces stochastic și poate doriți să rulați suficientă simulare sau destui loci independenți pentru a putea compara cu predicțiile din modelele matematice.


Teoria neutră a genelor

Teoria neutră a devenit un instrument util pentru detectarea selecției naturale darwiniene la nivel molecular. Acest lucru se datorează faptului că predicțiile cantitative sunt disponibile în cadrul teoriei. Diverse abateri de la predicție sunt interpretate ca apariție a selecției. O abordare foarte populară este compararea modelelor de substituții sinonime și nesinonime în evoluția regiunilor de codificare a proteinelor. De exemplu, substituțiile rapide de aminoacizi în comparație cu modificările sinonime sunt considerate semne de selecție pozitivă. De asemenea, sunt încercate diferite încercări de a examina modelele de polimorfisme ale ADN-ului. În special, comparațiile polimorfismelor dintr-o populație și divergențele dintre speciile strâns legate au fost efectuate pentru a afla efectele selecției naturale.


Cuprins

Charles Darwin a comentat ideea mutației neutre în lucrarea sa, ipotezând că mutațiile care nu oferă un avantaj sau dezavantaj pot fluctua sau se pot fixa în afară de selecția naturală. „Variațiile nici utile, nici vătămătoare nu ar fi afectate de selecția naturală și ar fi lăsate fie un element fluctuant, așa cum probabil vedem la anumite specii polimorfe, fie ar deveni în cele din urmă fixe, datorită naturii organismului și a naturii condiții. " În timp ce lui Darwin i se atribuie larg introducerea ideii de selecție naturală care a fost în centrul studiilor sale, el a văzut, de asemenea, posibilitatea unor schimbări care nu au beneficiat sau nu au afectat un organism. [1]

Părerea lui Darwin despre schimbare fiind condusă în mare parte de trăsături care oferă avantaj a fost larg acceptată până în anii 1960. [2] În timp ce cerceta mutațiile care produc substituții de nucleotide în 1968, Motoo Kimura a constatat că rata de substituție a fost atât de mare încât, dacă fiecare mutație ar îmbunătăți condiția fizică, decalajul dintre genotipul cel mai potrivit și cel tipic ar fi plauzibil de mare. Cu toate acestea, Kimura a explicat această rată rapidă de mutație sugerând că majoritatea mutațiilor au fost neutre, adică au avut un efect redus sau deloc asupra aptitudinii organismului. Kimura a dezvoltat modele matematice ale comportamentului mutațiilor neutre supuse derivei genetice aleatorii la populațiile biologice. Această teorie a devenit cunoscută sub numele de teoria neutră a evoluției moleculare. [3]

Deoarece tehnologia a permis o analiză mai bună a datelor genomice, cercetările au continuat în acest domeniu. În timp ce selecția naturală poate încuraja adaptarea la un mediu în schimbare, mutația neutră poate împinge divergența speciilor datorită derivei genetice aproape aleatorii. [2]

Mutația neutră a devenit o parte a teoriei neutre a evoluției moleculare, propusă în anii 1960. Această teorie sugerează că mutațiile neutre sunt responsabile pentru o mare parte a modificărilor secvenței ADN la o specie. De exemplu, insulina bovină și umană, deși diferă în secvența de aminoacizi, sunt încă capabile să îndeplinească aceeași funcție. S-a văzut că substituțiile de aminoacizi între specii sunt neutre sau nu au impact asupra funcției proteinei. Mutația neutră și teoria neutră a evoluției moleculare nu sunt separate de selecția naturală, ci se adaugă gândurilor originale ale lui Darwin. Mutațiile pot oferi un avantaj, pot crea un dezavantaj sau nu pot face nicio diferență măsurabilă în ceea ce privește supraviețuirea unui organism. [4]

O serie de observații asociate cu mutația neutră au fost prezise în teoria neutră, inclusiv: aminoacizii cu proprietăți biochimice similare ar trebui înlocuiți mai des decât aminoacizii diferiți din punct de vedere biochimic, ar trebui observate substituții de baze sinonime mai des decât substituții nesinonime intronii ar trebui să evolueze la aceeași viteză ca mutații sinonime în exoni și pseudogene codificatoare ar trebui, de asemenea, să evolueze la o rată similară. Aceste predicții au fost confirmate odată cu introducerea de date genetice suplimentare de la introducerea teoriei. [2]

Mutație sinonimă a bazelor Edit

Când se introduce o nucleotidă incorectă în timpul replicării sau transcrierii unei regiuni de codificare, aceasta poate afecta eventuala translație a secvenței în aminoacizi. Deoarece sunt folosiți mai mulți codoni pentru aceiași aminoacizi, o schimbare într-o singură bază poate duce în continuare la traducerea aceluiași aminoacid. Acest fenomen este denumit degenerare și permite o varietate de combinații de codoni care duc la producerea aceluiași aminoacid. De exemplu, codurile TCT, TCC, TCA, TCG, AGT și AGC toate codurile pentru aminoacidul serină. Acest lucru poate fi explicat prin conceptul de oscilație. Francis Crick a propus această teorie pentru a explica de ce moleculele specifice de ARNt ar putea recunoaște mai mulți codoni. Zona ARNt care recunoaște codonul numit anticodon este capabilă să lege mai multe baze interschimbabile la capătul său 5 'datorită libertății sale spațiale. O a cincea bază numită inozină poate fi, de asemenea, substituită pe un ARNt și este capabilă să se lege cu A, U sau C. Această flexibilitate permite modificări ale bazelor din codoni care duc la traducerea aceluiași aminoacid. [5] Schimbarea unei baze într-un codon fără schimbarea aminoacidului tradus se numește mutație sinonimă. Deoarece aminoacidul tradus rămâne același, o mutație sinonimă a fost considerată în mod tradițional o mutație neutră. [6] Unele cercetări au sugerat că există o părtinire în selectarea substituției bazei în mutația sinonimă. Acest lucru s-ar putea datora presiunii selective pentru îmbunătățirea eficienței traducerii asociată cu cele mai disponibile ARNt sau pur și simplu cu prejudecăți mutaționale. [7] Dacă aceste mutații influențează rata de translație sau capacitatea unui organism de a produce proteine, ele pot influența de fapt capacitatea de sănătate a organismului afectat. [6]

Proprietăți biochimice ale aminoacizilor Nonpolar Polar De bază Acid Încheiere: opriți codonul
Cod genetic standard
Primul
baza
A 2-a bază A treia
baza
T C A G
T TTT (Phe / F) Fenilalanină TCT (Ser / S) Serine TAT (Tyr / Y) Tirozină TGT (Cys / C) Cisteină T
TTC TCC TAC TGC C
TTA (Leu / L) Leucina TCA TAA Stop (Ocru) [B] TGA Stop (Opal) [B] A
TTG [A] TCG ETICHETĂ Stop (Chihlimbar) [B] TGG (Trp / W) Triptofan G
C CTT CCT (Pro / P) Proline PISICĂ (His / H) Histidina CGT (Arg / R) Arginină T
CTC CCC CAC CGC C
CTA CCA CAA (Gln / Q) Glutamină CGA A
CTG [A] CCG CAG CGG G
A ATT (Ile / I) Isoleucină ACT (Thr / T) Treonină AAT (Asn / N) Asparagine AGT (Ser / S) Serine T
ATC ACC AAC AGC C
LA UN ACA AAA (Lys / K) Lizină AGA (Arg / R) Arginină A
ATG [A] (Met / M) Metionină ACG AAG AGG G
G GTT (Val / V) Valine GCT (Ala / A) Alanină GAT (Asp / D) Acid aspartic GGT (Gly / G) Glicină T
GTC GCC GAC GGC C
GTA GCA GAA (Glu / E) Acid glutamic GGA A
GTG GCG GAG GGG G
A Codonul ATG ambele coduri pentru metionină și servește ca situs de inițiere: primul ATG dintr-o regiune de codificare a mARN este locul în care începe traducerea în proteină. [8] Ceilalți codoni de start enumerați de GenBank sunt rare în eucariote și, în general, coduri pentru Met / fMet. [9] B ^ ^ ^ Baza istorică pentru desemnarea codonilor stop ca chihlimbar, ocru și opal este descrisă într-o autobiografie de Sydney Brenner [10] și într-un articol istoric de Bob Edgar. [11]

Substituție de aminoacizi neutri Edit

În timp ce înlocuirea unei baze într-o zonă necodificatoare a unui genom poate avea puține diferențe și poate fi considerată neutră, substituțiile de bază în sau în jurul genelor pot avea impact asupra organismului. Unele substituții de baze duc la mutații sinonime și la nicio diferență în aminoacidul tradus așa cum s-a menționat mai sus. Cu toate acestea, o substituție de bază poate schimba și codul genetic, astfel încât să se traducă un aminoacid diferit. Acest tip de substituție are de obicei un efect negativ asupra proteinei care se formează și va fi eliminat din populație prin selecția purificatoare. Cu toate acestea, dacă schimbarea are o influență pozitivă, mutația poate deveni din ce în ce mai frecventă într-o populație până când devine o piesă genetică fixă ​​a acelei populații. Organismele care se schimbă prin intermediul acestor două opțiuni cuprind viziunea clasică a selecției naturale. O a treia posibilitate este că substituția aminoacizilor face puțină sau deloc diferență pozitivă sau negativă în proteina afectată. [12] Proteinele demonstrează o anumită toleranță la modificările structurii aminoacizilor. Acest lucru este oarecum dependent de locul în care are loc substituirea în proteină. Dacă apare într-o zonă structurală importantă sau în situsul activ, o substituție de aminoacizi poate inactiva sau modifica substanțial funcționalitatea proteinei. Substituțiile din alte zone pot fi aproape neutre și pot deriva aleator în timp. [13]

Mutațiile neutre sunt măsurate în populația și genetică evolutivă, adesea prin examinarea variației populațiilor. Acestea au fost măsurate istoric prin electroforeză pe gel pentru a determina frecvența alozimelor. [14] Analizele statistice ale acestor date sunt utilizate pentru a compara variația cu valorile prezise în funcție de mărimea populației, ratele de mutație și dimensiunea efectivă a populației. Observațiile timpurii care au indicat o heterozigoitate mai mare decât se aștepta și variația generală în cadrul izoformelor proteice studiate, au condus argumente cu privire la rolul selecției în menținerea acestei variații față de existența variației prin efectele mutațiilor neutre care apar și distribuția lor aleatorie datorită derivării genetice. [15] [16] [17] Acumularea datelor pe baza polimorfismului observat a dus la formarea teoriei neutre a evoluției. [15] Conform teoriei neutre a evoluției, rata de fixare într-o populație a unei mutații neutre va fi direct legată de rata de formare a alelei neutre. [18]

În calculele inițiale ale lui Kimura, mutații cu | 2 Ns| & lt1 sau |s| ≤1 / (2N) sunt definite ca neutre. [15] [17] În această ecuație, N este dimensiunea efectivă a populației și este o măsurare cantitativă a dimensiunii ideale a populației, care presupune constante, cum ar fi rapoarte egale de sex și fără emigrație, migrație, mutație și nici selecție. [19] În mod conservator, se presupune adesea că dimensiunea efectivă a populației este de aproximativ o cincime din dimensiunea totală a populației. [20] s este coeficientul de selecție și este o valoare cuprinsă între 0 și 1. Este o măsurare a contribuției unui genotip la următoarea generație în care o valoare de 1 ar fi complet selectată și nu ar aduce nicio contribuție și 0 nu este deloc selectată împotriva. [21] Această definiție a mutației neutre a fost criticată datorită faptului că dimensiunile foarte mari ale populației eficiente pot face ca mutațiile cu coeficienți de selecție mici să pară non-neutri. În plus, mutațiile cu coeficienți de selecție mari pot apărea neutre la populații foarte mici. [17] Ipoteza testabilă a lui Kimura și a altora a arătat că polimorfismul în cadrul speciilor este aproximativ cel care ar fi de așteptat într-un model evolutiv neutru. [17] [22] [23]

Pentru multe abordări de biologie moleculară, spre deosebire de genetica matematică, mutațiile neutre sunt în general presupuse a fi acele mutații care nu cauzează niciun efect apreciabil asupra funcției genetice. Această simplificare elimină efectul diferențelor alelice minore în fitness și evită problemele atunci când o selecție are doar un efect minor. [17]

Dovezile convingătoare timpurii ale acestei definiții a mutației neutre s-au arătat prin ratele mutaționale mai mici în părțile importante funcționale ale genelor, cum ar fi citocromul c versus părțile mai puțin importante [24] și natura interschimbabilă funcțională a citocromului mamifer în studiile in vitro. [25] Pseudogenele nefuncționale oferă mai multe dovezi ale rolului mutațiilor neutre în evoluție. Ratele de mutație la pseudogenele globinei mamiferelor s-au dovedit a fi mult mai mari decât ratele genelor funcționale.[26] [27] Conform evoluției neodarwiniene, astfel de mutații ar trebui să existe rareori, deoarece aceste secvențe sunt nefuncționale și selecția pozitivă nu ar putea opera. [17]

Testul McDonald-Kreitman [28] a fost folosit pentru a studia selecția pe perioade lungi de timp evolutiv. Acesta este un test statistic care compară polimorfismul în siturile neutre și funcționale și estimează pe ce fracțiune de substituții s-a acționat prin selecție pozitivă. [29] Testul folosește adesea substituții sinonime în gene care codifică proteinele ca componentă neutră, totuși, mutațiile sinonime s-au dovedit a fi sub selecție purificatoare în multe cazuri. [30] [31]

Ceasurile moleculare pot fi utilizate pentru a estima cantitatea de timp de la divergența a două specii și pentru plasarea evenimentelor evolutive în timp. [32] Pauling și Zuckerkandl, au propus ideea ceasului molecular în 1962 pe baza observației că procesul de mutație aleatorie are loc la o rată constantă aproximativă. S-a arătat că proteinele individuale au rate liniare de modificări ale aminoacizilor pe parcursul timpului evolutiv. [33] În ciuda controverselor unor biologi care susțin că evoluția morfologică nu va continua într-un ritm constant, s-a demonstrat că multe modificări ale aminoacizilor se acumulează în mod constant. Kimura și Ohta au explicat aceste rate ca parte a cadrului teoriei neutre. Aceste mutații au fost considerate neutre, deoarece selecția pozitivă ar trebui să fie rară, iar mutațiile dăunătoare ar trebui eliminate rapid dintr-o populație. [34] Prin acest raționament, acumularea acestor mutații neutre ar trebui influențată doar de rata mutației. Prin urmare, rata de mutație neutră în organismele individuale ar trebui să se potrivească cu rata de evoluție moleculară a speciilor pe parcursul timpului evolutiv. Rata mutației neutre este afectată de cantitatea de situsuri neutre dintr-o secvență de proteine ​​sau ADN față de cantitatea de mutație din situsurile care sunt restricționate funcțional. Cuantificând aceste mutații neutre în proteine ​​și / sau ADN și comparându-le între specii sau alte grupuri de interes, se pot determina ratele de divergență. [32] [35]

Ceasurile moleculare au provocat controverse din cauza datelor pe care le derivă pentru evenimente precum radiații explozive observate după evenimente de dispariție precum explozia cambriană și radiațiile mamiferelor și păsărilor. Există diferențe de două ori în datele derivate din ceasurile moleculare și din înregistrările fosile. În timp ce unii paleontologi susțin că ceasurile moleculare sunt inexacte din punct de vedere sistemic, alții atribuie discrepanțele lipsei de date fosile robuste și a prejudecății în eșantionare. [36] Deși nu fără constanță și discrepanțe cu înregistrările fosile, datele din ceasurile moleculare au arătat cum evoluția este dominată de mecanismele unui model neutru și este mai puțin influențată de acțiunea selecției naturale. [32]


Rezultate

Modelele componente fundamentale ale simulatorului nostru au fiecare aplicații în instruire, simulare și inferență. Propunem că o singură teorie evoluțională coerentă, care să cuprindă și să unifice aceste trei aspecte, este o caracteristică extrem de dorită pentru un cadru de simulare. Nu poate exista o simulare fără parametri, care trebuie măsurați folosind un anumit model, parametrii sunt măsurați în mod tipic din date care trebuie să fi fost aliniate sau adnotate, folosind un anumit model, iar instrumentele de aliniere sau adnotare trebuie să fie comparate folosind simulări, folosind și un anumit model. Trăsătura comună a tuturor abordărilor descrise aici, pe care încercăm să o subliniem, este că toate aceste trei modele sunt aceleași: estimarea parametrilor, adnotarea și simularea pot fi realizate folosind același model probabilistic.

Acum descriem trei instrumente. GSIMULATOR generează aleatoriu aliniere ale ADN-ului în evoluție neutră. SIMGRAM generează aleatoriu alinieri ale caracteristicilor structurate generic sub selecție și SIMGENOME le combină pentru a genera aleatoriu aliniile regiunilor sintenice din genomi, folosind un model rezonabil detaliat al peisajului caracteristicilor genomice. Toate aceste trei programe, împreună cu exemple de fișiere de intrare / ieșire și utilizarea liniei de comandă, pot fi accesate prin intermediul paginii web [37].

GSIMULATOR: un simulator bazat pe traductor pentru ADN-ul în evoluție neutră

Instrumentul GSIMULATOR simulează evoluția neutră a ADN-ului pe un copac filogenetic. De-a lungul fiecărei ramuri, mutațiile de substituție și inserție-ștergere sunt modelate, folosind un traductor dependent de context. Teoria traductoarelor este descrisă în lucrările anterioare [21] și rezumată în Materiale și metode (a se vedea „Eșantionarea din traductoare lexicalizate”). În esență, un traductor este o mașină cu stare finită, asemănătoare cu un model Pair ascuns Markov („Pair-HMM”), care mută o secvență prin introducerea de substituții aleatorii și indels.

Traductorul GSIMULATOR este dependent de context, ceea ce înseamnă că substituția și indel „ratele” sunt dependente de trecut K simboluri absorbite și emise, unde K este un parametru care poate fi configurat. Scopul permiterii dependenței de context este de a modela fluctuațiile locale dependente de secvență în ratele de substituție și indel, cum ar fi dezaminarea CpG dependentă de metilare, expansiunea și contracția microsateliților și „micro-duplicări” sau „micro-inversiuni” corespunzătoare K nucleotide sau mai puține.

GSIMULATOR permite, de asemenea, modelarea flexibilă a distribuțiilor de lungime a spațiului, permițând N stări de inserare și ștergere degenerate multiple, unde N este un parametru configurabil. În absența dependenței de context (adică când K = 0), rezultă un amestec de distribuții geometrice pentru lungimile golurilor. Distribuția este mai complicată pentru traductoarele dependente de context (K ≥ 1), deoarece lungimea spațiului depinde de secvența inserată (sau ștearsă).

Ca și în cazul tuturor metodelor descrise în această lucrare, GSIMULATOR este un simulator antrenabil, ceea ce înseamnă că parametrii pot fi evaluați direct din datele de aliniere în perechi (și nu trebuie să fie „ghimitați” de către un utilizator). Biasurile compoziționale ale predictorilor genetici, aliniatorilor, identificatorilor de motive și altor instrumente de adnotare pot fi extrem de sensibile la statisticile evolutive care stau la baza acestora, deci această caracteristică este extrem de importantă pentru un simulator robust.

SIMGRAM: un instrument de simulare filogramatică

Instrumentul SIMGRAM generează eșantioane de eșantioane dintr-o gramatică fără context contextual filogenetică specificată de utilizator, sau „filogramatică”. Spre deosebire de traductoarele utilizate de GSIMULATOR, filogramaticele SIMGRAM pot modela caracteristicile genomului sub selecție fin structurată, inclusiv coordonarea aspectului relativ al acestor caracteristici (în mod similar cu modul în care o gramatică în limbajul uman specifică aspectul a diferitelor părți ale vorbirii). Exemple de caracteristici care pot fi astfel modelate includ gene de codificare a proteinelor [1, 18], gene de ARN necodificatoare [22, 38, 39], site-uri de legare a proteinelor [40], domenii de proteine ​​[41, 42] și structura secundară a proteinelor [22, 43]. Teoria completă a filogramaticelor a fost descrisă pe scară largă (vezi aplicațiile citate mai sus de noi și de alții) o scurtă introducere poate fi găsită în Materiale și metode (a se vedea „Eșantionarea din filograme”).

Formatul de filogramatică utilizat de SIMGRAM este același cu cel al programului XRATE, un instrument descris anterior care folosește algoritmul de maximizare a așteptărilor pentru a estima parametrii de rată și probabilitate ai oricărui filogramatică personalizată [22]. Prin urmare, XRATE poate fi utilizat pentru a estima parametrii de simulare direct din datele de antrenament, apoi SIMGRAM poate fi utilizat pentru a genera date sintetice cu proprietăți similare, dar fără omologie directă (cu excepția cazurilor în care filogramatica în sine codifică informații despre omologie [41, 42]) . Aceasta reprezintă o nouă aplicație a XRATE: programul SIMGRAM nu a fost descris anterior, iar utilizările anterioare ale XRATE au implicat utilizarea filogramei pentru a adnota secvența sau pentru a măsura ratele de substituție care sunt ele însele de interes direct. (Astfel, ca și în cazul tuturor metodelor de aici, modelul generativ care stă la baza SIMGRAM poate fi ușor instruit cu privire la date și utilizat pentru adnotare și deducere.)

O caracteristică atrăgătoare a filogramaticelor pentru simularea unei secvențe de genom în evoluție și bogată în caracteristici este că este extrem de ușor să combinați mai multe sub-modele într-un model cuprinzător. În plus, formatul XRATE permite mai multe caracteristici utile pentru simulare. O astfel de caracteristică este modelele parametrice, unde ratele și probabilitățile sunt constrânse să aibă o anumită formă funcțională care depinde de un set de parametri mai mic. Acest lucru este util pentru a construi modele care au proprietățile de simetrie dorite, cum ar fi procesele de substituție simetrică a firelor sau K A/K smodele de codoni. Formatul XRATE poate aproxima, de asemenea, modelele ratei de substituție dependente de context, folosind tehnica [44]. În cele din urmă, există un limbaj macro puternic care poate fi folosit pentru a descrie în mod compact gramaticile cu multe stări sau pentru a modela parametrizări specifice liniei.

Pachetul DART (teste ADN, aminoacizi și ARN), în care este distribuit SIMGRAM, include reproduceri ale mai multor filograme publicate anterior care pot fi simulate folosind SIMGRAM. Exemplele includ modele de structură secundară pentru proteine ​​[43] și ARN [38], precum și o serie de modele de substituție punctuală pentru proteine ​​și secvențe de acid nucleic. Formatul filogramei este descris integral online [45].

SIMGENOME: o filogramatică bogată în caracteristici pentru alinierea genomului

Până în prezent, am descris programul GSIMULATOR pentru mutația ADN neutră dependentă de context și programul SIMGRAM pentru mutația caracteristicilor structurate aflate sub selecție. Fiecare model are punctele forte și punctele sale slabe: GSIMULATOR modelează bogat ADN-ul neutru, dar nu caracteristici sub selecție, în timp ce SIMGRAM are modele mai bune pentru astfel de caracteristici, dar lipsește ratele dependente de context ale GSIMULATOR sau modelul indel sofisticat. Acum descriem un program care combină aceste abordări, utilizând un cadru modular care poate fi ușor extins pentru a încorpora viitoare simulatoare de caracteristici specializate.

Simulatorul combinat, SIMGENOME, începe prin generarea unei alinieri multiple dintr-un șablon de filogramatică care include o gamă bogată de caracteristici ale genomului. Funcțiile sunt descrise mai jos mai detaliat, împreună cu schițe despre modul în care șablonul poate fi extins de un utilizator expert pentru a încorpora noi caracteristici.

În alinierea generată, anumite coloane sunt marcate ca intergenice. Programul SIMGENOME apelează apoi în mod repetat la GSIMULATOR pentru a genera alinieri de ADN neutru corespunzătoare acestor regiuni intergenice și le îmbină în alinierea principală. Acest proces este extensibil: șablonul de filogramatică poate fi editat pentru a adăuga noi caracteristici sau pentru a modifica parametrii de bază ai modelului. În plus, alte simulatoare de caracteristici externe pot fi specificate în gramatica șablonului, iar aliniamentele lor de ieșire vor fi îmbinate în alinierea principală exact în același mod în care este ieșirea GSIMULATOR.

Caracteristicile modelate de șablonul filogramatică includ gene care codifică proteinele (cu o aproximare aproximativă la structura regiunii exon-intron-netradusă (UTR) care include distribuții ale lungimii exonului), gene ARN necodificatoare, elemente conservate (cum ar fi factorul de transcripție) situsuri de legare), pseudogene și transpozone ADN cu repetări inversate terminale. Funcțiile pot apărea pe șiruri înainte sau înapoi. Toate caracteristicile sunt adnotate în alinierea generată, astfel încât recuperarea lor în repere automate poate fi evaluată.

Sub-modelele care generează aceste caracteristici utilizează parametri ai ratei de substituție, care au fost estimate direct din următoarele seturi de date experimentale, utilizând programul XRATE (și ar putea fi re-estimate din seturi de date alternative).

Modelul genei care codifică proteinele

Modelul genei care codifică proteinele folosește o matrice de rate empirică, complet reversibilă și altfel necontrolată de 61 × 61 asupra codonilor, estimată în lucrările anterioare [46]. Ștergerile de păstrare a cadrelor sunt permise. Setul de instruire pentru acest model a fost baza de date PANDIT formată din alinieri la nivel de ADN ale familiilor de domenii proteice [19].

Modelul genei ARN necodificatoare

Modelul de genă ARN necodificator tratează golurile ca pe un al cincilea caracter și, prin urmare, folosește (4 + 1) × (4 + 1) nucleotidă monocatenară și (4 + 1) 2 × (4 + 1) 2 bază dublu catenară- matrici de rate de pereche care sunt complet reversibile și altfel necontrolate și au fost estimate separat. Setul de instruire pentru acest model a constat în alinieri furnizate cu programul CONSAN [47], care, la rândul lor, au fost derivate din baza de date europeană de ARNr subunitate mare [48]. Distribuția inițială a probabilității pe perechi de baze în regiuni dublu-catenare a fost, de asemenea, utilizată pentru a genera repetări inversate terminale în transpoziții de ADN simulate (deși aceste caracteristici evoluează ulterior sub un model neutru, astfel încât acestea să nu afișeze mutațiile compensatorii caracteristice non- codificând genele ARN (ncARN) aflate în selecție).

Modelul de substituție neutră simetrică a catenelor

Modelul de substituție neutră simetrică a catenelor care stă la baza modelelor de pseudogenă și transpozon tratează lacunele ca pe un al cincilea caracter și a fost antrenat la 1% aleatoriu din aliniamentele de 12 Drosophila genomi [20, 49]. Alinierile în sine au fost făcute folosind programul PECAN [50]. Modelul a fost constrâns să fie simetric de șir și reversibil utilizând funcționalitatea de parametrizare a XRATE. O versiune mai lentă și neîntreruptă a acestei matrice a ratei de substituție este, de asemenea, utilizată pentru modelarea caracteristicilor conservate.

Modelul traductorului

Modelul traductor utilizat de GSIMULATOR pentru a simula secvența intergenică care evoluează sub un model neutru dependent de context a fost instruit pe un set de aliniamente perechi trase dintr-un subset de douăsprezece specii Drosophila alinieri, care au fost realizate folosind programul PECAN [50]. Subsetul a fost extras din aproximativ 5% din datele inițiale de aliniere multiplă, cărora li s-a aplicat un prag de identitate procentual minim de 95%.

Frecvențele și distribuțiile de lungime

Frecvențele și distribuția lungimii caracteristicilor genomice au fost estimate din Drosophila literatura genomului [20, 51, 52] folosind Minos ca model pentru transpozonii ADN [53, 54].

Șablonul de bază filogramatică este scris utilizând formatul XRATE documentat public și poate fi ușor editat. Parametrii de nivel înalt, cum ar fi frecvențele cu care apar genele sau alte caracteristici, sunt declarați la nivelul superior al gramaticii și pot fi ușor modificați. Invităm utilizatorii să încerce să adauge sub-modele care reprezintă noi caracteristici relevante pentru benchmarking (sau să se consulte cu noi în acest proces). Noile sub-modele pot fi parametrizate direct din date folosind XRATE, iar parametrii copiați și lipiți în gramatica SIMGENOME, acest lucru se aplică și modelelor existente, care pot fi re-antrenate și / sau parametrizate (de exemplu, pentru a reflecta codoni diferiți de utilizare sau conținut GC). Extensibilitatea și mai mare este oferită de arhitectura modulară de plug-in care permite utilizarea programelor terță parte pentru a genera caracteristici care nu pot fi simulate în prezent de filograme sau traductoare SIMGENOME, cum ar fi matrice de tandem, duplicări pe distanță lungă și așa mai departe.

Următoarele sunt exemple de caracteristici care nu sunt incluse în prezent în modelul SIMGENOME, dar sunt posibile folosind SIMGRAM și ar putea fi încorporate prin modificarea fișierului gramatical al SIMGENOME.

Frecvențe codon

Setul de date PANDIT care a fost utilizat pentru a estima modelul de codoni SIMGENOME se întinde pe un spectru larg de prejudecăți de utilizare a compoziției și codonilor. În ceea ce privește tiparele generale de conservare și ratele de mutație suprimate, considerăm că este un model general acceptabil de substituție a codonilor. De exemplu, dacă se compară un instrument de identificare a motivelor, regiunile de codificare generate de SIMGENOME vor avea ca rezultat mai mulți falsi pozitivi decât regiunile care nu codifică (deoarece nivelul de conservare este mai mare) și acest lucru poate fi adecvat în scopul acelei etape de referință . Cu toate acestea, în alte scopuri, s-ar putea dori un model parametric mai bogat, care să țină seama de efectele specifice genomului, cum ar fi părtinirea compozițională, părtinirea codonilor, raporturile de tranziție-transversie sau dezaminarea indusă de metilarea CpG. Utilizarea unui astfel de model parametric este un caz simplu de schimbare a matricei de rată relevante în fișierul de gramatică SIMGENOME. Modelul poate fi potrivit cu datele folosind XRATE în mod normal. Pregătim un manuscris care descrie o comparație directă între XRATE și PAML în aceste scopuri, inclusiv codul Perl pentru generarea unor astfel de modele parametrice mai bogate (Heger A, Ponting C și Holmes IH, în pregătire).

Parametrizări specifice liniei

Limbajul macro SIMGRAM permite parametri diferiți pe diferite ramuri ale arborelui. Nu am folosit această caracteristică în gramatica SIMGENOME, deoarece utilizarea sa este oarecum dependentă de clada filogenetică investigată: s-ar putea (de exemplu) să doriți să utilizați parametri diferiți pe fiecare ramură, pe o singură ramură internă sau într-o cladă specifică. Cu referire la documentația XRATE, este foarte posibil să se proiecteze gramaticile care folosesc această caracteristică, astfel încât modelul să poată (de exemplu) să utilizeze frecvențe de codoni diferite sau părtiniri compoziționale în diferite părți ale arborelui.

Mutații de pierdere a funcției

Formatul SIMGRAM permite, de asemenea, evoluția specifică liniei de caracteristici întregi, în modul programului DLESS de către Pollard și colab. [55]. Pe plan intern, am dezvoltat filograme care modelează mutații de pierdere a funcției în genele care codifică proteinele, pentru investigații ale evoluției pseudogenei. Nu am inclus mutații de pierdere a funcției în această primă versiune de SIMGENOME, dar ar fi destul de fezabil să o extindem fără a fi nevoie să scriem un cod nou.

Situri de îmbinare, metionine inițiale, UTR și alte aspecte ale structurii genei care codifică proteinele

Gramatica SIMGENOME include în prezent o aparență brută a structurii exonului și a intronului, pentru a reproduce fluctuațiile compoziționale largi asociate cu genele care codifică proteinele. Aceasta include în prezent o secvență cu evoluție lentă la limitele exon-intron, ca o machetă de conservare a sitului de îmbinare. Eliberarea actuală nu modelează, totuși, genele care codifică proteinele la un nivel suficient de detaliu pentru a fi utilizate ca exemple pozitive pentru un predictor al genelor care codifică proteinele. S-ar putea modifica cu ușurință simulatorul pentru a face acest lucru, modelând site-uri de îmbinare ca GT-AG și poate chiar perechi donator-acceptor AT-AC și incluzând alte caracteristici precum semnale poli-A, ATG-uri inițiale, cutii TATA și UTR-uri. Aceste tipuri de caracteristici sunt ușor de adăugat într-un cadru filogramatic.

Corelații de ordin superior între codoni

Distribuția SIMGRAM include exemple de filo-HMM care demonstrează corelații de ordin superior între aminoacizi în aliniamentele proteice. De exemplu, există o replică a filo-HMM cu 3 stări Thorne-Goldman-Jones pentru modelarea structurii secundare [43].Având în vedere o parametrizare adecvată pentru maparea de la matricele de viteză de substituție a aminoacizilor la matricile de codoni, ar fi simplu să se utilizeze ceva de genul acesta pentru a modela corelații la nivel superior între codoni în regiunile de codare-ADN ale SIMGENOME. S-ar putea întreba în mod rezonabil de ce nu am inclus astfel de dependențe de ordin superior în prima versiune a SIMGENOME, când am inclus dependențe de context de ordin superior în regiunile intergenice. Răspunsul este că caracteristicile intergenice (cum ar fi microsateliții sau regiunile conservate) pot include părtiniri compoziționale puternice care se extind pe zeci de baze și, astfel, în general contribuie mai mult la fluctuațiile densității informațiilor decât corelațiile dintre codoni, care sunt de obicei slabe, locale și adesea detectabile doar la nivelul aminoacizilor codificați [43].

Distribuții detaliate de lungime

Majoritatea distribuțiilor de lungime ale caracteristicilor modelate în SIMGENOME sunt geometrice (cel mai simplu tip de distribuție pe care o puteți modela cu un HMM). Înlănțuind stări duplicate, este posibil să se genereze distribuții de lungime mai complexe (realiste). De exemplu, fișierul gramatic SIMGENOME include funcții pentru generarea unei distribuții de vârf „binom negativ” (cunoscută și sub numele de distribuție „Pascal”) pentru lungimi de exoni, deci există exemple care pot fi utilizate pentru a face acest lucru. Cu toate acestea, s-a găsit empiric o distribuție geometrică simplă care se potrivește mai bine cu datele lungimii exonului Drosophila.

Evaluare: compararea unui predictor de ARN necodificat

Evaluarea unui instrument de simulare este o problemă ușor diferită de evaluarea unui instrument de adnotare, cum ar fi un predictor de gene. Când se compară un predictor, este de obicei interesat să minimizeze numărul de pozitive false pe care predictorul le găsește într-un set de date nul la un prag dat al limitei scorului de predicție. Se pot elimina toate pozitivele false prin stabilirea pragului de scor arbitrar ridicat, dar cu prețul lipsei, de asemenea, a tuturor genelor reale din setul de date. Pentru a obține o imagine reală a performanței predictorului, trebuie deci să luăm în considerare modul în care numărul de fals pozitivi variază în funcție de acest prag sau (mai semnificativ) în funcție de sensibilitatea predictorului, adică numărul de gene reale că detectează corect la un anumit prag de scor.

O practică obișnuită este de a utiliza un set de gene reale pentru a evalua sensibilitatea predictorului, dar de a utiliza date nule simulate pentru a evalua rata fals pozitivă. Motivul utilizării datelor nule simulate, în loc de o secvență ADN reală care nu conține nicio genă, este lipsa de adnotări negative: este dificil de demonstrat experimental că o anumită secvență de ADN definitiv nu conține nicio genă. Acest lucru este valabil mai ales pentru genele greu de identificat, cum ar fi genele ARN sau cadrele scurte de citire deschise.

Prin urmare, un simulator bun este unul care reproduce îndeaproape statisticile ADN-ului real și generează un număr similar de fals pozitivi cu ADN-ul real. În practică, al doilea criteriu se bazează pe alegerea statisticilor pentru primul criteriu. În general, dacă modelul statistic este prea simplist (de exemplu, omiterea regiunilor cu complexitate redusă), atunci complexitatea ADN-ului simulat va fi mai mult decât ar trebui fi, ducând la mai puține pozitive false. Pentru a lua un exemplu specific: un simulator de ADN care emite aleatoriu o secvență de simboluri independente, distribuite identic, se va potrivi cu compoziția la nivel de nucleotidă a ADN-ului real, dar nu va reproduce fluctuațiile pe termen scurt ale conținutului informațional care poate fi găsit în ADN-ul real , și așa va genera mai puțini falsi pozitivi într-un etalon de identificare a motivelor decât un simulator care include fenomene de fluctuație a complexității pe termen scurt (cum ar fi microsateliții).

Pe această bază, susținem că o măsură a unui simulator probabilistic bun este că modelul său (după ce se potrivește cu ADN-ul real) ar trebui să maximizeze numărul de fals pozitivi pentru un predictor dat la o sensibilitate dată. O altă modalitate de a pune acest lucru este că simulatorul ar trebui să ofere cel mai strict criteriu de referință posibil pentru predictorul genei, prin minimizarea ariei de sub curba caracteristicii de funcționare a receptorului (ROC) peste intervalul de interes (a se vedea mai jos).

Pentru a evalua motoarele de simulare descrise aici, le-am folosit pentru a estima rata fals pozitivă (FPR) pentru un ecran de calcul al întregului genom pentru genele ARN structurale conservate, efectuate utilizând XRATE [22]. Am ales predicția genei ARN ca caz de testare deoarece are un FPR extrem de ridicat, a cărui întindere reală este încă necunoscută [39, 56, 57] și pentru că FPR estimat pentru acest ecran este extrem de sensibil la proprietățile subiacente ale motorului de simulare, făcându-l un bun motivator al realismului sporit [12, 58].

Mai exact, FPR a fost estimat prin glisarea unei ferestre peste aliniamentele simulate și executarea XRATE pe fiecare fereastră folosind o gramatică de predicție a genei ARN. Această gramatică modelează modelele distincte de substituție a nucleotidelor în genele ARN, inclusiv covariația nucleotidelor asociate cu baze și este strâns legată de programul EVOFOLD pentru predicția comparativă a genelor ARN [39] și programul PFOLD pentru plierea ARN comparativ [59]. Gramatica în sine, împreună cu instrucțiuni detaliate pentru reproducerea ecranului, pot fi găsite online [37]. Justificarea completă care stă la baza dezvoltării gramaticii și evaluarea critică a acesteia (ca genefinder) și comparația cu gramaticile conexe vor fi descrise în altă parte (Bradley RK, Uzilov AV, Skinner M, Bendaña YR, Barquist L și Holmes I, prezentate) .

Graficele prezintă curbe ROC în care FPR este reprezentat grafic în funcție de sensibilitatea ecranului (măsurată utilizând adnotări de ncRNA cunoscute în Drosophila melanogaster), ca funcție parametrică a pragului de scor pentru ecran. Deoarece ncRNA reprezintă rezultate pozitive pentru acest ecran și avem deja un set de ncRNA curate cunoscute pentru D. melanogaster, am omis submodelul ncRNA din gramatica SIMGENOME pentru aceste teste.

Rezultatele sunt prezentate în Figura 1. Concluzia generală este că realismul crescut face ca un FPR mai mare. În cazul GSIMULATOR, fie în creștere N sau K crește radical FPR în cazul SIMGENOME, includerea caracteristicilor genomice conservate cu rate evolutive mai lente crește semnificativ FPR în raport atât cu un model de substituție punct pur, cât și cu cel mai realist model GSIMULATOR.

Curbe caracteristice de funcționare a receptorului (ROC) pentru doi algoritmi de predicție ARN necodificatori, ClosingBp (Bradley RK, Uzilov AV, Skinner M, Bendaña YR, Barquist L și Holmes I, prezentat) și EVOFOLD [39] (implementat utilizând XRATE), utilizând Modelele GSIMULATOR și SIMGENOME pentru a estima rata de descoperire fals pozitivă. Aceste curbe ilustrează principiul general că, cu cât un model de simulare este mai realist, cu atât este mai mare rata estimată de fals pozitiv (FPR). Această tendință este independentă de algoritmul de predicție genică utilizat. Panourile superioare arată rezultate pentru GSIMULATOR: se vede că distribuții mai complexe ale lungimii indel (N) și, în special, dependența de context (K) ambele cresc FPR. Panourile inferioare prezintă rezultate pentru SIMGENOME și modelele componente, în care FPR este crescut prin includerea lacunelor (care amplifică fluctuațiile conținutului informațional, datorită faptului că sunt tratate de obicei ca „informații lipsă”) și caracteristici genomice (dintre care unele evoluează mai lent rata decât secvența neutră). Motivul pentru care sensibilitatea asimptotică este mai mică de 1,0 este că standardul nostru de referință a folosit o abordare cu fereastră glisantă, prezicând cel mult un ARN necodificator (ncRNA) în fiecare fereastră. Setul nostru de ncRNA reale a fost preluat din multi-genom Drosophila alinieri produse de programul PECAN [50] în fiecare caz, pentru a asigura o comparație corectă, am luat o fereastră a aliniamentului PECAN care înconjoară ncRNA adnotat, cu dimensiunea acestei ferestre care se potrivește cu dimensiunea ferestrei glisante care a fost utilizată pe datele nule simulate. Unele dintre ncRNA-urile pozitive din aceste ferestre aliniate PECAN au un scor atât de slab în cadrul modelului de predicție genică - de exemplu, din cauza inexactităților din alinierea PECAN a acelei ferestre - încât ncRNA-ul prezis este plasat în mod constant în locația greșită din fereastră. Prin urmare, acești ncARN reali nu sunt detectați niciodată, oricât de scăzut ar fi pragul de notare, stabilind o limită superioară pentru sensibilitatea realizabilă.

De asemenea, am comparat metodele noastre de simulare, GSIMULATOR și SIMGENOME, cu DAWG [10], un program larg citat pentru simularea substituției neutre și a evenimentelor indel. Am ales DAWG deoarece exemplifică cel mai bine obiectivele pe care le-am identificat aici: se bazează în mod clar pe un model evolutiv subiacent și oferă instrumente pentru estimarea parametrilor modelului indel direct din datele secvenței. Se pare că este principalul simulator de uz general în momentul scrierii. Alte simulatoare (cum ar fi PSPE) sunt mai bogate, dar nu oferă funcționalitatea de estimare a parametrilor pe care o face DAWG.

Parametrii pentru DAWG au fost după cum urmează. Am folosit același model de substituție reversibilă în timp general (REV) pe care l-am estimat din aliniamentele PECAN Drosophila genomi. Modelul de indel „geometric” al DAWG (adică lungimi de indel distribuite geometric) a fost parametrizat folosind scriptul furnizat împreună cu programul. Deși modelul „putere-lege” pentru lungimile indel a dat o potrivire mai bună, a produs alinieri care erau în mare parte lacune. DAWG permite rate eterogene și site-uri invariante folosind Γ + Eu model pentru eterogenitatea ratei [60], oferind și câțiva exemple de parametri pentru acest model (γ = 1, ι = 0,1), pe care le-am folosit pentru aceste simulări.

Figura 2 compară DAWG cu cel mai realist GSIMULATOR și cele mai bogate modele SIMGENOME. Figura arată că dependența de context modelată de GSIMULATOR și caracteristicile genomice modelate de SIMGENOME au ca rezultat estimări fals pozitive mult mai stricte decât cele produse de DAWG.

Curbe ROC pentru doi predictori de ARN necodificatori, ClosingBp (Bradley RK, Uzilov AV, Skinner M, Bendaña YR, Barquist L și Holmes I, prezentat) și EVOFOLD [39] (implementat utilizând XRATE), comparând DAWG [10] cu cele mai bogate modele GSIMULATOR și SIMGENOME. Cele trei curbe pentru fiecare predictor de gene ilustrează clar că bogăția crescută a modelului (DAWG → GSIMULATOR → SIMGENOME) produce FPR estimat mai mare. Vedeți legenda din Figura 1 pentru o explicație a motivului pentru care sensibilitatea asimptotică este mai mică de 1,0.

Deoarece exemplul DAWG parametrii γ și ι (care determină eterogenitatea ratei și densitatea sitului conservat) au fost reglați manual de către autorul programului pentru alinieri om-cimpanzei, acestea pot fi o subestimare pentru Drosophila (unde elementele conservate sunt mai strâns distanțate decât la primate, datorită genomului mai mic și a unei rate de ștergere mai mari a ADN-ului nefuncțional). Cu toate acestea, nu există nicio metodă automatizată pentru setarea acestor parametri în DAWG și nici pachetul sau lucrarea DAWG nu oferă îndrumări explicite cu privire la modul de legare a acestor parametri la statistici directe privind densitatea caracteristicilor genomice. (În schimb, fișierul de gramatică SIMGENOME include comentarii care descriu derivarea parametrilor de distribuție a caracteristicilor publicate Drosophila adnotări.) În plus, programul GSIMULATOR include o procedură de instruire complet automatizată, nu are model pentru site-urile conservate sau heterogenitatea ratei (altele decât substituția dependentă de context) și generează totuși o rată de predicție falsă mai mare decât DAWG, chiar și atunci când aceste caracteristici sunt activate în DAWG (Figura 2). Deoarece SIMGENOME, la rândul său, generează un FPR mai mare decât GSIMULATOR (Figura 1), considerăm tranzitoriu că SIMGENOME este, de asemenea, un simulator mai realist decât DAWG. În mod corespunzător, observăm că acest argument motivează puternic un instrument automat (sau chiar un simplu euristic) pentru estimarea parametrilor γ și ι ai modelului de rată eterogenă al DAWG. Există o analiză considerabilă a acestor modele [60], deci acesta poate fi un obiectiv rezonabil al viitorului studiu.


Introducere

Una dintre cele mai mari provocări din domeniul geneticii medicale și populației moderne este determinarea consecințelor fenotipice și de fitness ale unei anumite mutații. Studiile de asociere la nivel de genom (GWAS) au implicat sute de loci din genom pentru controlul multor trăsături [1]. Cu toate acestea, găsirea variantelor cauzale la acești loci a rămas provocatoare datorită corelațiilor statistice dintre markeri (dezechilibru de legătură) și a faptului că majoritatea loviturilor GWAS cad în regiuni necodificate ale genomului cu funcție puțin evidentă [2]. Cunoașterea anumitor variante cauzale este un obiectiv important, deoarece va îmbunătăți predicția riscului și va permite o înțelegere mai detaliată a mecanismului biologic din spatele modului în care varianta influențează trăsătura. În genetica populației, există un interes extraordinar în a înțelege cât de mult din genom se află sub selecție și tipurile de mutații care stau la baza multor variații fenotipice și adaptare la diferite specii. Mai mult, studiile au urmărit să cuantifice cu precizie cantitatea de variație dăunătoare care separă populațiile pentru a evalua rolul istoriei populației la influențarea variației dăunătoare și pentru a determina dacă dimensiunea redusă a populației ar putea duce la o acumulare de variante dăunătoare, cauzând potențial o topire mutațională și dispariția [3-6].

O modalitate populară de a evalua care mutații dintr-un genom pot fi funcționale din punct de vedere biologic și pot afecta condiția fizică este de a examina măsura în care nucleotidele sunt conservate între taxonii la distanță evolutivă. Se consideră că site-urile care prezintă un deficit de substituții în mai multe linii sunt importante din punct de vedere funcțional și sunt supuse selecției purificatoare. Se consideră că site-urile care prezintă un număr mai mare de substituții evoluează la o rată neutră și ar fi mai puțin probabil să fie funcționale sau sub selecție purificatoare. Au fost dezvoltate o serie de abordări statistice pentru a găsi aceste situri în genom, care prezintă conservare la specii disparate [7-14]. În plus, acest concept a fost utilizat în mai multe instrumente de adnotare, cum ar fi scorurile SIFT, PolyPhen și CADD pentru a prezice care mutații sunt susceptibile de a fi dăunătoare [15-18].

O abordare genomică comparativă specială care a primit o utilizare pe scară largă este scorul de profilare a ratei evolutive genomice (GERP) [19,20]. Scorul GERP este definit ca reducerea numărului de substituții în alinierea secvenței multi-specii comparativ cu așteptarea neutră. De exemplu, un scor GERP de 4 ar însemna că există 4 substituții mai puține la un anumit sit decât ceea ce se așteaptă pe baza ratei neutre de evoluție a filogeniei. Ca atare, scorul GERP este o măsură a conservării secvenței la mai multe specii. Cu toate acestea, scorurile GERP au fost utilizate în mod obișnuit în studiile genomice evolutive ca măsură a puterii selecției care acționează asupra mutațiilor derivate care separă într-o specie. În aceste aplicații, se presupune că mutațiile care apar în siturile care sunt foarte conservate în multe specii sunt dăunătoare și, astfel, contribuie la încărcarea genetică a unei specii. Cantitativ, pentru fiecare sit de segregare din cadrul unei specii, scorul GERP este atribuit mutației derivate care segregează la acel sit. De exemplu, Schubert și colab. [3] au studiat tiparele de mutații dăunătoare la caii sălbatici și domestici. Au calculat o sarcină de scor GERP pentru fiecare cal, care a fost scorul mediu GERP pentru toate variantele derivate din acel individ. Ei au descoperit o creștere a sarcinii scorului GERP la calul domesticit, argumentând că domesticirea a dus la o creștere a variației dăunătoare. Henn și colab. [5] a folosit GERP pentru a evalua impactul asupra condițiilor fizice al mutațiilor care modifică aminoacizii la om. Ei au definit mutații cu scoruri GERP 4-6 pentru a avea efecte nocive „mari”, corespunzătoare unui coeficient de selecție de 10-3, și au observat o creștere a numărului acestor alele derivate dăunătoare la populațiile non-africane. Ei au raportat, de asemenea, că scorurile GERP rezumate pe toate site-urile dintr-un individ au fost mai mari într-un genom mediu maya nativ american comparativ cu un genom mediu african sub-saharian. Marsden și colab. [4] au folosit GERP pentru a identifica mutațiile care modifică aminoacizii dăunători la câini și lupi și au constatat o creștere a mutațiilor dăunătoare (GERP & gt4) la câini și că câinii au un scor GERP mai mare în toate variantele de schimbare a aminoacizilor în comparație cu lupii. În cele din urmă, Valk și colab. [6] au constatat că, într-o gamă largă de mamifere, speciile cu dimensiuni ale populației din punct de vedere istoric reduse și diversitatea genetică scăzută au un scor GERP mediu mai scăzut alelei derivate decât speciile cu dimensiuni mari ale populației, sugerând că purjarea alelelor dăunătoare reduce sarcina genetică în populațiilor pe termen lung.

În timp ce scorurile GERP au fost utilizate pe scară largă în domeniul geneticii medicale și al populației, rămân unele provocări. În primul rând, studiile descrise mai sus împart scorurile GERP într-un mod grosier pentru a reflecta coeficientul de selecție dăunător subiacent. Se presupune că acele mutații cu un scor GERP mai mare au un coeficient de selecție mai dăunător. Cu toate acestea, precizia atribuirii scorurilor GERP la anumite efecte de fitness rămâne neclară. Scorurile GERP nu pot furniza dovezi cantitative ale puterii selecției, deoarece orice mutații dăunătoare care au un coeficient de selecție scalat de Nes & lt -2 nu se vor acumula ca substituții [21-24]. Sub această valoare, nu se vor acumula mutații slab sau puternic dăunătoare ca substituții și, prin urmare, este posibil să nu fie posibilă distincția între ele folosind date genomice comparative [23]. În al doilea rând, majoritatea metodelor de detectare a conservării presupun presiuni constante de selecție în toate ramurile unei filogenii [8]. Orice tip de selecție specifică descendenței sau fluctuației secvenței funcționale (adică o secvență are un rol regulator specific într-o descendență, dar nu într-o altă descendență), ar putea fi ratată de aceste abordări genomice comparative. Dovezi recente au sugerat o cantitate echitabilă de cifră de afaceri a secvenței funcționale în regiunile necodificate ale genomului uman [25, dar vezi 26]. În cele din urmă, s-a arătat că puterea metodelor de genomică comparativă pentru a detecta secvențele sub selecție poate fi maximizată prin selectarea subseturilor optime dintr-un set mai mare de specii [27]. Cu toate acestea, subsetul optim de specii pentru a maximiza performanța în diferite scenarii de selecție și de rotație rămâne neclar. Acest lucru este preconizat în special în lumina proiectelor recente care vizează creșterea numărului de genomi secvențați între specii [28,29].

În cele din urmă, amploarea genomului uman sub selecție purificatoare a rămas sub o dezbatere puternică. Studiile genomice comparative timpurii au sugerat că cel mult 15% din genom se afla sub selecție [9,20,30-33].Cu toate acestea, studiile biochimice efectuate de ENCODE au sugerat că până la 80% din genom prezintă activitate în cel puțin un test biochimic [34]. Este posibil să se concilieze aceste estimări observând că acestea măsoară diferite procese - testele funcționale evaluează dacă nucleotida are activitate biochimică, dar această activitate nu poate fi neapărat legată de fitness [35,36]. Ca atare, mutațiile la siturile biochimic active pot să nu aibă un impact evolutiv și, prin urmare, ar putea părea neutre în abordările genomice comparative. Mai mult, au existat dovezi din studiile genomice comparative ale cifrei de afaceri a secvenței supuse selecției de purificare [25,30,33,37,38]. Acest lucru ar putea avea loc în mai multe moduri. În primul rând, secvențele pot avea o funcție biologică la unele specii și nu la altele datorită modificărilor arhitecturii de reglementare între specii [39]. În al doilea rând, chiar dacă regiunea de reglare își păstrează funcția biologică pe perioade evolutive îndelungate, coeficienții de selecție ai mutațiilor la anumite situri s-ar putea modifica în timp datorită efectelor epistatice cu alte mutații [40]. Rands și colab. sugerează că istoria evolutivă a genomului uman a fost extrem de dinamică, doar 25% dintre elementele aflate sub selecție purificatoare la om au menținut constrângerea la șoarece [25,30]. Alte studii au sugerat că cifra de afaceri evolutivă mai recentă a avut un impact redus asupra conținutului funcțional al genomului [26]. Astfel, rămâne o întrebare deschisă cu privire la cât de mult din genom se află sub selecție purificatoare și cantitatea de cifră de afaceri a secvenței funcționale care are loc.

Aici efectuăm simulări realiste sub modele genetice ale populației de selecție purificatoare pentru a evalua performanța scorurilor GERP în diferite scenarii. Mai întâi evaluăm dacă scorurile GERP pot furniza estimări fiabile ale coeficienților de selecție la mutațiile de codificare individuale. Apoi evaluăm măsura în care cifra de afaceri a secvenței afectează capacitatea GERP de a identifica secvențele selectate la site-uri necodificate. În cele din urmă, estimăm că cel puțin 4,51% din porțiunea necodificatoare a genomului uman se află sub selecție purificatoare și că mutațiile de la majoritatea acestor site-uri necodificate nu au fost sub selecție de-a lungul întregii evoluții a mamiferelor. Rezultatele noastre indică câteva limitări importante la utilizarea abordărilor genomice comparative pentru determinarea efectelor de fitness ale mutațiilor individuale și se adaugă la literatura în creștere, care susține utilizarea datelor polimorfismului pentru evaluarea cantităților actuale de selecție în cadrul speciilor [23].


Metode

Filtrarea secvenței SARS-CoV-2 GISAID

După eliminarea oricărei secvențe cu o gazdă neumană (de exemplu, liliac, pangolin), pentru a reduce impactul erorilor de secvențiere asupra analizei de selecție, datele de la GISAID au fost filtrate prin excluderea tuturor secvențelor care îndeplinesc oricare dintre următoarele criterii: orice secvență de lungime mai mică mai mult de 29.000 nucleotide orice secvențe cu nucleotide ambigue care depășesc 0,5% din genom orice secvențe cu mai mult de 1% divergență de cea mai lungă secvență eșantionată (Wuhan-Hu-1) și orice secvență cu codoni stop. Pentru a păstra alinieri de codoni și a păstra secvențe de aminoacizi corecte pentru analize de selecție, unice (secvențe identice sunt prăbușite pentru alinierea pentru a reduce cheltuielile de calcul „N” este tratat ca potrivind cu orice caracter rezolvat în acest proces) secvențele nucleotidice în cadru au fost traduse în amino- acizi și aliniați cu MAFFT [69]. Alinierea aminoacizilor este mapată înapoi la secvențele de nucleotide constitutive pentru a produce o aliniere la nivel de codoni. Numai haplotipurile unice sunt păstrate pentru analize filogenetice comparative, deoarece includerea copiilor identice nu este informativă pentru aceste tipuri de inferență.

SARS-CoV-2 selecție pozitivă

Am folosit inițial o metodologie care a încorporat filogenia completă și a găsit semnături înșelătoare, cum ar fi erori de secvențiere și recombinare bazată pe laborator pe ramurile terminale, care au confundat software-ul (textul S2). Aceste erori se reflectă în raportul dN / dS ridicat observat pe ramurile terminale (Fig 1). Ulterior am folosit metodele FEL [30] și MEME [36] pentru a deduce selecția pozitivă diversificând negativ și episodic, respectiv, folosind o implementare care a considerat doar ramuri interne pentru a deduce selecția pozitivă. MEME utilizează o metodologie de maximă probabilitate și efectuează un test al raportului de probabilitate pentru selecția pozitivă pe fiecare site, comparând modurile care permit sau nu permit selecția pozitivă de diversificare la un subset de ramuri (dN / dS & gt 1). FEL efectuează un test care presupune o presiune selectivă uniformă pe toate ramurile (selecție omniprezentă). Aceste analize de selecție au fost efectuate în pachetul software HyPhy v.2.5.14.

Sarbecovirusuri alinierea și recombinarea

Pentru a evita efectele confundante ale recombinării, am analizat fiecare ORF separat și am împărțit ORF1ab și Spike ORF în regiuni supuse non-recombinante, pe baza celor 7 puncte majore de recombinare prezentate la Boni și colegii [6]. Aceasta produce 5 regiuni nerecombinante pentru Orf1ab (regiunile A la E) și 5 regiuni pentru Spike (regiunile A la D și bucla variabilă - regiunea VL). Secvențele de proteine ​​din regiunile nerecombinante SARS-CoV-2, SARS-CoV-1 și 67 de virusuri strâns legate cu gazde neumane (lilieci și pangolini, tabelul S4), identificate pe baza similarității secvenței și recuperate din bazele de date online NCBI Genbank și GISAID , au fost aliniate folosind versiunea MAFFT 7 (L-INS-i) [69]. Corecțiile manuale ulterioare au fost făcute pe aliniamentele de proteine ​​și PAL2NAL (http://www.bork.embl.de/pal2nal) a fost utilizat pentru a le converti în alinieri de codoni. Filogeniile pentru fiecare aliniere a codonilor au fost deduse folosind RAxML cu un model de substituție a nucleotidelor GTR + Γ [70].

Sarbecovirus analiza selectiei

Am folosit o serie de metode de detectare a selecției pentru a examina dacă descendența care duce la SARS-CoV-2 a cunoscut episoade de selecție pozitivă diversificată. Fiecare regiune nerecombinantă a fost examinată separat. Am separat filogenia fiecărei regiuni într-o linie nCoV și non-nCoV / SARS-CoV-1. Clada nCoV include SARS-CoV-2 și virusurile care formează o monofilie cu aceasta, excluzând SARA-CoV-1 care conține clada soră. Acestea sunt virusurile infectante cu lilieci CoVZC45, CoVZXC21, RmYN02 și RaTG13 și virusurile infectante cu pangolin Pangolin-CoV și grupul P2V, P5L, P1E, P5E, P4L (tabelul S4). Notă, unele regiuni recombinante ale CoVZC45, CoVZXC21 și RmYN02 nu aparțin cladei nCoV și acestea au fost excluse din orice analiză a acestor regiuni.

Am testat pentru evidența selecției episodice de diversificare pe ramurile interne ale cladei nCoV folosind BUSTED [S], luând în considerare variația sinonimă a ratei, așa cum este descris în Wisotsky și colegii [71]. Am dezvoltat o extensie la BUSTED [S], care a inclus un HMM cu 3 categorii de rate pentru a descrie variația sinonimelor specifice site-ului [43]. Acest HMM permite încorporarea explicită a autocorelării în rate sinonime între codoni. O astfel de autocorelare ar fi de așteptat dacă selecția sau variația ratei mutației ar fi localizate spațial în ORF. Parametrul de comutare a ratei între codonii adiacenți ai HMM descrie amploarea autocorelației, cu valori sub 1 / N (N = numărul de clase de viteză) care sugerează autocorelația. Tehnicile HMM standard (de exemplu, calea Viterbi) aplicate acestor modele pot dezvălui unde apar comutările între diferite tipuri de rate, împărțind astfel secvența în regiuni cu constrângeri mai slabe sau mai puternice pe substituții sinonime.

Metoda aBSREL [41] a fost utilizată pe toate ramurile cladei nCoV pentru a determina care ramuri specifice determină inferența selecției. În cele din urmă, am examinat care site-uri specifice de codoni sunt sub selecție negativă în medie pe nadei nCoV utilizând FEL [30] și sub selecție pozitivă sau diversificată episodică pozitivă pe nadei nCoV folosind MEME [36]. P valori de ≤0,05 pentru testele raportului de probabilitate, specifice fiecărei metode, au fost luate ca dovadă a semnificației statistice. Majoritatea analizelor de selecție au fost efectuate în pachetul software HyPhy v.2.5.14, cu BUSTED ajustat pentru a permite substituții multiple de nucleotide folosind v2.5.24 [72].

Epuizarea CpG

Pentru a cuantifica supra / sub reprezentarea dinucleotidelor CpG în Sarbecovirus genomuri, am dezvoltat o versiune modificată a metricii sinonime a utilizării dinucleotidelor (SDU) [55], care reprezintă acum compoziția de bază părtinitoare. Valoarea originală SDU compară proporția observată de CpG sinonim, o pentru fiecare pereche de poziții de cadre într-o secvență de codare, h la cea așteptată în condiții de utilizare egală a codonilor sinonimi, e pentru fiecare aminoacid (sau pereche de aminoacizi) care poate avea codoni care conțin CpG (sau perechi de codoni), eu. Valoarea SDU este media acestor rapoarte ponderate cu numărul de aminoacizi (sau perechi) informative din secvență, n (Ec. 1).

Pentru a încorpora compoziția de bază părtinitoare și variabilă a SARS-CoV-2 și altele Sarbecovirusuri [47], aici am estimat utilizarea preconizată a codonilor pe baza compoziției nucleotidice a genomului întreg al fiecărui virus. Denumim această nouă metrică utilizarea sinonimă corectată a dinucleotidelor (SDUc). Folosim frecvențele de bază observate de la fiecare virus pentru a genera așteptarea nulă corectată a metricei, e ′, în loc să presupunem o utilizare egală (ecuația 1). Proporția așteptată, e ′, pentru fiecare pereche de aminoacizi / aminoacizi a fost estimată prin simularea aleatorie a codonilor pe baza proporțiilor nucleotidice unice ale întregului genom al fiecărui virus. Acest e ′ a fost apoi utilizat pentru toate calculele SDUc ale virusului corespunzător.

Deoarece această valoare este susceptibilă de eroare atunci când este utilizată pentru secvențe de codare scurte, am aplicat SDUc pe cel mai lung ORF (Orf1ab) dintre toți virusii. Pentru a estima amploarea independenței filogenetice între site-urile sinonime din punctele de date SDUc, am măsurat divergența sinonimă în perechi (Ks) între viruși. Valorile Ks în perechi au fost calculate utilizând pachetul R seqinr [73], care utilizează modelul de codoni ai lui Li (1993) [74], demonstrând independența parțială, dar nu completă, în cadrul celor două linii. Mediana și maximul Ks sunt 0,54 și 0,89 în cadrul nadei nCoV și respectiv 0,34 și 1,09 în cadrul nadei nCoV / SARS-CoV-1. Eq 1 (Eq 1, N = numărul total de aminoacizi informativi)

Analiza recombinării vârfurilor

Pentru a determina punctul de întrerupere a recombinării pe Spike ORF al virusului RmYN02, am folosit suita de metode RDP5 [75], implementând 7 metode: RDP, GENECONV, Chimaera, MaxChi, BootScan, SiScan și 3seq. Am efectuat mai întâi analiza asupra alinierii genomului întreg al Sarbecovirusuri și apoi a determinat punctul de întrerupere relevant în cadrul Spike ORF reanalizând metoda pe alinierea numai Spike. Punctul de întrerupere acceptat (poziția 24058 în genomul RmYN02) a fost apelat în mod constant de 6 din cele 7 metode testate (RDP, GENECONV, Maxchi, Chimaeara, SiSscan și 3seq). În mod similar, punctul de întrerupere precedent al regiunii non-nCoV a fost chemat la poziția 21248 a genomului RmYN02 (înainte de începerea Spike ORF).

Inferință bayesiană sub un ceas molecular local

Pentru a evalua relația filogenetică generală nerecombinantă dintre viruși, am folosit cele mai lungi regiuni genomice nerecombinante descrise în Boni și colegii [6], NRR1 și NRR2 și am adăugat RmYN02 la aliniamente. Pe baza punctelor de întrerupere de recombinare determinate mai sus pentru RmYN02, NRR2 a fost ajustat pentru a se încheia la poziția 21266 a Wuhan-Hu-1, în loc de 21753, corespunzător cu începutul regiunii recombinante RmYN02 (poziția RmYN02 21248). Istoriile evolutive măsurate în timp pentru NRR1 și NRR2 au fost deduse utilizând o abordare bayesiană, implementată prin cadrul Markov lanțului Monte Carlo (MCMC) disponibil în BEAST 1.10 [58]. Motivați de observarea unei divergențe mai mari de la rădăcină la vârf pentru clada nCoV (a se vedea textul S3 și fig. S5) și de deplasarea trăsăturii CpG pe ramura ancestrală la clada nCoV (a se vedea mai jos), am specificat un ceas local fix [76] care permite o rată diferită pe ramura care duce la clada nCoV. În absența unui semnal temporal puternic, am specificat o distribuție anterioară informativă normală (cu medie = 0,0005 și deviație standard = 0,0002) pe rata tuturor celorlalte ramuri pe baza estimărilor recente sub un ceas molecular relaxat [77]. Pentru a menține topologia generală sub modelul ceasului local, am restricționat virusul liliecian kenyan (KY352407) ca grup outgroup pentru virușii din China în NRR1 și kenianul, precum și un virus liliec bulgar (NC_014470) ca grup out pentru virusurile din China în NRR2. Am partiționat regiunile de codare ale NRR1 și NNR2 după poziția codonului și am specificat un model independent independent de substituție reversibil în timp (GTR) cu variație a ratei de distribuție gamma între site-uri pentru fiecare dintre cele 3 partiții. Am folosit un model coalescent de mărime constantă ca arbore anterior și am specificat un lognormal anterior cu medie = 6,0 și deviație standard = 0,5 pe dimensiunea populației. Trei analize MCMC independente au fost efectuate pentru 250 de milioane de state pentru fiecare set de date. Am folosit biblioteca BEAGLE v3 [78] pentru a crește performanța de calcul. Pentru completare, am efectuat aceeași analiză specificând modelul de ceas local pe întreaga cladă nCoV în loc doar de ramura care duce la aceasta. Acest lucru are un efect asupra ratei de substituție și estimările timpului nodului, cu toate acestea, nu putem distinge în mod formal care model se potrivește mai bine (text S3). Ne-am concentrat pe modelul de ceas local numai ramură, deoarece ne testează direct ipoteza că schimbarea ratei este legată de schimbarea adaptivă CpG specifică acelei ramuri. Fig. 3A prezintă filogenia pentru NRR2 deoarece este cea mai lungă regiune nerecombinantă intactă și ar putea produce estimări mai fiabile. Filogeniile NRR1 suplimentare și filogeniile modelului local al nadei Clade sunt prezentate în textul S3. Fișierele XML ale parametrilor BEAST sunt furnizate în datele S1 (alinierea secvenței este exclusă pentru a respecta restricțiile de partajare a datelor GISAID). Parametrii continui au fost rezumați și dimensiunile eficiente ale eșantionului au fost estimate folosind Tracer [79]. Estimările timpului de divergență pentru clada nCoV sunt rezumate în textul S3. Arborii au fost rezumați ca arbori de credibilitate maximă a cladei (MCC) folosind TreeAnnotator și vizualizați utilizând FigTree (http://tree.bio.ed.ac.uk/software/figtree/).

Identificarea schimbărilor în conținutul CpG

Schimbările în conținutul de CpG au fost identificate folosind o metodă comparativă filogenetică care deduce schimbări adaptive în trăsături corelate multivariate cu pachetul R PhylogeneticEM [56]. Această abordare modelează evoluția trăsăturilor pe filogenii folosind un proces Ornstein – Uhlenbeck (OU) și folosește o versiune tratabilă pe bază de calcul a modelului complet OU multivariant (OU scalar) pentru trăsături multivariate. Estimările pozițiilor de schimbare sunt obținute utilizând un algoritm de așteptare-maximizare (EM). Pozițiile de schimbare sunt estimate pentru diferite numere de schimbări necunoscute, iar o procedură de selecție a modelului lazo-regresiv identifică numărul optim de schimbări. Am aplicat procedura copacilor MCC pentru NRR1 și NRR2 cu valorile lor CpG SDUc respective (ln transformat după cum este cerut de abordarea filogenetică). Arborii au fost forțați să fie ultrametrici, așa cum este necesar pentru modelul UO scalar, prin extinderea tuturor marginilor externe ale copacilor pentru a se potrivi cu vârful cel mai recent eșantionat. Folosind această procedură, am identificat 3 și 2 schimbări CpG atât în ​​NRR1, cât și în NRR2, respectiv (S6 Fig), cu deplasarea pe ramură ancestrală la clada nCoV fiind singura consecventă identificată în ambele regiuni genomice (Fig 3A).


Rezultate

Evaluarea datelor simulate

Pe toate datele, am dedus cel mai probabil model, care a dat estimarea ωși apoi ne-am folosit abordarea pentru a estima dN, dS, și d N d S ⁠.

Cand ω se estimează cu un model staționar, scăderea conținutului de G + C de-a lungul arborelui are ca rezultat o supraestimare sistematică a ω, și creșterea conținutului G + C are ca rezultat o subestimare sistematică a ω (vezi fig. suplimentară S1 A, Material suplimentar online). Observăm părtiniri similare în estimările lui d N d S estimate cu mapare stocastică (fig. 1e).

Estimări de dN, dS, și d N d S cu un model staționar (stânga) și cu un model nestacionar (dreapta), pe date simulate cu conținut modificat de G + C și ω = 0,1 ⁠. θrădăcină: Frecvența G + C în secvența rădăcină. θechiv: Frecvența de echilibru G + C a modelului de simulare.

Estimări de dN, dS, și d N d S cu un model staționar (stânga) și cu un model nestacionar (dreapta), pe date simulate cu conținut modificat de G + C și ω = 0,1 ⁠. θrădăcină: Frecvența G + C în secvența rădăcină. θechiv: Frecvența de echilibru G + C a modelului de simulare.

Aceste subestimări sau supraevaluări pot duce la o interpretare calitativă falsă a selecției, întrucât selecția pozitivă dubioasă poate fi dedusă în cazul scăderii conținutului GC sau selecția negativă dubioasă în cazul creșterii conținutului GC (așa cum este ilustrat în simulările cu modele neutre și aproape neutre , vezi figurile suplimentare. S2 – S4, Material suplimentar online). De asemenea, am efectuat simulări în care conținutul G + C al unei poziții specifice de codoni a evoluat, iar celelalte două au rămas staționare cu 50% G + C. Din nou, am observat că modelele care presupun staționaritate conduc la estimări părtinitoare ale ratelor de substituție și ale ω (vezi figurile suplimentare. S6 – S8, Material suplimentar online). Interesant este că orientarea prejudecății este diferită indiferent dacă poziția de schimbare a G + C este a treia (adică cea mai sinonimă) sau nu. Prin urmare, diferite combinații de modificări specifice poziției G + C pot duce la diferite tipuri de părtiniri.

În ceea ce privește ratele de substituție (non) sinonime, presupunând o prejudecată de staționaritate atât estimările dN și dS în moduri similare (fig. 1a și c). Aceste valori sunt în mare parte subestimate atunci când echilibrul GC este foarte diferit de 0,5 și conținutul GC se modifică (fie în sus, fie în jos) (fig. 2). Astfel, în aceste cazuri, arborii deduși sunt prea scurți.

Raportul ratelor de substituție estimate cu modelul staționar față de ratele de substituție estimate cu modelul non-staționar. Secvențele au fost simulate cu modificarea conținutului G + C și ω = 0,1 ⁠. Stânga: dN. Dreapta: dS. θrădăcină: Frecvența G + C în secvența rădăcină. θechiv: Frecvența de echilibru G + C a modelului de simulare.

Raportul ratelor de substituție estimate cu modelul staționar față de ratele de substituție estimate cu modelul non-staționar. Secvențele au fost simulate cu modificarea conținutului G + C și ω = 0,1 ⁠. Stânga: dN. Dreapta: dS. θrădăcină: Frecvența G + C în secvența rădăcină. θechiv: Frecvența de echilibru G + C a modelului de simulare.

Pentru a verifica dacă aceste prejudecăți nu se datorează metodei noastre, am efectuat, de asemenea, aceleași estimări sub ipoteza staționară cu codeml (Yang și Nielsen 2000), iar rezultatele prezintă părtiniri similare (vezi fig. S9 suplimentar, Material suplimentar online).

Toate aceste prejudecăți sunt corectate atunci când se utilizează abordarea noastră cu un model nestatiar, atât pentru ω (fig. suplimentar S1 b, Material suplimentar online), dN, dS, și d N d S (fig. 1b, d și f), și chiar și atunci când nonstationarity diferă între pozițiile codonilor (vezi figurile suplimentare. S6 – S8, Material suplimentar online).

Interesant, am observat că estimările dN și dS scade cu conținutul GC de echilibru (fig. 1b și d). Acest lucru nu se datorează metodei noastre, deoarece pentru procesele staționare, estimările dN și dS calculat cu codeml prezintă o tendință similară (a se vedea linia punctată din fig. S9 suplimentar, Material suplimentar online). Această relație dintre dN sau dS iar echilibrul conținutului GC depinde de valoarea omega. Cand ω este scăzută, această corelație este negativă (fig. 1b și d), când ω este egal cu 1, corelația este nulă (fig. S3 suplimentară, Material suplimentar online) și devine pozitivă ca ω devine mai mare de 1 (de exemplu, pentru ω = 2 vezi fig. Suplimentar. S5, Material suplimentar online).

Trebuie remarcat faptul că, atunci când dinamica conținutului de GC este eterogenă, tendința nu este sistematic în aceeași direcție indiferent dacă crește (sau scade) GC, ci depinde și de CG de alte ramuri, deoarece o modelare staționară (deci omogenă) va estimează echilibrul său GC din toate ramurile. De exemplu, pe același copac, am luat în considerare un model cu GC staționar de la rădăcină la frunzele primatului și schimbarea GC pe ramurile care duc la câine și la rozătoare. Așa cum se arată în figura 3, estimările lui d N d S pe ramurile primatelor sunt părtinitoare cu ipoteza staționarității, chiar dacă procesul este într-adevăr staționar pe aceste ramuri. Dar non-staționalitatea de pe celelalte ramuri induce în eroare modelul staționar estimat.

Estimarea dN, dS, și d N d S pe ramurile primatelor cu un model staționar (stânga) și un model neomogen nestatiar (dreapta), pe date simulate cu conținut modificat de G + C pe ramurile câinilor și rozătoarelor și ω = 0,1 ⁠. θrădăcină: Frecvența G + C în secvențele rădăcină și primate. θechiv: Frecvența de echilibru G + C a modelului de simulare pe ramuri de câine și rozătoare.

Estimarea dN, dS, și d N d S pe ramurile primatelor cu un model staționar (stânga) și un model neomogen nestatiar (dreapta), pe date simulate cu conținut modificat de G + C pe ramurile câinilor și rozătoarelor și ω = 0,1 ⁠. θrădăcină: Frecvența G + C în secvențele rădăcină și primate. θechiv: Frecvența de echilibru G + C a modelului de simulare pe ramuri de câine și rozătoare.

Studiu asupra setului de date despre mamifere

Am efectuat două estimări diferite de probabilitate maximă a setului de date la mamifere: un model staționar omogen YN98 + F3X4 (21 parametri de ramură și model) și un model nestanțiar neomogen (31 de parametri suplimentari) cu trei modele omogene YN98, unul pentru clada primatelor, una pentru clada rozătoarelor și una pentru ramura câinelui. Am folosit trei modele pentru a se potrivi cu heterogenitatea conținutului GC de echilibru găsit între aceste clade Romiguier și colab. (2010). Am calculat dN (resp. dS) în clada primatelor prin însumarea mapării stochastice dN (resp. dS) din toate ramurile acestei clade.

Deoarece modelele sunt imbricate, am efectuat teste ale raportului de probabilitate la toate estimările și am corectat testele multiple folosind corecția Benjamini-Hochberg. Creșterea probabilității este semnificativă (folosind un test LRT cu 31 de grade de libertate) cu o valoare FDR de 1%, la 83,4% din gene (fig. 4).

Log10 al diferențelor în probabilitățile de logare între modelele staționare și non-staționare pe datele mamiferelor. Linia roșie reprezintă pragul de 5% FDR.

Log10 al diferențelor în probabilitățile de logare între modelele staționare și non-staționare pe datele mamiferelor. Linia roșie reprezintă pragul de 5% FDR.

Dacă comparăm estimările modelării staționare versus non-staționare, vedem că estimările dN sunt în mare parte mai mici, dar nu corelate cu evoluția conținutului GC la poziția a treia codon (GC3) (fig. 5). Dimpotrivă, vedem o influență a evoluției în GC3 asupra prejudecății în estimarea lui dS, și apoi o subestimare mai importantă a d N d S cu gene departe de staționaritate în GC3. După cum sa observat în secțiunea de simulare, tendința nu este corelată cu semnul schimbării în GC3, deoarece am efectuat o modelare neomogenă, iar tendința depinde și de evoluția conținutului de GC în celelalte ramuri. Cu toate acestea, efectul este destul de vizibil: eroarea relativă la ω estimarea este de cel puțin 10% pentru 59% din gene sau de cel puțin 33% pentru 13,4% din gene (fig. 6).

log2 al raporturilor estimărilor lui dN, dSși dN/ dS cu un model staționar peste estimări cu un model nestacionar, în funcție de modificarea GC3 conținut în clada primatelor.

log2 al raporturilor estimărilor lui dN, dSși dN/ dS cu un model staționar peste estimări cu un model nestacionar, în funcție de modificarea GC3 conținut în clada primatelor.

Histograma raporturilor în estimări de ω în modelul staționar față de modelul non-staționar pe date de mamifere. Liniile galbene, portocalii și roșii reprezintă 12,5%, 25% și 37,5% cuantile. Linia purpurie reprezintă mediana.

Histograma raporturilor în estimări de ω în modelul staționar față de modelul non-staționar pe date de mamifere. Liniile galbene, portocalii și roșii reprezintă 12,5%, 25% și 37,5% cuantile. Linia purpurie reprezintă mediana.


Metode și materiale

Prezentare generală a procesului de simulare

ALF generează un set de genomi de specii pornind de la o singură secvență de genom ancestral. Genomul ancestral poate fi reprezentat de secvențe biologice furnizate de utilizator sau generate aleator în conformitate cu specificațiile utilizatorului. Un arbore de specie poate fi, de asemenea, specificat de utilizator sau generat aleatoriu. În cursul simulării, ALF evoluează genomul rădăcinii de-a lungul arborelui, unde fiecare nod definește un eveniment de speciație. Genomurile emergente sunt expuse proceselor evolutive implementate în ALF.

Figura 1 oferă o prezentare grafică a conductei de simulare ALF. Substituțiile de caractere apar în funcție de matricea probabilității de substituție a unui model selectat de aminoacizi, codoni sau nucleotide pentru o lungime de ramură dată. Pot fi specificate diferite modele pentru simulare, de exemplu, un model de codon și unul de nucleotide ar putea fi utilizate pentru a distinge regiunile de codificare și respectiv necodificate. Rata substituției poate diferi în funcție de site-uri și gene. ALF permite fiecărei specii să aibă propriile sale frecvențe de bază de echilibru, de exemplu, pentru a simula derivarea către conținutul GC specific speciei.

Prezentare generală a procesului de simulare ALF. Un genom al rădăcinii este dezvoltat de-a lungul unui copac al speciei. Evenimentele de la locul, secvența și nivelul genomului sunt simulate iterativ.


Ratele de mutație a ADN-ului mitocondrial

O problemă este că ratele de mutație nu sunt cunoscute prin măsurare directă și sunt adesea calculate pe baza unor scale de timp evolutive presupuse. Astfel, toate aceste estimări de vârstă ar putea fi foarte greșite. De fapt, multe rate diferite de mutație sunt citate de diferiți biologi.

Nu ar trebui să fie foarte greu să se măsoare în mod explicit rata mutației ADN-ului mitocondrial pentru a obține o estimare mai bună a acestei vârste. De exemplu, din descendențele regale, s-ar putea găsi doi indivizi al căror cel mai recent strămoș matern comun a fost, să zicem, acum 1000 de ani. S-ar putea măsura apoi diferențele ADN-ului mitocondrial al acestor indivizi pentru a-și lega rata de mutație. Această schemă este atractivă, deoarece nu depinde de datarea radiometrică sau de alte ipoteze cu privire la ratele de evoluție sau mutație. Este posibil ca în 1000 de ani să existe o diferență prea mică de măsurat. Cel puțin acest lucru ne-ar oferi în continuare câteva informații utile.

(Un proiect pentru oamenii de știință de creație!)

Pe această linie, s-au făcut recent unele lucrări pentru a măsura în mod explicit rata de substituție în ADN-ul mitocondrial. Referința este Parsons, Thomas J. și colab., O rată ridicată de substituție observată în regiunea de control a ADN-ului mitocondrial uman, Nature Genetics vol. 15, aprilie 1997, pp. 363-367. Rezumatul urmează:

„Rata și modelul substituțiilor de secvențe în regiunea de control ADN mitocondrial (ADNmt) sunt de importanță centrală pentru studiile evoluției umane și pentru testarea identității medico-legale. Aici, raportăm o măsurare directă a ratei de substituție intergenerațională la om CR. Am comparat secvențe de ADN a două segmente hipervariabile CR de la rude materne apropiate, de la 134 linii de ADNmt independente care se întind pe 327 de evenimente generaționale. Au fost observate zece substituții, rezultând o rată empirică de 1/33 generații sau 2,5 / site / Myr. este de aproximativ douăzeci de ori mai mare decât estimările derivate din analizele filogenetice. Această diferență nu poate fi explicată doar prin substituții la punctele fierbinți mutaționale, sugerând factori suplimentari care produc discrepanța între ratele aparente de foarte scurtă durată și cele pe termen lung ale divergenței secvenței. datele indică, de asemenea, că segregarea extrem de rapidă a variantelor de secvență CR între generații este comună la om, cu un mtDN foarte mic Un blocaj. Aceste rezultate au implicații pentru aplicațiile criminalistice și studiile evoluției umane. "(Op. Cit. P. 363).

"Rata de substituție observată raportată aici este foarte mare în comparație cu ratele deduse din studiile evolutive. O gamă largă de rate de substituție CR au fost derivate din studii filogenetice, care se întind pe aproximativ 0,025-0,26 / sit / Myr, inclusiv intervale de încredere. Un studiu care produce unul dintre estimările mai rapide au dat rata de substituție a regiunilor hipervariabile CR ca 0.118 + - 0.031 / sit / Myr. Presupunând un timp de generație de 20 de ani, aceasta corespunde

1/600 de generații și o vârstă pentru ADNmt MRCA de 133.000 y.a. Astfel, observația noastră a ratei de substituție, 2,5 / sit / Myr, este de aproximativ 20 de ori mai mare decât s-ar prevedea din analizele filogenetice. Folosirea ratei noastre empirice pentru calibrarea ceasului molecular mtDNA ar avea ca rezultat o vârstă a ADNmt MRCA de numai

6.500 de ani, în mod clar incompatibil cu epoca cunoscută a oamenilor moderni. Recunoscând chiar că MRCA al ADNmt poate fi mai tânăr decât MRCA al oamenilor moderni, rămâne neverosimil să explicăm distribuția geografică cunoscută a variației secvenței mtDNA prin migrația umană care a avut loc doar în ultima

Un biolog a explicat estimarea vârstei tinere presupunând în esență că 19/20 dintre mutațiile din această regiune de control sunt ușor dăunătoare și în cele din urmă vor fi eliminate din populație. Acest lucru pare puțin probabil, deoarece această regiune tinde să varieze foarte mult și, prin urmare, are probabil o funcție mică. În plus, dezavantajul selectiv al acestor 19/20 dintre mutații ar trebui să fie de aproximativ 1/300 sau mai mare pentru a evita producerea unei divergențe mai mari în secvențe decât observată în mai mult de 6000 de ani. Aceasta înseamnă că unul din 300 de indivizi ar trebui să moară din cauza unor mutații în această regiune. Aceasta pare a fi o cifră ridicată pentru o regiune care pare să fie în mare parte fără funcții. Este interesant faptul că același biolog consideră că 9/10 dintre mutațiile regiunilor codificatoare ale ADN sunt neutre. Acest lucru face ca regiunile de codificare ale ADN-ului să fie mai puțin constrânse decât regiunea de control aparent fără funcție a ADN-ului mitocondrial!


Materialul suplimentar electronic este disponibil online la https://dx.doi.org/10.6084/m9.figshare.c.4084511.

Publicat de Royal Society. Toate drepturile rezervate.

Referințe

. 1991 Evoluția proteinelor adaptive la Adh locus in Drosophila . Natură 351, 652-654. (doi: 10.1038 / 351652a0) Crossref, PubMed, ISI, Google Scholar

Eyre-Walker A, Keightley PD

. 2009 Estimarea ratei evoluției moleculare adaptive în prezența mutațiilor ușor dăunătoare și a modificării dimensiunii populației. Mol. Biol. Evol. 26, 2097-2108. (doi: 10.1093 / molbev / msp119) Crossref, PubMed, ISI, Google Scholar

. 2013 Adaptare frecventă și testul McDonald – Kreitman. Proc. Natl Acad. Știință. Statele Unite ale Americii 110, 8615-8620. (doi: 10.1073 / pnas.1220835110) Crossref, PubMed, ISI, Google Scholar

. 2016 Evoluția adaptivă a proteinelor la animale și ipoteza efectivă a dimensiunii populației. PLoS Genet. 12, e1005774. (doi: 10.1371 / journal.pgen.1005774) Crossref, PubMed, ISI, Google Scholar

Tataru P, Mollion M, Glémin S, Bataillon T

. 2017 Inferența distribuției efectelor de fitness și proporția substituțiilor adaptive din datele polimorfismului. Genetica 207, 1103–1119. (doi: 10.1534 / genetics.117.300323) Crossref, PubMed, ISI, Google Scholar

. 2002 Schimbarea dimensiunii efective a populației și testul McDonald – Kreitman. Genetica 162, 2017–2024. PubMed, ISI, Google Scholar

. 1962 Despre probabilitatea fixării genelor mutante într-o populație. Genetica 47, 713–719. PubMed, ISI, Google Scholar

. 1999 Recolonizarea post-glaciară a biotei europene. Biol. J. Linn. Soc. 68, 87-112. (doi: 10.1111 / j.1095-8312.1999.tb01160.x) Crossref, ISI, Google Scholar

. 2001 Refugia criptică nordică și originile biotei moderne. Tendințe Ecol. Evol. 16, 608-613. (doi: 10.1016 / S0169-5347 (01) 02338-2) Crossref, ISI, Google Scholar

. 2006 Răspunsuri ecologice și evolutive la schimbările climatice recente. Annu. Pr. Ecol. Evol. Syst. 37, 637-669. (doi: 10.1146 / annurev.ecolsys.37.091305.110100) Crossref, ISI, Google Scholar

. 2017 SLiM 2: simulări genetice interactive flexibile, interactive. Mol. Biol. Evol. 34, 230-240. (doi: 10.1093 / molbev / msw211) Crossref, PubMed, ISI, Google Scholar

. 2002 Ratele de mutație la genomii de mamifere. Proc. Natl Acad. Știință. Statele Unite ale Americii 99, 803-808. (doi: 10.1073 / pnas.022629899) Crossref, PubMed, ISI, Google Scholar

Stapley J, Feulner PG, Johnston SE, Santure AW, Smadja CM

. 2017 Variația frecvenței de recombinare și distribuția între eucariote: modele și procese. Phil. Trans. R. Soc. B. 372, 20160455. (doi: 10.1098 / rstb.2016.0455) Link, ISI, Google Scholar

Eyre-Walker A, Woolfit M, Phelps T

. 2006 Distribuția efectelor de fitness ale noilor mutații nocive ale aminoacizilor la om. Genetica 173, 891–900. (doi: 10.1534 / genetics.106.057570) Crossref, PubMed, ISI, Google Scholar

Rousselle M, Mollion M, Nabholz B, Bataillon T, Galtier N

. Date din 2018 din: supraestimarea ratei de substituție adaptivă la populațiile fluctuante. Depozitul digital Dryad. (doi: 10.5061 / dryad.85qb2r1) Google Scholar


Priveste filmarea: Naturfag - Genetikk, arv og miljø (Ianuarie 2022).