Analisi di dati Microarray: Esercitazione...

Introduzione Analisi di una singola ibridazione Analisi di più ibridazioni PCA e Visualizzazione scientifica

Analisi di dati Microarray:Esercitazione Matlab

Laboratorio di Bioinformatica II

Pietro Lovato

Anno Accademico 2011/2012


Contenuti

1 IntroduzioneIntroduzione a MatlabBioinformatics ToolboxDNA Microarray

2 Analisi di una singola ibridazioneAmbiente MatlabSottrazione del backgroundScatter plot

3 Analisi di più ibridazioniEsempio: Serie temporali di lievitoGene SelectionCluster analysisValidare i risultati

4 PCA e Visualizzazione scientificaPrincipal Component Analysis


Introduzione a Matlab

Che cos’èMATrix LABoratory

Sistema interattivo per il calcolo tecnico e scientifico

Linguaggio di programmazione

Internamente: librerie numeriche altamente ottimizzate

Una caratteristica di Matlab è l’espansibilità delle sue funzioni,attraverso delle soluzioni denominate toolbox. Tali toolbox sonocollezioni complete di funzioni Matlab per risolvere particolaricategorie di problemi.


I toolbox

Statistica

Elaborazione di immagini

Elaborazione di segnali

Pattern recognition

Bioinformatica

Sistemi di comunicazione

Sistemi di controllo

Economia & Finanza

Calcolo parallelo

. . .


Applicazioni per la bioinformatica

1 Analisi di sequenzeAllineamento, primer design, identificazione di motivi

2 Analisi di dati microarrayAnalisi di immaginiNormalizzazione dei datiClustering & data mining

3 ProteomicaAnalisi di dati da spettrometria di massa

4 Systems BiologySimulazione di pathway metabolici

5 Filogenesi


DNA Microarray

Che cos’è

Tecnologia capace di misurare simultaneamentel’espressione di migliaia di geni.Sono analizzate differenti condizioni sperimentali:

Diversi tessuti (es. da individui sani/malati)Diverse condizioni di crescitaDiversi istanti di tempo

Enorme quantità di dati, da analizzare con strumentiinformatici.


DNA Microarray (2)


DNA Microarray (3)

ProblematicheDesign delle sondeProblematiche di Image Processing:

Segmentazione spotRimozione rumoreRilevamento della qualità degli spotQuantificazione del segnale

Analisi dei dati:Riduzione della dimensionalitàCluster analysisVisualizzazione scientifica

Interscambio di formati


DNA Microarray (4)

ProblemaA seconda del tipo di esperimento edella piattaforma utilizzata, i dati checi vengono forniti possono averediversi formati.

Qualche esempio

.cel (Affymetrix)

.chp (Affymetrix)

.gpr (GenePix)

.soft (Dataset GEO)

.txt...


L’interfaccia grafica

Sotto Applicazioni −→ Programmazione


Un tipico comando Matlab

data_mouse = gprread(

’mouse_a1pd.gpr’)

;



Risultato



;



Risultato



;

Nome del comando



Risultato



;

Nome del comando

Argomenti in input



Risultato



;

Nome del comando

Argomenti in input

NotaPer i comandi di visualizzazione non è necessario specificareun risultato.


Lavorare con una singola ibridazione

Esempio

Dati di espressione di topo in cui è stato indottofarmacologicamente il morbo di Parkinson.

Ogni chip contiene l’espressione in una regione (voxel) delcervello.

Dati grezzi: file in formato .gpr oppure soft (GEOaccession GSE30).

File .gpr già presente come demo Matlab.

Caricare i datipd = gprread(’mouse_a1pd.gpr’);


Strutture in Matlab

struct()

Tipo di dato che raggruppa più campi in un unico oggetto.

Un campo può memorizzare informazioni di vario tipo.

Accesso ai campi con variabile.nomeCampo.

NotaFormati di file microarraydiversi avranno campicon nomi diversi (mal’informazione contenutaè sostanzialmente lastessa).

Header: [1x1 struct]Data: [9504x38 double]

Blocks: [9504x1 double]Columns: [9504x1 double]

Rows: [9504x1 double]Names: {9504x1 cell}...


Cosa ci interessa

Campi

Il campo Data è una matrice che contiene i dati.

Il campo Names contiene i nomi dei geni (le righe).

Il campo ColumnNames contiene i nomi delle colonne.

All’interno del campo ColumnNames:

’F635 Median’ intensità del canale rosso (cy5,λ = 635nm).

’F635 Median - B635’, come sopra ma consottrazione del background.

’F532 Median - B532’ intensità del canale verde (cy3,λ = 532nm).

La mediana è uno stimatore robusto dell’intensità dei pixel.


Sottrazione del background

Background: fluorescenza che puòcontribuire all’intensità misurata manon è dovuta all’ibridazione del target.

Può dipendere da polvere, molecoleisolate, imperfezioni nel sistema diacquisizione, rumore, ...

Il background varia lungo il supporto eva stimato localmente.

B532 Median

