Studio di associazione Genome-Wide: Preprocessing e Selezione SNPs (PDF)

UNIVERSITA’ DEGLI STUDI DI PADOVA FACOLTA’ DI INGEGNERIA Corso di laurea magistrale in Bioignegneria Studio di associazione Genome Wide: Preprocessing e Selezione SNPs LAUREANDA : Giulia Bruscagin RELATORE: prof Barbara di Camillo CORRELATORE: dott. Francesco Sambo Anno Accademico 2010/2011 RINGRAZIAMENTI Colgo l’occasione per ringraziare i dottori Barbara di Camillo e Francesco Sambo, per avermi guidato nella realizzazione di tutto il lavoro, per i consigli e per la disponibilità. Ringrazio anche il dottor Emanuele Trifoglio per avermi aiutato a svolgere gran parte del lavoro di questa tesi e per la disponibilità dimostrata. Un ringraziamento alla mia famiglia che mi ha sostenuto e incoraggiato in tutti gli anni di studio, e alle mie compagne Emanuela, Sara, Marta e Valentina per i cinque bellissimi anni passati insieme. This study makes use of data generated by Wellcome Trust Case-Control Consortium. A full list of the investigators who contributed to the generation of the data is avaible from www.wtccc.org.uk. Funding for the project was provided by Wellcome Trust under award 076113 and 085475’. Indice Sommario Introduzione Capitolo 1: Polimorfismi a singolo nucleotide e test di associazione…………..……….1 1.1 DNA e Polimorfismi di un singolo nucleotide………………………………………1 1.1.1 Frequenza degli SNPs e Minor Allele Frequency…………………………3 1.2 Linkage Disequilibrium ……………………………………………………………..4 1.2.1 Variazioni di LD nel genoma………………………………………………5 1.2.2 Aplotipi e blocchi di aplotipi………………………………………………7 1.3 Equilibrio Hardy-Weinberg………………………………………………………….7 1.4 Studi di associazione...…………………………………………………………….10 Capitolo 2: Linkage Disequilibrium……………………………..…………………..…11 2.1 Misure di Linkage Disequilibrium..……………………………………………….11 2.1.1 Misure di LD Pairwise……………………………………………………11 2.1.2 Misure di LD Multilocus...………………………………………………17 2.2 Definizione di Aplotipi mediante misure di LD…………………………………...19 2.3 Algoritmi di LD e inferenza di Aplotipi...…………………………………………21 2.4 TagSNP...…………………………………………………………………………..28 Capitolo 3 : Stratificazione di Popolazione ……………………………………………35 3.1 Cause della stratificazione …………………………………………………………35 3.2 Identificazione della stratificazione ………………………………………………..36 3.2.1 Genomic Control (GC)...………………………………………………...38 3.2.2 Principal Component Analysis (PCA)……………………………………38 3.2.3 Mixed Models…………………………………………………………….40 Capitolo 4 : Il progetto internazionale HapMap …………….…………………………41 4.1 Cos’è HapMap……………………………………………………………………...41 4.2 Realizzazione del progetto HapMap….…………………………………………….42 4.3 Pubblicazione e Consultazione dei dati...………………………………………….44 4.3.1Esempio Ricerca di TCF7L2...…………………………...………………45 4.4 HapMart…………………………………………………………………………….53 4.5 Ricerca per malattia...……………………………………………………………...54 Capitolo 5: PLINK ……...……………………………………………………………...55 5.1 Formato dei dati PLINK……………………………………………………………55 5.1.1 PED files …………………………………………………………………56 5.1.2 MAP files ………………………………………………………………...57 5.1.3 Fileset trasposti …………………………………………………………..58 5.1.4 File binari ………………………………………………………………...59 5.1.5 Codifica cromosomi e alleli ……………………………………………...59 5.2 Data management ………………………………………………………………….59 5.3 Summary Statistic ………………………………………………………………….60 5.3.1 Missingness ………………………………………………………………60 5.3.2 Hardy Weinberg Equilibrium ……………………………………………61 5.3.3 Frequenza allelica...……………………………………………………...62 5.4 Analisi di stratificazione di popolazione …………………………………………..62 5.4.1 Definizione di similarità e distanza tra individui ………………………..63 5.4.2 Vincoli sul clustering …………………………………………………….63 5.4.3 Algoritmo di clustering …………………………………………………..64 5.4.4 Scaling Multidimensionale (MDS) ………………………………………65 5.4.5 Individuazione di individui outliers………………………………………66 5.5 Analisi di Associzione ……………………………………………………………..67 5.6 Stima IBD (identical by descend)………………………………………………….68 Capitolo 6: Obiettivi della Tesi ………………………………………………………..71 6.1 Motivazioni ………………………………………………………………………...71 6.2 Strategie ……………………………………………………………………………72 Capitolo 7: Dati e Preprocessing...….…………………………………………………75 7.1 Descrizione dei dati ………………………………………………………………..75 7.2 Preprocessing dei dati………………………………………………………………76 Capitolo 8: Definizione di Metavariabili basata sulla MI ……………………………..79 8.1 Calcolo di entropia, mutua informazione e definizione di metaviariabili………….80 8.2 Classificazione con Naive Bayes …………………………………………………..81 Capitolo 9 : Risultati …………………………………………………………………...83 9.1 Rete ricostruita dai dati di controlli e casi del diabete di tipo 1 …………………...84 9.2 Rete ricostruita dai dati di controlli e casi del diabete di tipo 2 ………………….128 9.3 Risultati della classificazione naive Bayes ……………………………………….163 Conclusioni …………………………………………………………………………...165 Bibliografia …………………………………………………………………………...166 Sommario Gli studi di associazione intervengono nello studio di dataset casi/controllo caratterizzati da un’importante mole di dati. Tipicamente un dataset è costituito da un numero di soggetti per classe dell’ordine delle migliaia e da un numero di variabili, i polimorfismi a singolo nucleotide, nell’ordine di 108 È necessario operare una riduzione sensata del numero di variabili iniziali per poter operare efficientemente una classificazione. Dopo un’attenta analisi dello stato dell’arte, e dopo aver evidenziato i limiti principali dovuti essenzialmente alla feature selection e alla riduzione di informazione che ne consegue, in questa tesi si propone un nuovo approccio basato sulla definizione di mutua informazione per la definizione di metavariabili, da utilizzarsi poi nella classificazione. La metodologia viene applicata con successo su due distinti dataset, isolando in prima battuta il pathway dell’insulina e procedendo alla classificazione delle metavariabili ricostruite. L’applicazione del metodo nella sua completezza prevede la costruzione di un classificatore aggregato dei singoli classificatori costruiti su diversi pathway biologici. Introduzione I recenti sviluppi nelle tecnologie del sequenziamento del genoma hanno permesso l’applicazione a livello di popolazione su larga scala dei test di associazione genetica sfruttando in particolare uno specifico tipo di marcatori: i polimorfismi a singolo nucleotide, o SNPs, i quali contribuiscono a regolare la predisposizione a una determinata patologia. I test in questione vengono denominati Test di Associazione Genome Wide, in quanto, come base di partenza, considerano la totalità di dati a disposizione dal sequenziamento degli SNPs (circa 10 milioni nel genoma umano). L’obiettivo che questi studi si propongono è quello di effettuare una classificazione tra campioni di soggetti che si dividono in controlli (soggetti sani) e casi (soggetti affetti da una specifica patologia che ha base genica) sulla base delle differenze di frequenza allelica con cui gli SNPs si presentano. Nella prima parte della tesi, dopo un breve inquadramento della tematica da un punto di vista biologico, viene proposto lo stato dell’arte che riguarda i test di associazione facendo particolare attenzione a come i ricercatori si pongono davanti alla grande mole di dati che si hanno a disposizione in questo genere di studi. Se da una parte un dataset estremamente informativo può costituire un enorme vantaggio, dall’altra, un enorme numero di variabili costituisce un limite all’applicazione di qualsiasi algoritmo di classificazione. È’ necessario ridurre il numero di variabili in gioco. Dallo stato dell’arte, si osservano due procedure alternative (individuazione di tagSNPs, o selezione univariata e successivamente multivariata) che sicuramente portano a una riduzione delle variabili da considerare, ma allo stesso tempo riducono pesantemente il contenuto informativo del dataset iniziale. L’idea della tesi nasce proprio da questa osservazione e si pone come obiettivo quello di ridurre il numero di varaibili ma mantenendo la quantità globale di informazione intatta. Per farlo, si ricorre alle definizioni di entropia e mutua informazione che andranno a stabilire il criterio con cui costruire le metavariabili. La tesi, rispetto allo stato dell’arte, propone un’altra novità: infatti, a differenza delle nomali procedure che considerano separatamente i cromosomi nello studio di determinate regioni geniche, qui viene considerato un pathway nel suo insieme. Si ritiene infatti che, in questo modo, sia possibile individuare, se esiste, una rete di regolazione tra diversi geni di un unico pathway che hanno la caratteristica di intervenire nella regolazione di un unico prodotto. Questa strategia è stata applicata a due dataset, rilasciati dalla WTCCC, costituiti da controlli e casi di diabete di tipo 1 e 2 rispettivamente. La tesi si compone di 8 capitoli. Nel primo capitolo vengono brevemente forniti alcuni concetti base relativamente alla biologia degli argomenti trattati. Nel secondo capitolo si approfondisce il fenomeno del Linkage Disequilibrium, le misure ad esso associato e lo stato dell’arte dei principali algoritmi disponibili per la definizione di aplotipi e tag SNP. Nel terzo capitolo viene presentato in generale il problema della stratificazione di popolazione, ricorrente nella maggior parte degli studi di associazione, e vengono riportati i principali metodi per la sua risoluzione. Il quarto capitolo è dedicato alla descrizione di uno dei più importanti database disponibile gratuitamente in rete per la consultazione degli SNPs: HapMap. Viene brevemente riportata la storia del progetto, le sue finalità e i metodi di consultazione. Il quinto capitolo descrive il programma PLINK utilizzato per l’analisi preliminare sui dati WTCCC, per la mappatura dei marcatori e il preprocessing. Ne vengono brevemente illustrate le principali funzionalità. Il sesto, settimo e ottavo capitolo sono infine dedicati alla descrizione dei dati disponibili, al metodo implementato e ai risultati ottenuti dalla mappatura della rete, ricostruita con le metavariabili, e dalla classificazione ottenuta. Capitolo 1 Polimorfismi a singolo nucleotide e test di associazione In questo capitolo si introducono alcuni concetti fondamentali oggetto dello studio della tesi. L’intenzione di questa sezione non è di presentare in modo esaustivo tutti gli argomenti sulla genetica e sulla biologia, per la quale si rimanda ai testi specifici, ma richiamare i fondamenti di biologia cellulare e genetica necessari per la comprensione degli argomenti di cui si discuterà in seguito. 1.1 DNA e Polimorfismi di un singolo nucleotide Il DNA è il veicolo per l’immagazzinamento e la trasmissione dell’informazione genetica, codificata nella sequenza lineare dei nucleotidi, utile alla sopravvivenza della cellula nell’organismo. Questa informazione è trasferita in maniera opportuna alla cellula, affinché questa possa essere operativa. La molecola deputata al trasferimento di questa informazione è l’RNA e il segmento codificante per un prodotto biologico operativo è il gene. Il gene rappresenta un carattere ereditario (mendeliano o monofattoriale) definito da un segmento di DNA. In ciascun cromosoma, i geni hanno un ordine preciso e ciascuno occupa una posizione specifica detta locus. Le forme alternative di un gene di dicono alleli; in un individuo, i due alleli occupano sui cromosomi omologhi lo stesso locus. È da sottolineare che i geni rappresentano la porzione codificante del genoma e corrispondono solo a circa il 2% dell’intera sequenza. Il resto è costituito da sequenze non codificanti (ripetizioni, sequenze introniche, regioni intrageniche,..) e per gran parte di esse non è ancora chiara la funzione. Il gene è l’unità fondamentale, fisica e funzionale, dell’informazione genetica e contiene le istruzioni per l’assemblaggio di proteine e RNA funzionali. L’espressione dell’informazione codificata in un gene avviene in due stadi: la trascrizione e la 1 traduzione. Nella trascrizione il filamento di DNA fa da stampo per produrre RNA messaggero (mRNA) che viene poi opportunamente modificato. L’mRNA servirà poi a sua volta da stampo per l’assemblaggio di proteine. La salvaguardia del materiale genetico richiede meccanismi estremamente precisi sia di duplicazione che di riparazione. Nonostante ciò, possono avvenire nel DNA di una cellula delle variazioni casuali della normale sequenza nucleotidica (mutazioni). Le malattie genetiche sono causate proprio da mutazioni del genoma. Si parla di polimorfismo quando una variazione genetica ha una frequenza maggiore dell’1% nella popolazione, dove con variazione si intende sostituzione, aggiunta o delezione di un nucleotide. Il fenomeno può avvenire lungo l’intera sequenza, quindi sia in corrispondenza di regioni codificanti che non codificanti. Il più semplice dei polimorfismi è il risultato della mutazione di una singola base in cui il nucleotide viene sostituito con un altro. I polimorfismi di un singolo nucleotide (Single Nucleotide Polimorphism, SNP, Figura 1.1.a, ) vengono definiti come differenze di una singola base in sedi specifiche. Gli SNP costituiscono fino al 90% delle differenze delle sequenze tra individui umani e nel genoma umano si verificano a diversi intervalli di basi (mediamente ogni 200 basi nucleotidiche). La presenza degli SNP, e degli alleli con cui si presentano, è verificata attraverso il sequenziamento di un campione di DNA. Lo specifico set di alleli osservato su un singolo cromosoma, o su parte di esso, è chiamato aplotipo (Figura 1.1.b ). Nuovi aplotipi sono formati da ulteriori mutazioni o da eventi di ricombinazione che hanno luogo nel momento in cui i cromosomi materno e paterno si scambiano rispettivamente i corrispondenti segmenti di DNA per dare un cromosoma che risulta un “mosaico” degli aplotipi parentali. È stato empiricamente osservato che il numero degli aplotipi nella popolazione è limitato, in quanto gli SNP e i loro alleli mostrano un forte livello di associazione (fenomeno conosciuto come linkage disequilibrium che sarà descritto in dettaglio in seguito). L’elevata frequenza di questi marker nel genoma li rende uno strumento ideale per la precisa mappatura genica di differenti malattie. Sebbene la maggior parte di essi non abbia effetti diretti sulla funzione cellulare (per esempio non cambiano la sequenza proteica, né alterano la trascrizione), molti sono strettamente associati, attraverso il linkage disequilibrium (paragrafo 1.2), con alleli o geni causa di malattie e spesso controllano il modo in cui il soggetto risponde a specifici farmaci. 2 Figura 1.1: SNPs e aplotipi. (a): Sono mostrati quattro segmenti di DNA da quattro versioni della stessa regione cromosomica di persone diverse. La maggior parte della sequenza è identica in questi cromosomi, ma si evidenziano tre basi in cui si verifica una variazione. Ogni SNP ha due possibili alleli; il primo allele mostrato nel pannello a presenta gli alleli C e T. (b): Aplotipi. Un aplotipo è composto da una particolare combinazione di alleli di SNPs vicini. Sono mostrati qui i genotipi osservati per 20 SNPs che si estendono per 6,000 basi di DNA. Vengono mostrate solo le basi variabili, inclusi i tre SNPs che sono mostrati nel pannello a.(c): sono evidenziati i tagSNPs che sono gli SNPs più rappresentativi dell’aplotipo, fortemente correlati con gli altri. 1.1.1 Frequenza degli SNPs e Minor Allele Frequency La frequenza di un allele a in una popolazione è definita come = Si propone di seguito un esempio di calcolo della frequenza di un allele. Si suppone di avere una popolazione costituita da 8 marker che presentano le seguenti coppie alleiche {AA, Aa, AA, aa, Aa, AA, AA, Aa}. Si vuole ora calcolare la frequenza dell’allele A; considerando che ognuno dei marker ha una coppia di alleli, quindi si ha un totale di 16 alleli nel locus. Risulta: 2+1+2+0+1+2+2+1 = = 0.6875 16 3 Sono verificate anche le seguenti relazioni: 1 = + = 2 1 = + = = 1 − 2 All’intero di una popolazione è possibile determinare una minor frequenza allelica (minor allele frequency, MAF), definita come il rapporto tra la frequenza della variante più rara e quella più comune di un determinato SNP. Lo studio della MAF consente di distinguere gli SNP monomorfici (MAF = >K − >C è quindi la misura della deviazione della sequenza osservata dallo stato di linkage equilibrium atteso. Per permettere confronti tra diversi set di loci si è proposta la versione scalata di ∆H compresa tra [0,1) indicata con ε: ∆> >C M= =1− >K >K ε è la misura per il LD denominata Nomalized Entropy Difference. Per definizione, la misura ε permette di considerare un numero illimitato di loci contemporaneamente. Comunque, nel caso pratico, il numero è limitato dalle dimensioni del campione. Alcuni aplotipi che possono essere presenti con frequenze rilevanti in alcune popolazioni, potrebbero non essere presenti nei campioni osservati, dovuto anche semplicemente al 18 fatto che il campione è troppo piccolo per poterli osservare. In questo caso, la ε può portare a una sovrastima del LD. La misura è sensibile sia al numero di aplotipi osservati che alle loro frequenze. ε è pari a 0 se e solo se una sequenza è in stato di equilibrium, aumenta al diminuire del numero di aplotipi presenti e aumenta all’aumentare della deviazione delle frequenze dalla condizione di equilibrium. La misura distingue anche tra vari gradi di LD oltre che l’assenza di più di un aplotipo. Sia le misure Pair-Wise che Multilocus, sono influenzate fortemente nella loro accuratezza della dimensione del campione su cui viene fatto il calcolo di LD. In particolare, campioni più grandi possono minimizzare l’errore di campionamento e produrre una valutazione più adeguata del LD. Inoltre i soggetti appartenenti a gruppi famigliari possiedono più informazione rispetto a campioni di individui scorrelati tra loro. Segue che la valutazione di LD su nuclei famigliari è più accurata. Risulta evidente, da quanto detto, che i pattern di LD osservati nella popolazione sono il risultato di una complessa interazione tra i fattori genetici e la storia demografica della popolazione stessa. In particolare, la ricombinazione gioca un ruolo chiave nel formare e modellare i pattern di LD in una popolazione. Quando avviene un evento di ricombinazione tra due loci, questo tende a ridurre la dipendenza tra gli alleli che quegli stessi loci portano e quindi riduce il LD. Anche se gli eventi di ricombinazione in una singola meiosi sono relativamente rari, se si considerano piccole regioni, il numero totale di meiosi che hanno luogo ad ogni generazione ha un sostanziale effetto cumulativo sui pattern di LD. A causa della complessità del fenomeno, molti modelli per l’interpretazione e l’analisi del LD non tengono conto di questo aspetto. 2.2 Definizione degli aplotipi mediante misure di LD Gli aplotipi, o blocchi di aplotipi, rappresentano regioni ereditate senza sostanziale ricombinazione dagli antenati nella popolazione moderna. La storia della ricombinazione tra le coppie di SNPs può essere stimata con l’uso delle misure definite precedentemente. Se si utilizza la misura D’, gli aplotipi vengono quindi definiti come regioni in cui la misura standard del disequilibrium D’ è considerata pari a 1 (o comunque molto vicina al valore unitario, facendo opportunamente riferimento a degli intervalli di confidenza), in assenza di ricombinazione, per tutte (o quasi) le coppie di marker nella regione. Poiché i valori di D’ tendono ad essere sovrastimati, nel caso in 19 cui si abbiano a disposizione pochi campioni o si lavori con alleli rari, si preferisce fare riferimento a intervalli di confidenza piuttosto che alle stime puntuali della misura. A separare gli aplotipi, ossia regioni in cui la ricombinazione è al livelli considerevolmente ridotti, sono le variazioni locali, posizionate lungo il cromosoma, denominate punti caldi, o hotspots, dove al contrario si registra un tasso di ricombinazione elevato e una caduta significativa del valore di LD. Si possono definire hotspots quei segmenti di lunghezza mediamente inferiore alle 10 kb, dove la misura di r2 tra due marker vicini non eccede mai il valore di 0.10. I punti caldi, in sostanza, rappresentano interruzioni del LD, costituiscono meno del 10% delle sequenze, e si è stimato che in queste brevi regioni avvenga circa il 50% degli eventi di ricombinazione. Queste variazioni nel tasso di ricombinazione spiegano almeno in parte le recenti osservazioni del grado di eterogeneità nei valori di LD tra i marker SNPs lungo il genoma. Considerando nuovamente la popolazione CEU rappresentata prima (in Figura 2.1, sul LD), si propone ora la mappa dei punti caldi lungo l’intero genoma in Figura 2.4. Gli hotspots sono individuati mediante la colorazione rossa. Figura 2.4: riepilogo per ogni cromosoma degli intervalli tra i marker (in rosso) dove il LD decade molto rapidamente. Intervalli in cui il LD decade in modo meno pronunciato sono evidenziate in azzurro. In ascissa è riportata la distanza in Mb, in ordinata il numero identificativo del cromosoma [Albert VS, Sequence features in regions of weak and strong linkage disequilibrium, Genome Res 2005]. 20 2.3 Algoritmi di LD e inferenza di aplotipi Gli algoritmi per il calcolo del linkage sono derivati dall’iniziale problema del mappaggio genico in base alla frazione di ricombinazione, così da risalire alla posizione del locus malattia nel genoma rispetto al marcatore. In molti casi gli aplotipi non sono “letti” direttamente nel momento in cui il DNA umano viene sequenziato (fatta eccezione di alcuni casi come il sequenziamento del cromosoma Y o tecniche di conversione diploide-aploide), ma vengono dedotti dai dati di genotipo. Un problema fondamentale con la ricostruzione di mappe di linkage nell’uomo è che alcuni dati importanti spesso mancano. In questo modo non è possibile semplicemente contare le ricombinazioni negli incroci, proprio perché la mancanza di tutte le informazioni non permette di capire senza ambiguità dove queste sono avvenute. Nonostante l’indiscussa necessità di capire i pattern di LD nel genoma, soprattutto a causa dell’impatto che questo aspetto ha nel design e nell’analisi del mappaggio dei geni “malattia” nell’uomo, i metodi più comunemente usati per l’interpretazione e l’analisi del LD soffrono almeno di una delle seguenti limitazioni: 1. Sono basati su una valutazione delle misure di LD definite solo per coppie di loci, piuttosto che considerare tutti i loci contemporaneamente; 2. Assumono una struttura a blocchi per schematizzare i pattern di LD, il che potrebbe non essere appropriato per tutti i loci; 3. Non mettono direttamente in relazione i pattern di LD con i meccanismi biologici di interesse, come ad esempio il tasso di ricombinazione. In letteratura si possono trovare un’ampia gamma di algoritmi che permettono di risolvere il problema di inferenza di aplotipi. È importante ricordare che questi algoritmi sono pensati per applicazioni su dati di genoma diploide, quindi con patrimonio genetico che presenta due copie di ogni cromosoma. Mentre gli aplotipi rappresentano le informazioni degli alleli di SNPs su un cromosoma, i genotipi rappresentano le informazioni combinate di alleli di SNPs su due cromosomi. Se ho n set di genotipi, dove ognuno riporta l’informazione di m SNPs, il numero di aplotipi (massimo) sarà 2m. È necessario sottolineare che il problema di inferenza di aplotipi non è di immediata e diretta soluzione, dovuto all’ambiguità di risoluzione. L’ambiguità nasce dalla presenza di SNPs eterozigoti: in particolare, quando si è in presenza di c (dove c > 1) SNPs eterozigoti nel genotipo, ci sono 2c-1 coppie di aplotipi che possono risolvere il genotipo. 21 Quindi, il genotipo non può essere univocamente risolto senza l’aggiunta di vincoli o considerazioni di tipo biologico. Un tipo di approccio, uno tra i primi implementato, fa riferimento a modelli di parsimonia. Questi metodi assumono che una data popolazione target condivida un numero relativamente ridotto di aplotipi in comune dovuto al linkage disequilibrium. Gli algoritmi si basano sulla costruzione incrementale della soluzione mediante una applicazione iterata di una regola di inferenza, detta regola di Clark (dal nome del ricercatore che per primo ha introdotto il principio nel 1990). L’algoritmo identifica inizialmente i genotipi che contengono solo alleli omozigoti o al più un unico allele eterozigote. Questi genotipi possono essere risolti in maniera univoca e la corrispondente coppia di aplotipi viene memorizzata in un set di aplotipi già identificati, denotato con I. Per i rimanenti genotipi “ambigui”, si esamina l’insieme I per vedere se in esso sono contenuti aplotipi che siano compatibili con il dato genotipo. Se viene trovato un aplotipo tale da soddisfare questa condizione, il rispettivo genotipo viene etichettato come “risolto”. Il procedimento è iterato fintanto che tutti i genotipi sono stati risolti o non sono trovati nuovi aplotipi. In questi modelli, l’algoritmo ricerca per ogni passo una soluzione localmente ottima (algoritmo di greedy), cercando una soluzione locale che incrementi il meno possibile il numero di aploitipi, senza vincoli relativi alla soluzione globale. L’algoritmo di Clark è semplice e intuitivo, tuttavia presenta numerosi svantaggi: molti dei genotipi “ambigui” possono rimanere non risolti e un ordine differente di iterazione può portare alla costruzione di aplotipi diversi. La mancanza di un modello di soluzione globale implica che la soluzione trovata può essere diversa, ad ogni applicazione dell’algoritmo, a seconda delle scelte effettuate ad ogni passo, quindi non tutte le sequenze di applicazione della regola di inferenza consentono di raggiungere una soluzione al problema. Un’altra formulazione del problema si basa sul principio di pura parsimonia, proposto da Gusfield. In questo caso si ricerca l’insieme di aplotipi di minima cardinalità che risolve l’insieme dei genotipi dati. Tutti i metodi basati sui modelli di parsimonia assumono che il numero di aplotipi distinti in una popolazione sia più piccolo del numero possibile di aplotipi distinti che si avrebbero in condizione di linkage equilibrium (cioè in totale assenza di linkage disequilibrium). Perciò, quando un dataset non soddisfa questa condizione, le performance di questi modelli diventano insoddisfacenti. 22 I modelli coalescenti, invece, si basano su assunzioni biologiche relative alla storia evolutiva delle mutazioni. In particolare, secondo questi modelli, le mutazioni (cioè le transizioni da un allele all’altro in alcuni polimorfismi) avvengono una sola volta e, per questo, sono condivise da un insieme di individui con il medesimo antenato. Ne segue che un cromosoma sprovvisto di tale mutazione, non può essere discendente del medesimo antenato che invece presenta quella mutazione. Inoltre è assunta l’assenza di ricombinazione; segue da ciò che una sequenza target non viene considerata come il risultato di un evento di ricombinazione tra i due cromosomi parentali, ma è come se venisse ereditata da un singolo progenitore. Un approccio di questo tipo assume che gli aplotipi di una popolazione evolvano secondo un modello genetico identificato da un grafo ad albero che descrive la storia evolutiva di un set di sequenze di DNA. Il termine con cui si identifica l’albero è filogenesi perfetta. Se si hanno 2n aplotipi, dove ogni aplotipo è costituito da m SNPs, l’albero sarà composto da 2n foglie, esattamente una per aplotipo, e ognuno degli m SNP andrà a costituire gli m archi. In Figura x.5 è riportato un esempio di filogenesi perfetta per un set di 4 aplotipi. Figura x.5: esempio di filogenesi perfetta. In alto è proposto lo spettro dei 4 aplotipi composti dai 3 SNPs. Mentre la figura sottostante è la rappresentazione mediante grafo ad albero della filogenesi perfetta. Anche se le perfomance dei metodi basati sulla filogenesi perfetta sono state migliorate, questi soffrono ancora della stretta conformità al modello coalescente: è possibile che non esista filogenesi perfetta per un dato genotipo. Nel caso reale, spesso i dati non soddisfano i requisiti del modello coalescente e violano le assunzioni biologiche alla 23 base del modello stesso. Le principali ragioni possono essere dovute a errori nel sequenziamento o assenza dell’informazione relativa alla ricombinazione. Per mantenere la filogenesi perfetta si è proposto di eliminare un numero di genotipi dai dati originali in modo che i restanti potessero soddisfare i requisiti di filogenesi perfetta o assegnare un valore arbitrario agli alleli mancanti in modo che i genotipi così ricostruiti potessero ancora soddisfare la filogenesi perfetta. Le soluzioni proposte ancora non hanno dato risultati adeguati e in molti casi manca una soluzione euristica. Un approccio più realistico è offerto dai metodi di filogenesi imperfetta: questi assumono che la maggior parte dei genotipi (non tutti) soddisfino il modello della filogenesi perfetta. Quindi, considerano un modello meno “rigido” del precedente che permetta di inserire un certo numero di mutazioni ricorrenti e di ricombinazioni. Tra più soluzioni candidate che soddisfano questo modelli, quella che soddisfa la maximum-likelihood, dato il genotipo, è scelta come soluzione. Tuttavia, la gestione di un numero esponenziale di soluzioni candidato possibili rimane ancora un problema non risolto. I metodi basati sui modelli di parsimonia e sui modelli coalescenti propongono un approccio diretto nella risoluzione del genotipo con ogni coppia di aplotipo. Viceversa, i metodi basati su modello statistico si basano su un approccio più indiretto. Gli algoritmi proposti nei modelli precedenti richiedono un insieme iniziale di aplotipi risolti e la loro soluzione dipende fortemente dall’ordine in cui vengono risolti. L’approccio proposto dai modelli statistici diminuisce l’importanza di tale dipendenza, rendendo la procedura più affidabile. L’idea principale è che gli aplotipi hanno una distribuzione di probabilità non nota nella popolazione in esame e che i genotipi osservati di ogni individuo sono semplicemente combinazione di due aplotipi presi a caso dalla popolazione. L’obiettivo dell’approccio statistico è quindi di stimare la frequenza degli aplotipi in modo che sia possibile calcolare facilmente gli aplotipi di ogni individuo, basandosi sulla distribuzione di probabilità e su considerazioni di tipo biologico. I metodi di risoluzione sono di due tipologie: la Maximum Likelihood Iinference e il Bayesian Frequencies Haplotype Inference problem. Gli elementi principali di un modello di tipo statistico per generare dati relativi agli aplotipi e ai genotipi sono: G = (g1, g2, …, gm) è l’insieme dei genotipi osservati; 24 H = (H1, H2, …, Hm) sono le corrispondenti coppie di aplotipi incognite, dove Hi = (hi1, hi2) per i = 1,…,m; ϴ = (θ1, …, θv) indica il valore delle v frequenze non note degli aplotipi (dove v è la cardinalità dell’insieme di tutti gli aplotipi); Comunemente, la stima delle frequenze aplotipiche di una popolazione è fatta partendo da sequenze geniche di individui non correlati tra di loro, campionati casualmente, sotto l’assunzione dell’equilibrio Hardy-Weinberg. La formulazione Maximum Likelihood Inference (o formulazione ML) ha come obiettivo quello di massimizzare la distribuzione degli aplotipi nella popolazione, massimizzando la verosimiglianza (likelihood) dei dati di genotipo. In input viene dato l’insieme dei G genotipi: a differenza delle frequenze del genotipo, che possono essere direttamente calcolate dal dataset, le frequenze aplotipiche sono non note e vanno stimate. In uscita si ha l’insieme delle frequenze degli aplotipi {h1, …, hv} (dove v è il numero di tutti i possibili aplotipi) che massimizza la funzione di somiglianza (likelihood) per l’insieme di genotipi osservato. Usando le frequenze stimate, ogni genotipo può essere risolto dalla coppia di aplotipi con la frequenza massima tra tutte le coppie compatibili con il dato genotipo. L’algoritmo implementato più diffusamente e proposto per la risoluzione di un problema di Maximum Likelihood è l’Expectation Maximization (E-M) algorithm di Excoffier e Slatkin (1995), e viene proposto per stimare gli aplotipi e le frequenze che massimizzano la funzione di somiglianza dell’insieme dei genotipi. Sia ϴHt l’insieme delle frequenze degli aplotipi e Gt l’insieme delle probabilità di tutti i genotipi al passo t. L’algoritmo EM assegna un valore iniziale alle frequenze aplotipiche ϴH0 (un insieme iniziale possibile di frequenze è quello che corrisponde all’assunzione che tutti i possibili aplotipi siano equiprobabili). Basandosi su ϴH0 , può essere facilmente calcolato il valore atteso di un genotipo di G e può quindi essere calcolato il valore di G1 (Expectation step) I valori attesi dei genotipi in G1 vengono usati successivamente per stimare nuovamente le frequenze degli aplotipi, attraverso lo step di massimizzazione, il che porta a calcolare ϴH1. L’algoritmo consiste nell’iterare i due passi (l’expectation e il maximization) fino alla convergenza, cioè finché la differenza tra ϴHt e ϴHt+1 non sia più piccola di un valore predefinito. Ad ogni iterazione, il valore di ϴHt viene migliorato massimizzando la funzione di somiglianza dell’insieme G. Diversi autori riportano il passo di Expectation dell’algoritmo mediante la seguente formulazione: 25 ∑UVW6 OP (>)P, P! 4NOP (>6 )|R6 )S = ∑UVW6 P, P! Dove δh(H) = 0,1, o 2 è il “dosaggio aplotipico” ossia il conteggio del numero di copie di h contenute nel vero (ma in generale non nota) coppia di aplotipi H portata dall’individuo i-esimo. La stima di δh(H) è calcolata condizionata ai dati di genotipo G per ogni soggetto e considerando il set delle frequenze p come se fossero note. La sommatoria ∑UVW6 indica la sommatoria sulle coppie di aplotipi ordinate H=(h1, h2), con frequenze ph1 e ph2 rispettivamente, che sono compatibili con i dati di genotipo osservato. Poiché si assume l’equilibrio Hardy-Weinberg, la coppia H ha probabilità pari al prodotto ph1 ph2. Questo algoritmo risulta di utile applicazione nei modelli con variabili nascoste: tipici esempi di variabili nascoste sono i dati mancanti o non osservabili. Si procede, infatti, andando ad assegnare un valore atteso alle variabili nascoste come se queste fossero note. La complessità temporale dell’algoritmo EM è O(m2k) dove m è il numero di genotipi e k è il numero massimo di SNPs eterozigoti nei genotipi. Il limite maggiore dell’algoritmo risiede nell’incremento esponenziale del numero di aplotipi all’aumentare del numero di SNPs eterozigoti. La soluzione più comune al problema è la divisione dell’intero set di SNPs in sottoinsiemi più piccoli e contigui, definiti anche pseudo-blocchi, procedendo poi alla combinazione degli aplotipi selezionati dai singoli sottoinsiemi attraverso un approccio bottom-up. Poiché l’algoritmo si basa sul’assunzione dell’equilibrio H-W, e dato che questa condizione è quasi sempre soddisfatta se si considerano grandi gruppi di individui, allora l’algoritmo può essere applicato con risultati soddisfacenti su campioni di dati molto grandi. Studi di simulazione hanno comunque dimostrato la validità dell’algoritmo anche quando l’assunzione dell’HWE è violata. Sono state inoltre valutate le performance dell’algoritmo quando si è in presenza di errori di sequenziamento e/o dati mancanti: se si è in condizioni di moderato o forte LD tra gli SNPs, l’assenza fino al 30% di dati non sembra influire sull’accuratezza del risultato. Tuttavia, se si è in condizioni di basso LD, gli errori di sequenziamento influiscono pesantemente sul risultato diminuendone l’accuratezza, pertanto in tale situazione è preferibile considerare i dati incerti come non noti. Come già anticipato, l’algoritmo EM è diffusamente impiegato e si può trovare implementato anche nel software Haploview (2003-2006 Broad Institute of MIT and Harvard). 26 Una seconda formulazione è la Bayesian Frequencies Haplotype Inference problem: come la Maximum Likelihood, si propende anche in questo caso per un approccio di tipo statistico. Mentre nel Maximum Likelihood si andava a trovare un set di parametri del modello che andassero a massimizzare la probabilità dei dati di genotipo G dato il modello, l’approccio Bayesiano va a calcolare la distribuzione a posteriori dei parametri del modello dati i dati di genotipo G, ossia P(ϴ|G). In input sono passati un insieme G di genotipi e una distribuzione a priori delle frequenze dei genotipi; in uscita è riportata la distribuzione a posteriori delle frequenze degli aplotipi dato G. Rispetto ai modelli ML, questi richiedono più tempo computazionale ed è più difficile raggiungere la convergenza del risultato. Gli approcci statistici sono i metodi più popolari e più largamente implementati. L’accuratezza di questi metodi è infatti in qualche modo migliore di quella risultante dall’applicazione delle altre soluzioni. Inoltre, i metodi di parsimonia e i modelli coalescenti spesso presentano soluzioni multiple, rendendo difficile il confronto delle loro performance con altri metodi. Per ultimo, i metodi statistici risultano sicuramente più robusti perché applicabili anche nel caso di genotipi mancanti o ambigui. Questa condizione non è infrequente, in quanto dovuta ad errori di sequenziamento o alleli mancanti. La mancanza dell’informazione sull’allele aumenta notevolmente la complessità computazionale dell’algoritmo, mentre l’errore di sequenziamento rende più difficile la risoluzione del problema, dal momento che non si conosce quale allele sia sbagliato. È da sottolineare inoltre che, nonostante i risultati di questi approcci siano molto promettenti, possono comunque cadere in difetto quando i livelli di LD decrescono. Questa scarsa accuratezza si presenta soprattutto quando consideriamo un numero consistente di SNPs, dal momento che il LD tende a decrescere con l’aumentare della distanza tra SNPs. Scarsa accuratezza si registra anche in presenza di rari aplotipi, in quanto il razionale alla base di molti algoritmi è la condivisione del maggior numero di aplotipi possibile nella popolazione. In conclusione, alla luce delle questioni illustrate, la ricerca sugli algoritmi di ricostruzione di LD e aplotipi dovrebbe concentrarsi sul miglioramento delle performance di questi in dataset con basso LD, in cui sono presenti errori di sequenziamento, alleli mancanti e aplotipi rari. 27 2.4 Tag SNP In regioni ad alto LD, si può selezionare un ridotto set di SNP per identificare efficientemente gli aplotipi comuni. Un ridotto sottoinsieme di SNP è preferibile sia per ridurre le esigenze di sequenziamento del genoma, sia per eliminare la ridondanza di informazione data dal considerare gli SNPs nella loro totalità. Quindi, per ovviare alla grande quantità di dati che è necessario per risolvere problemi legati al DNA, molti scienziati si sono dedicati alla ricerca di sottoinsiemi minimi di SNP che permettano comunque di rappresentare le associazioni tra malattie e geni. Questo processo viene chiamato identificazione di tagSNP e in generale gli SNP selezionati vengono indicati come tagSNP e quelli non selezionati come taggedSNP. Per risolvere un problema di selezione di TagSNP è necessario trovare un sottoinsieme ottimale di SNPs, di dimensione minore rispetto l’originale, valutando come questo rappresenti bene i dati di genotipo rispetto a tutti gli altri possibili sottoinsiemi ottenibili. La motivazione della selezione di TagSNP ha origine, come già accennato sopra, nel linkage disequilibrium, concetto che peraltro è alla base dei test di associazione gene-malattia. Quando è presente un alto LD tra due SNPs, l’informazione portata dai loro alleli è pressoché identica. Quindi, possiamo selezionare uno tra gli SNPs ridondanti in modo che, anche un sottoinsieme dell’originale gruppo di SNPs, mantenga la stessa informazione. Quale sia la migliore strategia che permetta la selezione del sottoinsieme ottimale, rimane ancora un problema aperto. Vengono di seguito proposte vari approcci possibili al problema finora sviluppati. Come già osservato nel capitolo introduttivo, la struttura a blocchi del genoma umano dimostra che il genoma può essere partizionato in blocchi discreti tali che, all’interno di ciascuno, la maggior parte della popolazione condivide un piccolo numero di aplotipi comuni (all’incirca 3.5). Basandosi su questa assunzione, i primi algoritmi sviluppati miravano a trovare il sottoinsieme di SNPs che meglio catturava la maggior parte della limitata diversità tra aplotipi dai dati originali. Questi metodi sono definiti Haplotype Diversity-based Methods. Per definire e quantificare la diversità tra gli aplotipi, sono state proposte diverse misure. Alcune usano, come misura di diversità, il numero di aplotipi che sono distinguibili in modo univoco dal sottoinsieme T’ di TagSNPs. Si sceglie quindi il sottoinsieme con la misura di diversità più alta. Un’altra soluzione è definire la diversità non catturata dal sottoinsieme T’ di tagSNPs come il numero di alleli diversi tra ogni coppia di aplotipi nello stesso gruppo basato su T’. Se il 28 sottoinsieme T’ partiziona correttamente tutti i distinti aplotipi in gruppi diversi, allora la diversità aplotipica residua sarà nulla. Si andrà a selezionare in questo caso il sottoinsieme T’ che mostra diversità residua più piccola o nulla. Un’altra misura popolare della diversità aplotipica è basata sull’entropia di Shannon. Sia n’ il numero di aplotipi distinti nel dataset D di aplotipi, e pi siano le frequenze relative dell’i-esimo aplotipo distinto. La diversità di D può essere calcolata come entropia H: : >(#) = − ? 6 log ! 6 6D, Come negli altri metodi, gli aplotipi sono partizionati in gruppi in modo che quelli nello stesso gruppo condividano gli stessi alleli. L’entropia del dataset D è misurata basandosi su questa partizione. Gli aplotipi che sono posizionati nello stesso gruppo sono considerati identici. Maggiore è il numero dei sottoinsiemi candidato T’, maggiore sarà l’entropia del dataset basato su quella partizione. Il sottoinsieme con la misura di entropia più alta sarà selezionato come soluzione. I metodi finora introdotti esaminano esaustivamente tutti i sottoinsiemi del set originale di SNP, limitando la loro applicazione solo a piccoli insiemi di SNPs. Per ovviare a questa limitazione sono state proposte diverse soluzioni di tipo euristico o algoritmi greedy. I metodi basati sulla misura della diversità sono intuitivi e diretti. Tuttavia, per assicurarsi che la diversità tra aplotipi sia effettivamente limitata e per velocizzare la performance dell’algoritmo, si esegue una divisione in blocchi adiacenti della regione in esame, e la selezione dei TagSNP viene fatta blocco per blocco. Una limitazione a questo modo di operare in blocchi sta nella possibilità che l’unione dei set di TagSNP definiti sui blocchi possa non essere il sottoinsieme ottimale di TagSNP per l’intera regione. Per di più, tra i blocchi possono esistere, con alta probabilità, regioni a basso LD che vanno ad aumentare il numero di aplotipi diversi rendendo inapplicabili tali metodi. Risulta evidente quindi che, con questo approccio, il particolare partizionamento della regione target influenza pesantemente la selezione dei TagSNP. E’ possibile effettuare la scelta dei TagSNP affidandosi anche all’algoritmo EM descritto nei paragrafi precedenti in merito all’inferenza degli aplotipi. Infatti la stessa procedura può essere applicata alle situazioni in cui dobbiamo scegliere il particolare set di TagSNP. Si utilizza nuovamente lo stimatore, E{δh(H)|Gk}, del numero di copie δh(H) con l’osservato genotipo G, assumendo di nuovo che le frequenze aplotipiche 29 siano note e valido l’equilibrio Hardy-Weinberg. Si procede quindi con il calcolo della correlazione R2h tra E{δh(H)|Gk} e δh(H) come : X Y4NOP (>)|RSZ X Y4NOP (>)|RSZ 9P! = = X OP (>) 2P (1 − P ) Dove ph è la frequenza dell’alotipo h, e la varianza dell’aspettazione a nominatore è calcolata come sommatoria su tutti i possibili valori del genotipo G come: X Y4NOI (>)|RSZ = ? 4NOI (>)|RS! (R) − 4I! W Per ogni sottoinsieme di SNPs si può formalmente calcolare R2h usando la formula riportata. Per il ridotto set di SNP ci saranno più coppie aplotipiche che saranno compatibili per il genotipo basato solamente sul ridotto set di SNP piuttosto che sul genotipo basato sul data set con tutti gli SNPs. Ne consegue una varianza minore e quindi un minor valore per R2h. Il calcolo dei TagSNP prevede innanzitutto l’applicazione dell’algoritmo EM per identificare gli aplotipi comuni all’interno di regioni ad alto linkage disequilibrium. Per un blocco comprendente n SNPs ad alto LD, si definisce il miglior set di m TagSNP (m < n) come quegli SNPs che massimizzano il minimo valore di R2k calcolato per ogni aplotipo comune. Il calcolo di R2k per ogni aplotipo richiede la generazione del’intero set di copie di aplotipi, H, per un dato set di frequenze aplotipiche non nulle e un valore di δh(H) compatibile con ognuno dei possibili SNP sequenziati. Questo deve quindi essere fatto per ogni aplotipo comune h, e per ogni set candidato di TagSNP. In dipendenza del numero di SNPs, il procedimento può risultare piuttosto oneroso in molti casi. Per ottimizzare la scelta, invece di optare per una ricerca esaustiva di tutte le ! ( − )! ! scelte dei m TagSNPs è stato implementato un metodo di inclusione passo per passo: si sceglie come TagSNP candidato, da un set di k SNP, quello che contribuisce al maggior incremento del valore di R2h calcolato sui rimanenti k-1 SNPs. Il processo viene ripetuto per vedere se R2h può essere aumentato ulteriormente per sostituzione con ognuno degli altri k-1 SNP con ogni altri SNP non selezionato come Tag. 30 Una seconda categoria di metodi proposti sono i Pairwise Association-based Methods. Questo approccio si basa sull’idea che un set di TagSNP dovrebbe essere il più piccolo sottoinsieme di SNPs disponibili in grado di predire il locus malattia di un aplotipo. Il locus malattia è l’obiettivo degli studi di associazione, quindi generalmente non è noto a priori e si deve procedere con una stima. Il set di tagSNP è selezionato in modo tale che tutti gli SNPs di un aplotipo siano altamente associati con uno dei tagSNP selezionati. In questo modo, anche se uno SNP che svolge un ruolo rilevante in una malattia non viene selezionato come tag, la sua associazione alla malattia può essere indirettamente dedotta dai tagSNP, in quanto sono scelti secondo un criterio di forte associazione con esso. Nella maggior parte degli studi, il LD è utilizzato come misura per la stima dell’associazione pairwise. Sono stati proposti due algoritmi per la soluzione del problema. Il primo ricorre a una partizione del set originale di SNPs mediante clustering gerarchico, dove gli SNPs all’interno dello stesso cluster hanno un valore di LD con almeno uno degli altri SNPs maggiore di un livello scelto a priori. Inizialmente, ogni SNP costituisce un cluster a sé stante. Le fusioni tra i cluster Ci e Cj si basano sulla seguente definizione di distanza: #5]6 , ]7 8 = min (#Je ()) ∀_∈(ab ∪ad ) dove Dmax(s) è la massima distanza, definita dal valore di LD (ad esempio r2), tra lo SNP s e tutti gli altri SNPs nei due cluster. I due cluster più vicini, basati sulla distanza D(Ci,Cj), sono quindi uniti iterativamente. L’unione dei cluster si ferma quando la più piccola distanza tra due cluster è più grande di una certa soglia. Lo SNP s che definisce la distanza di ogni cluster unito, è scelto come rappresentativo del cluster. Generalmente, viene scelto uno SNP da ogni cluster basandosi su considerazioni di tipo pratico come la facilità di sequenziamento, l’importanza della collocazione del locus o significatività portata dalla mutazione dello SNP. Il secondo algoritmo proposto è un algoritmo di greedy: questo innanzitutto esamina tutte le relazioni di LD tra tutte le coppie di SNPs, e per ogni SNP conta il numero di SNPs con i quali ha un valore di LD maggiore di una certa soglia. Lo SNP con il maggior numero di conteggi è clusterizzato insieme agli SNPs associati ad esso e diventa il tagSNP del cluster. Questa procedura è iterata con i rimanenti SNPs fintanto che tutti gli SNPs non sono stati clusterizzati. Gli SNP che invece presentano livello di LD minore della soglia fissata, sono clusterizzati singolarmente. Alternativamente, è stata proposta per i metodi di Pairwise Association 31 la scelta di tagSNP tali che il loro valore di LD sia come prima maggiore di un certo livello fissato, ma in questo caso il valore di LD deve essere soddisfatto non solo per uno, ma per tutti gli SNP nel cluster. Anche per questo criterio sono applicabili i due algoritmi appena descritti. Tutti i metodi di associazione pairwise hanno complessità O(cnm2), dove il numero di cluster è c, il numero di aplotipi è n, e il numero di SNPs è m. Quindi, in generale, sono veloci dal punto di vista computazionale rispetto ai metodi basati sulla definizione di distanza, e non richiedono una partizione in blocchi della regione target. Il principale svantaggio è rappresentato dal fatto che non riescono a identificare dipendenze multiple tra più di due SNPs; tendono inoltre a selezionare più tagSNP rispetto agli altri metodi. Una possibile soluzione agli aspetti negativi presentati per i metodi Pairwise è rappresentata dai metodi Tagged SNP Prediction-based Methods. Questi metodi considerano la selezione di TagSNP come un problema di ricostruzione dell’aplotipo originale usando solo un sottoinsieme di SNPs. Il loro obiettivo è quindi selezionare un set di SNPs che possano predire gli SNPs non selezionati (i tagged) con il minimo errore. È stato inizialmente proposto di selezionare i tagSNP sulla base della loro accuratezza nel predire i tagged. Sono state quindi definite misure di informatività per trovare il sottoinsieme ottimale di SNP, utilizzate poi in programmazione dinamica. A differenza dei metodi pairwise, questo approccio permette di identificare associazioni multiple tra SNP e il numero di TagSNP è generalmente minore rispetto a quelli selezionati dai precedenti. Inoltre, il ricorso a tecniche di programmazione dinamica garantisce di trovare una soluzione che sia un ottimo globale. Tuttavia, i maggiori limiti di queste tecniche sono il tempo computazionale di ordine esponenziale e il confinamento che viene imposto alla regione di ricerca dei tagSNP nell’intorno di un locus fisico. L’ultima categoria è rappresentata dai metodi Phenotype Association-based Methods. Questo approccio assume la disponibilità dell’informazione del fenotipo e, sulla base di questa, cercano di trovare un set di SNPs che possa distinguere gli individui portatori di malattia (i casi) dagli individui sani (controlli). Sotto questa prospettiva, la selezione dei tagSNP è una sorta di feature selection, che mira a selezionare un set di features che distingue tra due classi (casi/controlli) con un piccolo errore. Uno dei classificatori più implementati è il classificatore Bayesiano. Assume che l’allele di uno SNP è indipendente da quello degli altri fenotipi dati, e classifica ogni aplotipo come caso o 32 controllo basandosi sulla sua probabilità di appartenere a una di queste classi. Il sottoinsieme T’ di SNPs con la miglior accuratezza di classificazione è selezionato come set di tagSNP. La limitazione più forte sta nell’assunzione di indipendenza tra SNP alla base del classificatore. Nel caso reale, infatti, tra gli SNPs esiste un’associazione non casuale, data dal linkage disequilibrium. È stato successivamente proposto un altro metodo che non solo classifica correttamente i dati, ma garantisce anche che le sue performance siano statisticamente significative. Questo algoritmo è basato su una tecnica di bootstrap nel quale si eseguono circa 1000 permutazioni su n coppie di aplotipo-fenotipo che costituiscono il dataset originale. Questi metodi sono direttamente correlati all’obiettivo principale dell’analisi di aplotipi e dei testi di associazione gene-malattia. la principale limitazione di questo approccio consiste nel la necessità dell’informazione sull’aplotipo, che potrebbe non essere disponibile in anticipo. Inoltre, generalmente, il numero di aplotipi usati per la selezione di tagSNP è relativamente piccolo, quindi i tag selezionati che classificano un campione ridotto molto bene, potrebbero non rappresentare altrettanto bene un campione più grande. Questo ha conseguenze che possono condizionare i test di associazione gene-malattia. La praticità della selezione dei TagSNP è stata empiricamente dimostrata in numerosi studi di simulazione. I risultati dimostrano che lavorare con i tagSNP porta a una perdita minima della potenza dei test di associazione (dove con potenza dei test di associazione intendiamo a probabilità che il test rifiuti l’ipotesi nulla). Tuttavia, rimangono numerosi problemi aperti: 1. È necessario assicurarsi della qualità dei database di SNP e la loro applicabilità a una popolazione più ampia di quella campionata. 2. La maggior parte degli algoritmi di selezione dei tagSNP si concentrano su aplotipi comuni o SNPs comuni piuttosto che quelli rari. Le variazioni comuni sono sicuramente di grande interesse in quanto molti disturbi e malattie comuni possono essere spiegate proprio da variazioni comuni di DNA piuttosto che da quelle rare. Tuttavia, è ancora aperta la questione su quali variazioni, se quelle comuni o quelle rare, influenzino la predisposizione a malattie comuni e complesse. 3. Diversi algoritmi richiedono dati di aplotipo piuttosto che di genotipo. Quando sono disponibili solo dati di genotipo, deve essere calcolato l’aplotipo da questi dati e lavorare su questa informazione. Si sa però che molti algoritmi sviluppati 33 per il calcolo degli aplotipi possono produrre più di una soluzione e dare quindi un certo grado di incertezza. Finora, nessun metodo di selezione di TagSNP considera l’incertezza derivante dalla ricostruzione di aplotipi. 4. Tutti gli algoritmi presentati assumono che il set di tagSNP selezionato da un campione possa dare buoni risultati anche per un altro campione proveniente dalla stessa popolazione. Per assicurarsi che le performance possano essere generalizzate, si dovrebbero campionare un numero sufficiente di individui ed evitare situazioni di overfitting. 5. È importante definire i limiti o confini dei blocchi di LD all’interno dei quali concentrare lo spazio di lavoro per la selezione dei tagSNP. Sebbene l’idea di trovare TagSNP non dipenda in sé dall’esistenza o dalla qualità dei blocchi di LD, il pattern di LD influenza la proprietà dei TagSNP, incluso il vantaggio che può derivare dal loro uso e il grado di trasferibilità su un campione di popolazione più ampio. 6. I metodi presentati non tengono in alcun modo conto del tasso di ricombinazione e degli aspetti demografici (ad esempio collo di bottiglia) che possono modificare la struttura dei pattern di LD. 34 Capitolo 3 Stratificazione di popolazione Le recenti proposte tra le tecnologie del sequenziamento e l’aumentata disponibilità di marker genici hanno aperto la strada agli studi di associazione genome wide su larga scala. Un potenziale problema che nasce da ogni studio basato su popolazione è la presenza non identificata di una stratificazione di popolazione che può simulare un segnale di associazione spurio e di conseguenza portare a falsi positivi e/o alla non individuazione dei reali effetti [5,8,9,10]. Quando casi e controlli hanno frequenze alleliche diverse attribuibili a diversi background di popolazione, e non ad eventuali espressioni fenotipiche, allora si dice che lo studio è caratterizzato da stratificazione di popolazione. La stratificazione di popolazione è probabilmente la ragione più spesso citata della mancata riproducibilità dei risultati dei test di associazione. 3.1 Cause della stratificazione La stratificazione può essere dovuta principalmente a tre cause che sono: la presenza di una struttura (cioè di sottogruppi) all’interno di una popolazione (population structure, Figura 3.1), presenza di legami famigliari tra i soggetti campionati (family structure), e presenza di correlazioni, anche lontane, tra i soggetti con legami famigliari di grado superiore al primo (cryptic relatedness). Tipicamente, gli studi di associazione genome-wide evitano di campionare, all’interno dello stesso studio, individui provenienti da diverse popolazioni con l’intenzione di evitare la stratificazione di popolazione. In generale, per ogni studio caso-controllo, il pool di popolazione da cui i casi vengono campionati dovrebbe essere lo stesso dal quale anche i controlli vengono campionati. 35 Figura 3.1: Effetti della struttura di popolazione su un locus di uno SNP. Se lo studio di popolazione coinvolge sottopopolazioni che differiscono geneticamente, e se anche la prevalenza della malattia differisce tra queste sottopopolazioni, allora le proporzioni di casi e controlli campionati da ogni sottopopolazione tenderà ad essere diversa, dal momento che le frequenze alleliche o genotipiche tra casi e controlli ad ogni locus saranno divere per ogni sottopopolazione. La Figura mostra un esempio di questo scenario con due popolazioni dove i casi sono costituiti da un eccesso di individui provenienti dalla popolazione 2 e la popolazione 2 ha una frequenza minore dell’allele A rispetto a quella rilevata nella popolazione 1. In questo esempio, la struttura simula il segnale di associazione data dalla significativa differenza nelle frequenze alleliche e genotipiche tra casi e controlli causata non dal fenotipo, ma dalle caratteristiche della popolazione stessa. 3.2 Identificazione della stratificazione I metodi, di seguito brevemente esposti, propongono una possibile soluzione al problema di stratificazione. Questi richiedono un numero sufficiente di SNPs (preferibilmente maggiore di 100) che siano stati sequenziati nei casi e nei controlli in aggiunta agli SNPs candidati nel test di associazione. 3.2.1 Genomic control (GC) Il metodo più comunemente utilizzato, soprattutto nell’identificazione della presenza di stratificazione, è il Genomic Control (GC). È basato sull’osservazione che la population 36 structure cambia la distribuzione nulla della statistica χ2 di un fattore moltiplicativo λ, che viene calcolato basandosi su una collezione di marker come il rapporto: (g6 ) f= 0.456 a numeratore viene calcolata la mediana della statistica S2= λ χ2 , ed è divisa per la media teorica calcolata in ipotesi nulla. In assenza di stratificazione, l’associazione tra le varianti genetiche e la malattia dovrebbe seguire una distribuzione χ2 a 1 grado di libertà (1 df), pertanto il valore di λ risulta circa uguale a 1. In presenza di stratificazione, il valore di λ è maggiore di 1, condizione che si verifica in quanto che la stratificazione tende a incrementare la statistica di un fattore costante. Generalmente, valori di λ minori di 1.05 vengono ancora tollerati come non critici, e si può considerare pressochè assente la stratificazione. Il valore di λ varia al variare della dimensione del campione, sia in funzione del numero di loci considerati nel calcolo del GC che del numero di individui campionati [8,9,10,11]. La dipendenza di λ dalle dimensioni del campione (intese come numero di individui) è mostrata in Figura 3.2. Risulta evidente che, a parità di condizioni di stratificazione, all’aumentare delle dimensioni del campione, aumenta anche il valore di λ. Dalla Figura si nota anche che il valore di λ risulta unitario in assenza di stratificazione e tende invece ad aumentare man mano la presenza di stratificazione diventa più marcata, in accordo con quanto detto prima. Infine, quando per identificare la stratificazione vengono usati solo un piccolo numero di loci (da 50 a 100) allora il test GC è spesso non conservativo, il p-value calcolato è più grande di quello dato dalla distribuzione χ2 e ancora risultano falsi positivi. E’ stato dimostrato che il GC non corregge adeguatamente la stratificazione di popolazione se per stimare il fattore di correzione sono utilizzati troppo pochi loci (quindi pochi marker). Al contrario, quando il numero di loci è molto grande, il test è molto conservativo e tende ad eliminare i falsi positivi dal risultato, ma ne consegue una perdita di potere statistico in alcune applicazioni. La scelta del numero di SNPs necessari per la corretta soluzione al problema di stratificazione dovrebbe tenere conto anche dell’intensità degli effetti genetici causati dal marker in studio. 37 Figura 3.2: effetti della stratificazione sugli studi di associazione. La stratificazione influisce sulla statistica χ2 di un fattore λ, che cambia in dipendenza delle dimensioni del campione. Lo scenario 1 corrisponde a una condizione di evidente stratificazione, lo scenario 4 corrisponde a assenza di stratificazione, gli scenari 2 e 3 corrispondono a condizioni intermedie, dove la stratificazione è rispettivamente più e meno marcata. In ordinata si ha il valore di λ, mentre in ascissa si ha la dimensione del campione per uno studio in cui si suppone che il numero dei casi sia uguale al numero dei controlli. λ1000 fa riferimento alla situazione in cui il numero di casi e controlli è uguale a 1000. 3.2.2 Principal component analysis (PCA) La Principal Component Analisys (PCA) rappresenta un approccio multivariato che permette di rilevare differenze in termini di strutture genetiche tra sottogruppi di individui. Se applicata a dataset con dati campionati da individui provenienti da diverse sottopopolazioni, allora l’analisi è in grado di dare una interpretazione geografica ai dati, individuando graficamente in modo distinto i gruppi relativi alle sottopopolazioni. I soggetti sono caratterizzati in uno spazio vettoriale come punti distinti, in modo tale che gli individui geneticamente simili formino “nuvole” compatte (cluster), mentre gli individui geneticamente più dissimili si posizionino in regioni più distanti dello spazio (Figura x.3). La presenza di cluster può aiutare a individuare sottopopolazioni di individui omogenei tra loro che possono essere considerati separatamente nelle analisi di associazione aumentandone il potere statistico. 38 Figura 3.3: rappresentazione delle prime due componenti principali dell’analisi PCA (rispettivamente in ascissa e ordinata) che mettono a confronto le caratteristiche di background genetico dei diversi soggetti con diversi progenitori. Ogni punto rappresenta un individuo analizzato, mentre il colore fa riferimento al diverso gruppo etnico. Il metodo della PCA utilizza le componenti principali come covariate per correggere la stratificazione negli studi GWA: un modello semplice lineare rappresenta il fenotipo Y come funzione degli effetti fissi X: h = ij + M Qui, X=[Gi PC1i …PCmi] indica il genotipo del marker candidato (Gi) in aggiunta a covariate opzionali, come ad esempio il sesso o l’età, B denota invece i coefficienti degli effetti fissi e ε è il termine di rumore che tiene conto della variazione di Y non spiegata. La PCA applica una correzione più grande ai marker che presentano forti differenze alleliche. A differenza delle implementazioni iniziali, la PCA è diventata un’analisi computazionalmente trattabile anche su grandi dataset. Approcci correlati a questo, come il multidimensional scaling (MDS) e il matching genetico sono stati implementati da diversi software come PLINK [riferimenti]. Una limitazione del metodo è l’incapacità di tener conto nel modello della presenza di family structure e cryptic relatedness. 39 3.2.3 Mixed models I modelli misti (mixed models) possono modellare tutte e tre le cause della stratificazione di popolazione, cioè population structure, family structure e cryptic relatedness. L’approccio alla base è la modellizzazione dei fenotipi (Y) considerando contemporaneamente effetti fissi (X) e effetti random (u). Il modello lineare è dato dal h = ij + + M Gli effetti fissi comprendono lo SNP candidato e covariate opzionali (B), come sesso o età, mentre gli effetti random sono basati su una matrice delle covarianze fenotipiche, che a sua volta è modellata come la somma di variazione random ereditabile e non ereditabile. Nel modello, u denota una componente della varianza del rumore u+ε che è distribuita secondo la matrice K di Kinship. X () = k ! l Qui, u rappresenta la componente ereditabile della variazione random mentre ε rappresenta la componente non ereditabile. La matrice di Kinship K è definita secondo la similarità genotipica tra coppie di individui, quindi la sua struttura è sicuramente influenzata dalla struttura di popolazione, family structure e cryptic relatedness. Il parametro k ! esprime quanto gli individui, geneticamente simili, sono anche fenotipicamente simili. Più sono simili i due individui, più è alto il valore della correlazione. 40 Capitolo 4 Il progetto internazionale HapMap 4.1 Cos’è HapMap HapMap (Haplotype - Map) è un database, disponibile gratuitamente in internet, catalogo delle variazioni genetiche comuni che si osservano negli esseri umani e ne descrive i pattern comuni. Il progetto ha permesso l’emergere e lo sviluppo degli studi di associazione genome wide, e si ritiene possa essere una risorsa chiave per i ricercatori nel tentativo di identificare geni che possono avere ripercussioni sulla salute, malattie, risposte ai farmaci e fattori ambientali. Il progetto nasce da una collaborazione cominciata nell’Ottobre del 2002 tra diversi paesi quali: Giappone, Inghilterra, Canada, Cina Nigeria, e Stati Uniti. L’obiettivo che i ricercatori si sono proposti con questo progetto è determinare i pattern comuni nelle variazioni della sequenza di DNA nel genoma umano attraverso la caratterizzazione di queste stesse varianti (individuate negli SNPs), la stima della frequenza con cui compaiono e il grado di correlazione tra queste, utilizzando campioni di DNA provenienti da popolazioni di discendenza africana, asiatica ed europea. Malattie comuni come le patologie cardiovascolari, cancro, obesità, diabete, malattie psichiatriche o infiammatorie sono causate dall’azione combinata di fattori genetici e ambientali. Le cause di queste malattie non sono imputabili all’azione di un unico gene, bensì viene riconosciuta una predisposizione genetica, a carico degli aplotipi, che si combina poi con fattori ambientali. Le variazioni del genoma possono quindi servire da marker genici per determinare l’associazione tra una particolare regione del genoma e la malattia. Un approccio di questo tipo può essere velocemente applicato grazie ad HapMap a qualsiasi gene candidato nel genoma, o a qualsiasi regione che può essere individuata da studi di linkage familiare, o, in ultimo, all’intero genoma per l’individuazione di fattori di rischio. HapMap propone una importante scorciatoia 41 nell’individuazione dei geni candidati e negli studi di linkage e di associazione basati su l’intero genoma. Nel suo scopo, l’International HapMap Project ha molto in comune con il Human Genome Project, il quale sequenzia tutto il genoma umano. Ma mentre quest’ultimo progetto si occupa di sequenziale l’intero genoma, incluso il 99.9% di genoma che tutti abbiamo in comune, il progetto HapMap si occupa di caratterizzare i pattern comuni all’interno del 0.1% che ci differenzia l’uno dall’altro. 4.2 Realizzazione del progetto HapMap Il progetto è di fatto divenuto possibile grazie al confluire delle seguenti conoscenze: la disponibilità della sequenza del genoma umano (si è fatto riferimento a diversi database già esistenti, come ENCODE), un database di SNPs comuni (dbSNP, successivamente arricchito dal progetto stesso), l’intuizione sul linkage disequilibrium (LD) e lo sviluppo di tecnologie accurate e relativamente a basso costo per il sequenziamento high- throughtput degli SNPs. Il lavoro ha avuto effettivamente inizio nel 2002 con la realizzazione di una prima fase che ha interessato 269 campioni di DNA prelevati da individui provenienti da Nigeria, Utah (USA), Cina e Giappone, dai quali sono stati sequenziati 1.007.329 SNPs. La seconda fase ha consentito di sequenziale complessivamente 4.6 milioni di SNPs in ognuno dei campioni HapMap. Attualmente è in corso una terza fase che coinvolge un numero di campioni ancora maggiore, proveniente da diversi gruppi etnici (inclusi paesi come il Messico, altri stati degli USA, il Kenya e anche l’Italia). L’elenco completo delle popolazioni coinvolte nel progetto è riportato in Tabella 4.1. I campioni hanno solo un’etichetta che ne identifica la provenienza etnica e il sesso del donatore; non è possibile in nessun caso risalire, attraverso le informazioni disponibili, all’individuo. Le varie fasi si differenziano non solo per la quantità di campioni e donatori coinvolti nel progetto, ma anche nella densità con cui vengono sequenziati gli SNP, la stima della minor allele frequency (MAF), e i pattern di linkage disequilibrium (LD). Infatti nella fase due, gli SNPs sono sequenziati più densamente e sono stati inclusi anche quelli che presentano in media una MAF minore rispetto a quella della fase 1: questo ha contribuito a un significativo miglioramento nella rappresentazione di variazioni rare rispetto a quanto fatto inizialmente. Le fasi del progetto e le relative caratteristiche vengono riassunte schematicamente in Tabella 4.2. 42 Etichetta Campione di popolazione # campioni QC campioni America (USA) del sud ovest con antenati ASW (A)* 90 71 Africani Residenti in Utah (USA) con antenati dal nord e CEU (C)* 180 162 ovest europa CHB (H) pop. Han di Pechino, Cina; 90 82 CHD (D) cinesi in Metropolitan Denver, Colorado 100 70 GIH (G) indiani Gujarati residenti a Houston, Texas 100 83 JPT (J) giapponesi di Tokyo, Japan 91 82 LWK (L) Luhya in Webuye, Kenya 100 83 MEX(M)* Pop. con antenati messicani a Los Angeles (Ca) 90 71 MKK(K)* Maasai in Kinyawa, Kenya 180 171 TSI (T) Toscani, Italia 100 77 YRI (Y)* Yoruban in Ibadan, Nigeria 180 163 1,301 1,115 Tabella 4.1: elenco delle popolazioni coinvolte nel progetto(aggiornato alla fase 3). La colonna Etichetta si riferisce all’abbreviazione utilizzata come riferimento alla popolazione. La presenza dell’asterisco indica che i campioni sono stati relativi a nuclei famigliari composti da 3 individui (madre, padre e figlio). La colonna QC riporta il numero di campioni che vengono mantenuti dopo il filtraggio Quality Control- Le cifre nell’ultima riga si riferiscono rispettivamente al totale dei campioni e dei campioni filtrati. Fasi Descrizione Fasi Raccolta campioni da 4 popolazioni (CEU, YRI, CHB, JPT) per un Fase 1 totale di 269 individui; #SNPs: circa 1 milione, approssimativamente circa 1 ogni 5 kb, MAF>0.05 Raccolta dati intensificata nelle stesse popolazioni, #SNPs: vengono Fase 2 aggiunti più di 3.1 milioni di SNPs per un totale di 4.6 milioni, 1 ogni 1 kb approssimativamente, MAF>0.05 Raccolta dati nelle popolazioni che in questa fase arrivano a 11. Si Fase 3 contano più di 1.6 milioni di SNPs nuovi. Tabella 4.2: in questo schema vengono riportate le fasi principali del progetto con una breve descrizione nella colonna a destra 43 Con l’avanzare del progetto non solo vengono aggiunti al database nuovi SNPs, ma vengono anche ridefiniti SNP la cui posizione e/o funzione prima risultava inesatta. I gruppi di ricerca che si occupano del sequenziamento del DNA seguono un uguale e rigido protocollo per il controllo della qualità dei dati e per l’identificazione degli SNPs. Si stima che l’accuratezza per fase di sequenziamento sia in media del 99.5% per i vari gruppi; tuttavia c’è un alto tasso di dati mancanti e discrepanze nel sequenziamento stesso. 4.3 Pubblicazione e consultazione dei dati Il progetto è impegnato in un rilascio rapido e completo dei dati, e si assicura che i risultati siano di dominio pubblico a nessun costo per l’utenza. Tutti i dati relativi ai nuovi SNPs, alle frequenze alleliche e genotipiche sono consultabili e scaricabili dal sito dell’ HapMap Data Coordination Centre (http://www.hapmap.org). La ricerca dei contributi genetici alle malattie umane generalmente si concentra su geni candidati identificati sulla base di studi di linkage e/o associazione, o sulla base di pathway che si presume siano coinvolti in un particolare aspetto della malattia in esame. Nello studio dei geni candidati, un ricercatore vorrà conoscere se ci sono SNPs comuni nelle immediate vicinanze, quali alleli hanno, e quali sono le relative frequenze alleliche nella popolazione. La sezione Data del sito consente un accesso interattivo al database attraverso un browser (Generic Genome Browser) che permette agli utenti di ricercare un particolare gene, o una regione di interesse di piccole o medie dimensioni, all’interno del genoma e quindi di visualizzare la distribuzione di SNPs e modelli di variazione comune nella regione stessa. La ricerca viene effettuata utilizzando il nome di una sequenza, di un gene, locus o altri punti di riferimento. I principali risultati su cui HapMap permette agevolmente di indagare sono: Identificazione SNPs e stima delle loro proprietà; Distribuzione della frequenza allelica dei campioni provenienti da un gruppo etnico; Frequenza allelica degli SNPs tra le varie popolazioni; Aplotipi condivisi dai diversi gruppi etnici; Valutazione del LD nel genoma umano e variazioni nel tasso di ricombinazione; Selezione di tagSNP per studi di associazione. 44 Viene di seguito presentato un esempio di ricerca che mostra un possibile utilizzo del Genomic Browser. Dopo l’accesso al sito www.hapmap.org e alla sezione Data, si entra nel browser vero e proprio. La query si effettua andando a digitare, nel campo Cerca, il termine selezionato. Può essere usato indifferentemente uno di questi termini: il nome del cromosoma (ad esempio “Chr10”); la posizione in intervallo di basi sul cromosoma nel formato Cromosoma:start…stop (ad esempio “Chr10: 114,700,201…114,916,051”); il nome dello SNP utilizzando l’identificativo “rs” (ad esempio “rs081062”); il nome del gene secondo la nomenclatura del NCBI RefSeq (ad esempio “NM 153254”); il nome comune del gene (ad esempio “BRCA2”); la banda cromosomica (ad esempio “5q31”). Una volta eseguita la ricerca, verrà mostrata la pagina con i risultati. Se si hanno più corrispondenze alla query effettuata, allora la pagina mostrerà tutti i risultati specificandone la locazione sul genoma. Per default, il browser propone i risultati aggiornati alla data più recente; in alternativa è possibile selezionare l’origine dei dati, aggiornati in corrispondenza delle varie fasi del progetto. In particolare il menù consente di consultare i dati relativi alla fase 2, alla fase 3 o entrambe, con differenti date di rilascio dei dati stessi. Questa scelta non è banale, in quanto condiziona le opzioni di scaricamento e ricerca. Si vuole indagare, a scopo illustrativo, sul fattore di trascrizione Trascription Factor 7- like2 ,TCF7L23. 4.3.1 Esempio: Ricerca di TCF7L2 Nel caso in esame è sufficiente scrivere il nome TCF7L2; in alternativa si può indicare il numero del cromosoma e l’intervallo di basi che ne individuano la posizione. In questo esempio vengono consultati i dati relativi alla fase 3 aggiornati al febbraio 2009. Una volta inviata la richiesta, il DB individua automaticamente la posizione del fattore di trascrizione, riportando il numero del cromosoma in cui si trova e l’intervallo di basi 3 TCF7L2 è un fattore di trascrizione comune noto anche come TCF4. Influenza la trascrizione di diversi geni, quindi svolge molteplici di funzioni all’interno della cellula. È anche implicato in diverse malattie; diversi SNPs sono associati in particolar modo al diabete di tipo 2. Nella popolazione europea è stato identificato come il principale fattore di rischio per questa patologia. 45 in cui il fattore è codificato. Il progetto permette di applicare metodi di analisi nuovi o già esistenti per l’analisi e la visualizzazione dei dati. Di seguito sono presentati i risultati grafici e sperimentali di tale ricerca. Panoramica Il DB indica che il fattore si trova sul cromosoma 10. La regione occupata dal fattore di trascrizione in esame è evidenziata in giallo in Figura 4.1 e interessa la posizione 114,700,201 - 114,916,051 individuabile sul righello che suddivide l’intero cromosoma in Mb. Il grafico gt’d SNPs/500Kb riporta il numero di SNPs sequenziati da HapMap ogni 500Kb sul cromosoma. Figura 4.1: panoramica della regione TCF7L2 Regione Il riquadro sottostante (Figura 4.2) ripresenta in scala diversa la porzione evidenziata in Figura 4.1, concentrandosi sulla zona che interessa il fattore di trascrizione scelto. La sezione del grafico gt’d SNPs/20Kb riporta il numero di SNPs individuati ogni 20 Kb. Figura 4.2: regione del cromosoma 10 occupata dal fattore di trascrizione 46 Dettagli Nella sezione Dettagli vengono proposti diversi tipi di analisi; per default inizialmente viene mostrata solo una parte di tutta l’informazione disponibile. La sezione più importante è la Genotype SNPs (Figura 4.3). I vari SNPs sequenziati dal progetto sono indicati con un triangolo equilatero in corrispondenza della loro posizione sul cromosoma (Figura 4.3). Il grado di dettaglio con cui vengono visualizzati cambiano a seconda dello zoom imposto. Figura 4.3: visualizzazione grafica della posizione degli SNPs individuabile per mezzo del righello posizionato sopra Andando a cambiare scala è possibile individuare con maggior dettaglio gli SNPs e leggerne il nome e le frequenze alleliche stimate per ogni popolazione. Se si sceglie ad esempio di indagare sull’intervallo di posizione 114,730,000 - 114,800,000 troviamo il seguente risultato (Figura 4.4). Posizionando il cursore su ogni singolo SNP o cliccando su esso, si apre una tabella che riporta le frequenze genotipiche e alleliche. I dati sono suddivisi in base alla popolazione di provenienza per le quali vengono utilizzate le seguenti sigle: Selezionando, ad esempio, lo SNP rs081062 si ottengono le informazioni nella tabella riportata in Figura 4.5. 47 Figura 4.4: visualizzazione più dettagliata degli SNPs compresi all’interno della regione 114,730,000-114,800,000 Figura x.5: tabella che si visualizza selezionando lo SNP rs7981062 Si può visualizzare la percentuale di contenuto di G/C nel DNA nella regione selezione (Figura 4.6). 48 Figura 4.6: contenuto di G e C nella sequenza del fattore di trascrizione Nella pianificazione di uno studio di associazione è essenziale la conoscenza dell’estensione del linkage disequilibrium (LD) nella regione target per la riduzione del numero di SNPs da sequenziare. La determinazione dei pattern di LD è stato uno dei principali obiettivi del progetto HapMap. I dati possono essere scaricati in blocco dal sito o consultati interattivamente utilizzando il browser. Per quanto riguarda la seconda opzione, nella sezione Dettagli è possibile vedere graficamente la distribuzione dei pattern di LD nella regione selezionata. I parametri chiave sono il tipo di misura che si intende usare per il calcolo del LD (parametri D’, r2 e il LOD score) [ref capitolo], l’orientazione del triangolo con cui si visualizza il pattern (se con il vertice in alto o in basso), lo schema di colori., e se le dimensioni dei box nel plot devono essere proporzionali alla distanza genomica tra i marker o di dimensione uniforme. I pattern di LD possono essere visualizzati per una o più popolazioni. In Figura x.7 si visualizza il Plot del LD ricavato con i soli dati dei campioni della popolazione CEU e utilizzando le due misure più comuni: D’ e r2. In alternativa è possibile visualizzare il plot per ognuna delle altre popolazioni. Il passo successivo alla determinazione dei pattern di LD, è la scelta dei tagSNP. Per piccole regioni è possibile selezionare i tagSNPs “a mano”, utilizzando i supporti grafici e numerici generati sopra. Tuttavia il miglior risultato è garantito sempre dall’utilizzo di un algoritmo che scegla i tagSNPs massimizzando formalmente il numero di SNPs in LD catturati nel tag-set. Non esiste un unico insieme di tagSNP che soddisfa le esigenze del probelma. Saranno i ricercatori a selezionare quali SNPs lavorano meglio con un particolare sistema di genotipizzazione, e a effettuare il compromesso tra costi di sequenziamento di uno studio di popolazione e la forza del livello di associazione che possono identificare. 49 Figura 4.7: visualizzazione dei plot di LD utilizzando i soli dati relativi al gruppo CEU. Il riquadro superiore è ottenuto plottando i valori di LD calcolati con r 2 , mentre il riquadro inferiore è ottenuto plottando i valori di LD misurati tramite D’. Il plot a triangolo è costruito congiungendo ogni coppia di SNPs su una linea orientata di 45° rispetto l’orizzontale. Il colore più intenso indica un grado maggiore di LD, mentre le zone colorate in grigio indicano dati mancanti. Per questo motivo, il sito HapMap non offre un set di tagSNPs preselezionati, ma invece offre ai ricercatori un mezzo per selezionare interattivamente i tag basandosi su criteri scelti dall’utente. La lista di tagSNPs è generata da algoritmi supportati dal programma Tagger. (http://www.broad.mit.edu/mpg/tagger/, de Bakker et al. 2005). Nella sezione Scaricamento, Ricerca e altre operazioni si sceglie Annota TagSNP Picker per identificare i TagSNP. Se si seleziona Configura, è possibile settare le diverse opzioni che includono: la scelta di una popolazione e di un algoritmo (Tagger Pairwise o Tagger Multimaker), l’inclusione o esclusione di una lista di SNPs dai tag da selezionare, inclusione di una lista di “punteggi” che pesano diversamente i maker, selezione di soglie di cutoff sui valori di LD accettabili e le frequenze alleliche per gli SNPs da includere. Il metodo Tagger Pairwise, sviluppato da Carlson et al, seleziona uno SNP come tag se presenta alta correlazione con un altro. Il metodo Tagger Multimaker invece utilizza un approccio basato non sul confronto a coppie come per il caso precedente, ma su predittori multi-marker e per questo risulta avere una maggiore efficienza. Oltre al tipo di algoritmo, è possibile settare anche il valore minimo del coefficiente r2 con il quale selezionare lo SNP (se impostato a 1 si ottiene un set di SNP 50 non ridondante) e il valore della minor frequency allele (MAF). Il risultato viene mostrato graficamente nel riquadro inferiore di Figura x.8, dove nella sezione tagSNP Picker, vengono riportati i tagSNPs del fattore di trascrizione TFC7L2. Nel riquadro superiore di Figura 4.8, è presente invece il riquadro relativo alla OMIM disease associations. Posizionando il cursore sullo SNP evidenziato e cliccando su esso, si accede alla pagina web OMIM (Online Mendelian Inheritance in Man) relativa al fattore di trascrizione TFC7L2 dove si può leggere non solo quanto riguarda il fattore stesso, ma come quel particolare SNP è implicato nella predisposizione genica verso determinate malattie. Nel caso in esame l’articolo riporta brevemente i risultati più importanti conseguiti dai gruppi di ricerca con i relativi riferimenti bibliografici. Figura 4.8: visualizzazione di OMIM disease association e i TagSNP Per la visualizzazione degli aplotipi, il sito HapMap si appoggia al programma PHASE versione 2.1 (Stephens e Donnely, 2003). In Figura 4.9 sono riportati gli aplotipi relativi al fattore di trascrizione ricostruiti statisticamente utilizzando il software PHASE. Nel plot vengono mostrati tutti i 120 cromosomi con gli alleli colorati in giallo o in blu. 51 Figura 4.9: visualizzazione degli aplotipi di due popolazioni. Nel riquadro superiore è rappresentato l’aplotipo del gruppo CEU, mentre quello inferiore è relativo al gruppo CHB. Durante il phasing, ogni allele in un genotipo è assegnato a uno o all’altro cromosoma parentale usando un algoritmo di maximum likelihood che usa l’informazione sui nuclei familiari di tre persone (genitori e un figlio) nei gruppi della popolazione HapMap. Se questa informazione non è disponibile (non è presente cioè informazione sui gruppi familiari, ma solo su individui singoli), si fittano i dati in un modello che minimizza il numero di crossover nella popolazione. Gli aplotipi phased sono visualizzabili graficamente mediante un codice a due colori. Ogni cromosoma è rappresentato come una linea di altezza 1 pixel e ogni allele dello SNP è arbitrariamente colorato in blu o giallo. Una regione ad alto LD apparirà come una regione nella quale ci sono lunghe file di SNPs in cui gli alleli hanno il medesimo colore, indicando una bassa ricombinazione. Una regione a basso LD apparirà come un’area dove invece i segmenti sono più corti e più frammentati. L’ordine dei cromosomi è determinato da un algoritmo di clustering gerarchico che raggruppa i cromosomi che condividono che condividono aplotipi simili. In Figura x.9 vengono messi a confronto gli aplotipi dei gruppi CEU e CHB. 52 4.4 HapMart L’analisi condotta finora ha consentito di visualizzare i dati prevalentemente da un punto di vista grafico. Per avere listati gli SNPs di una determinata regione di interesse e/o di un determinato gruppo etnico, con le relative caratteristiche (frequenza allelica, genotipica..), si accede ad HapMart. HapMart è una versione modificata di BioMart, un sistema di data management orientato alle query. BioMart è stato sviluppato in modo che i ricercatori potessero eseguire query anche complesse, consultando i maggiori database di sequenze biomolecolari, pathway, e di annotazione, come Ensembl, Uniprot, Reactome HGNC, Wormbase e PRIDE. HapMart utilizza la stessa interfaccia grafica di BioMart e la stessa modalità di esecuzione della query, ma quest’ultima è limitata al database di HapMap, pertanto la ricerca risulta circoscritta ai soli SNPs del progetto. Nella sezione di HapMart è possibile selezionare in dettaglio i criteri con cui effettuare la query. In particolare : Tipo di database (anche se di fatto la scelta è limitata a un'unica risorsa); popolazione (posso considerare tutte le popolazioni o un gruppo soltanto); il valore della MAF (%); supporto per il sequenziamento dei campioni: Perlegen amplicon-based platform (che sequenzia gli SPNs da frammenti di DNA amplificati con PCR), Affimetrix GeneChip Mapping Array, Illumina HumanHap100, MIP,..; la regione genica di interesse (numero cromosoma e/o intervallo di basi); inclusione o meno di determinati SNP. Si selezionano infine i dettagli, nella sezione Attributes, che si vogliono inclusi nel report finale: ID, cromosoma di appartenenza, posizione, alleli, allele di riferimento,..; codice della popolazione, genotipo; frequenza allelica e genotipica dello SNP. Nel risultato finale vengono riportati per default i primi 10 risultati; è possibile scaricarne una versione anche in formato Excel. La stessa ricerca può essere effettuata con BioMart, selezionando il database relativo ad HapMap. 53 4.5 Ricerca per malattia Se nel campo ricerca del browser si inserisce il nome della malattia, ad esempio diabetes, il database provvede a dare come risultato tutte le regioni su ogni cromosoma che risultano coinvolte con questa malattia. In particolare, per ogni singolo risultato sono indicati il numero del cromosoma e l’intervallo di basi che interessano la regione. Se disponibile viene riportata una breve descrizione che spiega cosa codifica quella sequenza, e come è coinvolta nella patologia. Nel caso in esame del diabete, la query produce 401 risultati. Andando ad indagare per ognuno di essi, si trovano collegamenti ad altri database (Reactome, BioXRT,..) che permettono di studiare come e in quali processi sono coinvolte le sequenze. Il database permette quindi di selezionare ogni singolo risultato e ottenere un’analisi uguale a quella condotta in precedenza per il fattore di trascrizione TCF7L2. 54 Capitolo 5 PLINK Gli studi di associazione sull’intero genoma (WGAS) hanno comportato una nuova sfida computazionale e analitica per i ricercatori. Molti supporti già esistenti per l’analisi genetica non sono stati progettati per maneggiare un così ampio data set in modo pratico ed efficiente e non riescono a spiegare la complessità di indagine che deriva dall’utilizzo dei dati dell’intero genoma. Uno degli strumenti più diffusi per analizzare questo tipo di dati è PLINK (di Shaun Purcell, http://pngu.mgh.harvard.edu/~purcell/plink/), uno strumento open-source, che permette di condurre le analisi di routine in modo computazionalmente più efficiente, e offre la possibilità di introdurre nuovi metodi che sfruttano al meglio le potenzialità di data set così grandi. Le principali funzionalità di PLINK WGAS permettono di : provvedere a un modo semplice per gestire grandi set di WGAS; stimare gli errori dovuti al problema della stratificazione della popolazione [ref] e ai genotipi errati [ref]; operare una varietà di test di associazione standard in modo efficiente su grandi data set (su una popolazione, su famiglie, con o senza covariate, test su aplotipi,…). 5.1 Formato dei dati PLINK I dati utilizzati per l’analisi con PLINK sono tipicamente di due tipi: mydata.ped e mydata.map. 55 5.1.1 PED files Un file.ped (Tabella 5.1) è caratterizzato dalla presenza di 6 colonne obbligatorie, dove ogni individuo è rappresentato da una riga, caratterizzate dai seguenti campi: 1. Family ID 2. Individual ID 3. Paternal ID 4. Maternal ID 5. Sex (1 = maschio, 2 = femmina, other = non specificato) 6. Phenotype. Gli ID sono alfanumerici; la combinazione degli ID dei campi relativi all’appartenenza a un gruppo famigliare e all’individuo (rispettivamente prima e seconda colonna) deve identificare univocamente una persona. La sesta colonna, Phenotype, deve essere caratterizzata dalla presenza di un solo codice relativo a un preciso fenotipo che può assumere uno dei seguenti valori: 1 corrisponde alla presenza del fenotipo, 2 indica l’assenza del fenotipo, 0, -9 o missing genotype identificano, in forma alternativa, la mancanza di informazione sulla presenza o assenza del fenotipo. Quest’ultimo può essere una caratteristica, uno stato o una malattia: PLINK individua automaticamente lo status dell’individuo basandosi sul valore osservato in corrispondenza di questo campo. Nel caso in cui il sesso dell’individuo non sia noto, allora viene identificato con un carattere diverso da 1 e 2. Se nel data set non sono presenti informazioni sul fenotipo o sul sesso dell’individuo, questo sarà automaticamente escluso dalle analisi che fanno uso di queste stesse informazioni. Ad esempio, verrà dato un messaggio di errore nel caso in cui si sta eseguendo un’indagine nella quale si richiede di ricostruire il nucleo famigliare e uno degli individui, di cui manca il sesso, deve essere identificato come padre o madre. In caso in cui sia disponibile anche l’informazione sul genotipo, questa viene riportata dalla settima colonna in poi, e può essere denotata con un codice di caratteri qualsiasi (ad esempio 1,2,3,4 o A,C,T,G) ad eccezione dello 0 che indica il genotipo mancante. I c

Studio di associazione Genome-Wide: Preprocessing e Selezione SNPs (PDF)

Document Details

Tags

Related

Summary

Full Transcript

Upgrade to continue