Esercitazioni di bioinformatica e genomica funzionale

ES 1 – Identificazione di Trascription Factor Binding Sites (TFBS)

Calcolare il coefficiente informativo di ogni posizione di un sito di legame partendo dall'allineamento di 60 sequenze nella matrice A del file 'alignment'. Dato che dice che il fattore di trascrizione considerato si può legare ad entrambi gli strend dobbiamo considerare anche le sequenze complementari (che vanno aggiunte alla matrice A).

Procedimento

load ('alignment.mat');

% Aggiungo alla matrice A le sequenze complementari
A2 = A;
for j = 1:60
    A2 = [A2; seqrcomplement(A2(j,:))];
end

% Usare A2 per calcolare la weight matrix
wm = zeros(4,21); % Dimensioni: 4 basi in 21 allineamenti
for i = 1:21
    wm(1,i) = length(strfind(A2(:,i), 'a'));
    wm(2,i) = length(strfind(A2(:,i), 'c'));
    wm(3,i) = length(strfind(A2(:,i), 'g'));
    wm(4,i) = length(strfind(A2(:,i), 't'));
end
wm = wm/120;

% Uso la weight matrix per calcolare il

contenuto informativo (dopo aver calcolato l'entropia);

NB. dato che devo calcolare il log in base 2 delle frequenze non posso avere valori uguali a 0!

wm(wm==0)=1;

H=-wm.*log2(wm);

H=sum(H);

Commentato [MC29]: Sommo sulle 4 basi per ottenere l'entropia totale

I=2-H;

• Eseguo e ottengo un vettore I di 21 elementi, più il contenuto è alto più sono certa di trovare quella base in quella posizione.

Trovo quali posizioni hanno I>0.8

find(I>0.8)

• Ottengo 4 posizioni che sono la 4, la 8, la 14 e la 18. Visualizzo I ad esempio con stem.

stem (I)

• Visualizzo graficamente la significatività delle lettere e il contenuto informativo di ogni singola posizione

seqlogo(wm)

ES 2 – sequence logo

Utilizzare una funzione specifica per i siti di legame ovvero sequence logo con seqlogo.

Procedimento

seqlogo(wm)

ottengo per ogni posizione dell'allineamento tutte le basi che si possono trovare in una posizione e l'altezza complessiva di una

ES 3 - Find TFBSPrendiamo una sequenza nucleotidica che corrisponde al promotore in genere e andiamo a cercare in questa sequenza il sito di legame del nostro fattore di trascrizione. Teniamo wm, scorriamo la frequenza e calcoliamo segnale per ogni posizione; confrontiamo la frequenza in esame con la wm e più sono simili più è alto il segnale complessivo: S(n)=SUM(I(l)*f(b,l)La frequenza ci viene data nel file fna che è un file fasta e non mat quindi lo leggo con una

• Ci chiede di calcolare il segnale da 400 basi a monte a 100 basi a valle e la base da cui inizia la trascrizione è la 3395 quindi vado ad isolare la parte che ci interessa

• Controllare 400 basi a monte e 100 a valle rispetto al fattore di trascrizione (#3395)

• Eseguo e cerco le posizioni in cui il segnale è maggiore di 4.5

• Trovo 4 posizioni che sono la 64, 108,189, 233 e posso plottare il

segnalex=1:500-21+1;plot(x,S,’b’);hold onplot(x,4.5*ones(size(x)),’r’)

LABORATORIO 5Esercizi su espressione genica quindi dati da microarray.

ES 1 – Preprocessing di microarray a cDNA per visualizzarne i due canali rosso e verde

Procedimentoload (‘microdata.mat’)

• Rimuovere dai canali del rosso e del verde tutto quello che è dovuto al rumore di fondo: valorinegativi e di bassa qualità (-50)
• R = Rsig-Rbkg;G = Gsig-Gbkg;neg_valueR = find(R<0);neg_valueG = find(G<0);neg_value = union(neg_valueR,neg_valueG);neg_quality = find(flags==-50);datogliere = union(neg_value,neg_quality);R(datogliere) = []; Commentato [MC31]: Rimuovo ponendo gli indici deglielementi da togliere uguali ad un vettore vuotoG(datogliere) = [];
• Fare scatter plot di log2(R) vs log2(G) per vedere se serve una normalizzazioneplot(log2(G),log2(R),’.’)xlabel(‘log2G’)ylabel(‘log2R’)hold onplot([2 16],[2 16],

` per creare una lista puntata per il procedimento.

‘r’)Vedo che i punti sono molto sbilanciati verso il basso della diagonale, quindi, hanno bisogno di una normalizzazione.

• Visualizzo anche con MA plot

figure
plot(1/2*(log2(R)+log2(G)),log2(R./G),'.')

Vedo che i punti sono sbilanciati verso il basso della linea dello zero, quindi, hanno bisogno di una normalizzazione. Un altro modo per ottenere l’MA plot posso usare una funzione matlab:
[x,y] = mairplot(R,G)

• Normalizzo con la normalizzazione lowess

ysmooth = malowess(x,y);
hold on
plot(x,ysmooth,'rx')

La linea rossa (rx) è quella che dopo la normalizzazione dovrà diventare la linea orizzontale sullo zero; vado quindi a raddrizzare questa linea:

ynorm = y-ysmooth;
figure
plot(x,ynorm,'g.')

Dal nuovo plot si vede che i dati sono distribuiti in maniera più uniforme intorno all’asse orizzontale. In alternativa si poteva usare la funzione mairplot che apre una finestra da cui è possibile visualizzare

```html

lasmooth curve e normalizzare direttamente col bottone normalize.

ES 2 – array oligonucleotidi con un solo canale

Studio di array affymetrix con 42 chip (campioni) e 140k probes (sonde). Sono date intensità di sonde con match perfetto e mismatch.

Procedimento

1. Rimuovere il background e visualizzare l'istogramma della distribuzione dell'intensità
pms_bg = rmabackadj(probeData.PMIntensities, ‘showplot’, 1);

Restituisce istogramma del primo campione, da cui si vede in rosso la stima della densità di background

2. Utilizzare il quantilenorm per normalizzare le intensità di PM
pms_bgnorm = quantilenorm(pms_bg);
3. Plottare le intensità grezze, dopo la rimozione del background e dopo la normalizzazione; uso subplot così li visualizza insieme
samples = probeData.CELNames; figure subplot(3,1,1)
```

deicampioni con valore di intensità più alto abbiamo più trascritto.

ES 3 – Studio espressione genica da dati già processati: confronto differenziale

Confronto differenziale tra soggetti sani e patologici per comprendere quali geni vengono espressi maggiormente in patologia e per capire cosa succede nelle cellule patologiche. (sano/malato, due stadi diversi della patologia, due diversi tessuti)

Useremo il ptest che dice se la media di due gruppi di valori è statisticamente significativa.

Procedimento

• Scaricare i dati di espressione: vado su ncbi e cerco nella barra GSE5847, scarico il file txt
• Importare i dati su Matlab: gseData = geoseriesread('GSE5847_series_matrix.txt')
• Dividiamo i dati tra tessuto stromale ed epiteliale per trovare i geni più espressi nello stromale rispetto all'altro: sampleSources = unique(gseData.Header.Samples.source_name_ch1); stromaIdx = strcmpi(sampleSources{1},

gseData.Header.Samples.source_name_ch1);

nStroma = sum(stromaIdx);

• T test per Analisi differenziale

[pvalues] = mattest(gseData.Data(:,stromaIdx),gseData.Data(:,~stromaIdx));

• Fare il volcano plot e correggere il pvalue con correzione di Bonferroni per ridurre in numero di falsipositivi

mavolcanoplot(gseData.Data(:,stromaIdx),gseData.Data(:,~stromaIdx),pvalues)

pvalues_Bonf = pvalues*22283; 22283 è il numero di test

mavolcanoplot(gseData.Data(:,stromaIdx),gseData.Data(:,~stromaIdx),pvalues_Bonf)

Nel secondo grafico si vede che tutti i valori sono stati schiacciati verso il basso; siamo più sicuri che i genitrovati non sono falsi positivi dovuti solo al gran numero di campioni. Analogo di bonferroni: invece dimoltiplicare pvalue correggo la soglia di cutoff, alzandola.

mavolcanoplot(gseData.Data(:,stromaIdx),gseData.Data(:,~stromaIdx),pvalues, 'PCutoff', 0.05/22283)

LABORATORIO 6

Esercizi su espressione genica e analisi dati su clustering e PCA (per ridurre

complessità dati) + geneontologyES 1 – clustering gerarchico

Procedimento

Importare i dati per un esperimento a tempo (timedata.mat). Identificare 16 cluster distinti utilizzando un raggruppamento gerarchico:

load('timedata.mat')

Ottenere 16 cluster con una strategia di clustering gerarchico

corrdist = pdist(yestvalues, 'corr')

clustertree = linkage(corrdist,'avarage')

I dati sono già normalizzati e ogni gene rappresenta una dimensione dello spazio su cui fare clustering

cluster_labels = cluster (clustertree, 'maxclust',16)

Plottare i clusters clustering: metrica di linkage --> come considero la distanza tra due cluster? se h