Provare a visualizzare il bg del canale verde con il comandomaimage(pd, ’B532 Median’);


Salvare in variabili separate i campi di interesse

Sintassi

cy5Data = magetfield(pd,’F635 Median - B635’);cy3Data = magetfield(pd,’F532 Median - B532’);genes = pd.Names;


Scatter plot: R vs G

Dati “osservati”: {(R,G)i}


Scatter plot: log2 R vs log2 G

Dati “osservati”: {(log2 R, log2 G)i}


Scatter plot: I vs R

Dati “osservati”:{(

log10(R · G), log2

(RG

))

i

}


In Matlab

Per creare uno scatter plot del tipo log R vs log G:

maloglog(cy5Data,cy3Data,’labels’,genes);

Per creare uno scatter plot del tipo I vs R:

mairplot(cy5Data,cy3Data,’labels’,genes);


In Matlab

Per creare uno scatter plot del tipo log R vs log G:

maloglog(cy5Data,cy3Data,’labels’,genes);

Per creare uno scatter plot del tipo I vs R:

mairplot(cy5Data,cy3Data,’labels’,genes);

EsercizioProvare a disegnare uno scatter plot “a mano”:

Dato un vettore x , il logaritmo delle sue componenti sicalcola con il comando: variabile = log10(x);Un qualsiasi scatter plot (non microarray-specifico) sidisegna con il comando: scatter(x,y);


Ricapitolando...

» pd = gprread(’mouse_a1pd.gpr’);

» maimage(pd, ’B532 Median’);

» cy5Data = magetfield(pd,’F635 Median - B635’);

» cy3Data = magetfield(pd,’F532 Median - B532’);

» genes = pd.Names;

» maloglog(cy5Data,cy3Data,’labels’,genes)

» mairplot(cy5Data,cy3Data,’labels’,genes)


Idea: combinare dati provenienti da diverse ibridazioni


Cluster analysis

Passi (a grandi linee)

Identificare geni differenzialmente espressi:t-test se il problema è a due classi.Altri indici per selezionare geni “rilevanti” (varianza,entropia, ...).

Definire una distanza tra profili di espressione.

Scegliere un algoritmo di clustering.

Validare il risultato.


Esempio: Serie temporali di lievito

ProblemaInvestigare come variano nel tempo i livelli di espressionedel lievito (Saccharomyces cerevisiae).

Si è interessati alla riprogrammazione metabolica nelpassaggio dalla fermentazione alla respirazione.

IdeaCercare di raggruppare geni con profili di espressione“simili”.

Cercare di caratterizzare ciascun gruppo attraverso ilprocesso biologico cui i geni prendono parte.

Caricare i datiload yeastdata;


Plottare un profilo di espressione

plot(

times, yeastvalues(

15 , :)

)

;



Matrice dei dati

plot(

times, yeastvalues(

15 , :)

)

;



Matrice dei dati

plot(

times, yeastvalues(

15 , :)

)

;

Seleziono la quindicesima riga=gene



Matrice dei dati

plot(

times, yeastvalues(

15 , :)

)

;

Seleziono la quindicesima riga=gene

Seleziono tutte le colonne=tempi


Gene Selection

Scartare i geni “non informativi”

Geni a bassa varianza.

Geni con basso valore assoluto in ogni istante temporale.

Geni a bassa entropia informativa.

Sintassi» [mask, yeastvalues, genes] =genevarfilter(yeastvalues, genes);

» [mask, yeastvalues, genes] =genelowvalfilter(yeastvalues, genes, ’absval’, log2(3));

Quanti geni rimangono dopo la selezione?


Introduzione: Cluster Analysis

Una possibile definizione

Il clustering rappresenta l’organizzazione di un insieme di entitàin gruppi (clusters) sulla base della similarità.

Il concetto di cluster è vago:

Dipendentemente dalle misure di similarità utilizzatecambia il risultato.La scelta della misura di similarità è cruciale

Dovrebbe essere fatta in modo da inglobare la maggiorquantità possibile di informazione a priori.

Il risultato può cambiare anche a seconda dellametodologia utilizzata per fare clustering.


Clustering partizionale vs gerarchico

Clustering partizionale

Il risultato è una singola partizione dei dati.

Tipicamente il numero di cluster deve essere dato a priori.

Esempio: K-means.

Clustering gerarchico

Il risultato è una serie di partizioni innestate (undendrogramma).

Non è necessario settare a priori il numero di cluster.

Più informativo del clustering partizionale ma improponibileper dataset molto grandi.

Esempi: Complete link, single link, ...


Clustering gerarchico

Funzione clustergram

clustergram(data, ’Linkage’, ’tipo’, ’rowlabels’, genes);

Dove:data: matrice di espressione.’tipo’: Algoritmo di clustering gerarchico.

’single’: Algoritmo Single Link’complete’: Algoritmo Complete Link’average’: UPGMA...

genes: identificativi dei geni.


Clustering partizionale

K-meansAlgoritmo più famoso di clustering partizionale.

Ogni gene è un punto in uno spazio n-dimensionale (nnumero di esperimenti).

Ogni cluster è rappresentato dalla sua media.Idea:

L’algoritmo è iterativo: si “assesta” in più passi.Si parte da una clusterizzazione iniziale, e ad ogniiterazione si assegna ogni gene alla media più vicina.Si riaggiornano le medie.Si ripete fino a convergenza.


K-means clustering


K-means clustering (2)

Funzione kmeans

[cidx, ctrs] = kmeans(data, K, ’dist’, ’tipo’);

Input:

data: matrice di espressione.

K : numero di cluster. Per questo esempio K=16.’tipo’: Distanza da utilizzare.

’sqEucl’: Distanza euclidea’corr’: Distanza di correlazione statistica...

Output:

cidx: indica, per ogni gene, il cluster a cui appartiene.

ctrs: centri dei cluster – profilo “medio” di espressione.


Visualizzare i risultati

Funzione maplot

Per visualizzare i profili appartenenti ad ogni cluster:

maplot(times, yeastvalues, cidx)

Per visualizzare il profilo medio di espressione in ognicluster:

maplot(times, yeastvalues, ctrs)

Quale/i cluster rappresenta geni che si esprimono solo durantela respirazione cellulare?


Validare i risultati: GOstat

Una volta individuato un cluster di interesse x , esportare inun file di testo i geni che lo compongono:cluster = genes(cidx==x);dlmwrite(’cluster.txt’, char(cluster), ’’);

Carico il file su GOstat:Nel campo “Gene-association database” mettere sgd(lievito)(Opzionale) Nel campo “Subset of GO hierarchy” metterebiological_process (siamo interessati al pathway cui igeni nel cluster prendono parte)

In output ho i termini GO che meglio descrivono i geniall’interno del mio cluster.

Che informazioni riesco a dedurre?


Ricapitolando...

» load yeastdata

» plot(times, yeastvalues(15,:))

» [mask, yeastvalues, genes] =genevarfilter(yeastvalues, genes);

» [mask, yeastvalues, genes] =genelowvalfilter(yeastvalues,genes,’absval’,log2(3));

» clustergram(yeastvalues,’Linkage’,’complete’, ...’rowlabels’,genes)

» [cidx, ctrs] = kmeans(yeastvalues, 16, ’dist’,’corr’);

» maplot(times, yeastvalues, cidx)

» maplot(times, yeastvalues, ctrs)

» cluster = genes(cidx==1);

» dlmwrite(’cluster.txt’, char(cluster), ’’);


Principal Component Analysis (PCA)

MotivazioniRidurre la dimensione dei dati (qui ogni dato è un genecaratterizzato da 7 numeri)...

per visualizzare il dataset;per ridurre il carico computazionale delle tecnicheapplicate;per eliminare possibile ridondanza;

Mettere in evidenza particolari strutture presenti nei dati.


PCA (2)

IdeaPresentare i dati nel riferimento che evidenziamaggiormente la loro struttura.

Estrarre le direzioni di massima varianza dei dati.


PCA (3)

Idea: eliminare le direzioni a bassa varianzaIn questo modo si comprimono i dati minimizzando laperdita di informazione.

Per trasformare ogni punto: formule e nozioni derivantidall’algebra lineare (autovalori e autovettori della matrice dicovarianza...).


PCA (4)

Vantaggi

Migliore tecnica di compressione dei datiMigliore in senso di errore quadratico medio

Svantaggi

Non è detto in tutti i casi che le direzioni a varianzamaggiore siano le direzioni ottimali;

Si perde l’interpretabilità delle grandezze rappresentate.


PCA in Matlab

Funzione princomp

[pc, comp] = princomp(data);

Input:

data: matrice di espressione.

Output:

pc: coefficienti della trasformazione.

comp: coordinate dei punti nello spazio trasformato,ordinate per importanza decrescente.


PCA in Matlab (2)

Plottare le prime due o tre componenti per visualizzare i datitrasformati.

Il comando è sempre scatter:scatter(comp(:,1),comp(:,2));

o, in alternativa:scatter3(comp(:,1),comp(:,2),comp(:,3));


Esercizio

Effettuare un clustering k-means sui dati trasformati,prendendo solo le prime due o tre componenti principali;

Ad es, per prendere solo le prime due componenti

data_pca = comp(1:2,:);

K = 2;

Vengono separati i due cluster?


Ricapitolando...

» [pc, comp] = princomp(yeastvalues);

» scatter(comp(:,1),comp(:,2));

» scatter3(comp(:,1),comp(:,2),comp(:,3));

» pca_data = comp(:,1:2);

» [pca_idx, pca_ctrs] = kmeans(pca_data, 2,

’dist’, ’sqEucl’);

» scatter(pca_data(:,1), pca_data(:,2), 4,

pca_idx);

NotaMatlab offre un help molto semplice e dettagliato!

Analisi di dati Microarray: Esercitazione...

Documents

Transcript of Analisi di dati Microarray: Esercitazione...