Biomedical Engineering - Bioinformatica e Genomica Funzionale

Exercises

Complete course

1 2 %%%%%%%% LAB0 3 4 %%%%%%%% LAB1 5 6 %%%%%%%% LAB2 - filogenetica 7 8 %% ESERCIZIO 1 9 10 % Compute the probability of the tree depicted in figure as the sum of the probability of individual trees corresponding to possible reconstructions (mutually exclusive events). Assume the following background frequencies and substitution model (contenuti nel file 'ex1.mat'): 11 12 % fb = [0.3 0.2 0.2 0.3] 13 14 % A C G T 15 % A | 0.4 0.2 0.2 0.2 | 16 % M = C | 0.2 0.4 0.2 0.2 | matrice con le probabilità di 17 % G | 0.2 0.2 0.4 0.2 | sostituzione delle diverse basi. 18 % T | 0.2 0.2 0.2 0.4 | 19 20 % dobbiamo far variare le basi nei nodi intermedi (chiamati i,j,k) e calcolare la probabilità di tutte le possibili combinazioni di basi. Poi sommiamo la probabilità delle singole combinazioni e otteniamo la probabilità complessiva per ottenere questo singolo albero. 21 22 load ('ex1.mat' ) 23 24 P=0; 25 26 for i=1:4 % stiamo considerando tutte le possibili combinazioni di basi che portano al nostro albero finale CCAT 27 for j=1:4 28 for k=1:4 29 p=fb (i)* M(i,j)* M(j,2)* M(j,2)* M(i,k)* M(k,1)* M(k,4); %probabilità di una singola combinazione 30 % i nodo centrale 31 % j nodo sinistra che collega C e C 32 % k nodo destra che collega A e T 33 % M(j,2) indica ramo che collega j alla colonna 2 di M, ovvero base C nota 34 P=P+p; % sommiamo poichè sono eventi mutualmente esclusivi. 35 end 36 end 37 end 38 39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 40 41 %% ESERCIZIO 2 - buiding a phylogenetic tree for the hominidae species 42 43 % The mitochondrial D-loop is one of the fastest mutating sequence regions in animal DNA, and therefore, is often used to compare closely related organism. 44 % 1) Given the pre-aligned sequences contained in primates.mat, compute pairwise distances using the 'Jukes-Cantor' formula and the phlogenetic tree with the 'UPGMA' distance method. 45 % 2) You can explore the phylogenetic tree by considering the nodes (leaves and branches) within a patristic distance of 0.3 from the 'European Human' entry and reduce the tree to the sub-branches of interest by pruning away non-relevant nodes. 46 47 48 load ('primates.mat' ) 49 50 distances = seqpdist (primates ,'Method' ,'Jukes-Cantor' ,'Alpha' ,'DNA' ) 51 % parametro = method = Jukes-Cantor 52 % sequenza in ingresso = alpha = DNA (nucleotidi di DNA) 53 % distanze come vettore: [(2,1),(3,1),...,(M,1),(3,2),...,(M,M-1)] Distanze tra la seconda sequenza e la prima, terza e la prima, quarta e la prima,..., ultima con la prima. Poi seconda coppia: terza con la seconda, quarta con la seconda, quinta con la seconda,...fino ad arrivare all'ultima con la penultima. 54 55 % 1) costruisco l'albero 56 UPGMAtree = seqlinkage (distances , 'UPGMA' , primates ); %eseguo e vedo che nel workspace appare UPGMAtree 57 % posso visualizzare l'albero con la funzione plot 58 h = plot (UPGMAtree , 'orient' ,'top' ); 59 title ('UPGMA Distance Tree of Primates using Jukes-Cantor model' ); 60 ylabel ('Evolutionary distance' ); % titolo dell'asse ylabel 61 62 % 2) FILTRARLO --> voglio mantenere solo le sequenze con una distanza entro 0.3 rispetto all'organismo di riferimento (european human). 63 % Salvo i nomi delle foglie dell'albero nella variabile "names" 64 names = get (UPGMAtree , 'LeafNames' ); 65 66 % Seleziono i nodi che ci interessano, quelli che hanno una distanza utilizzo la funzione prune 74 leaves_to_prune = ~ h_leaves ; 75 pruned_tree = prune (UPGMAtree ,leaves_to_prune ); 76 77 % Plotto questo nuovo albero ancora con orientamento verticale 78 h = plot (pruned_tree , 'orient' , 'top' ); 79 title ('albero filtrato' ) 80 view (UPGMAtree , h_leaves ); % vedo in rosso queste foglie sull'albero totale 81 82 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 83 84 %% ESERCIZIO 3 - bootstrapping phylogenetic trees 85 86 % The topology of every bootstrap tree is compared with that of the original tree. Any interior branch that gives the same partition of species is counted. Since branches may be ordered differently among different trees but still represent the same partition of species, it is necessary to get the canonical form for each subtree before comparison. In a canonical tree the leaves are ordered alphabetically and the branches are ordered first by their width and then alphabetically by their first element. 87 % 1) Apply bootstrap analysis with 100 replicates. 88 % The confidence information associated with each branch node can be stored within the tree by annotating the node names. Thus, you can create a new tree, equivalent to the original primates tree, and annotate the branch names to include the confidence levels computed above. 89 % 2) Select the branch nodes with a confidence level greater 0.9. 90 91 % Quanto siamo sicuri di questa relazione? Assegnamo un valore di confidenza --> BOOTSTRAP 92 93 num_boots = 100 ; % numero di iterazioni che vogliamo utilizzare 94 num_seqs = length (primates ); % numero di sequenze nel nostro allineamento 95 seq_len = length (primates (1). Sequence ); % lunghezza dell'allineamento. Ho preso la prima sequenza tanto sono tutte uguali. 96 97 % 1) Definiamo un array di celle che contenga gli allineamenti fittizi 98 boots = cell (num_boots ,1); 99 100 for n = 1:num_boots % itero sulle 100 iterazioni e sulle 12 sequenze 101 reorder_index = randsample (seq_len , seq_len , true ); % campiono l'allineamento originale. Metto true perchè voglio che ci sia rimpiazzo. Parametri: 102 % 1° --> insieme di campionamento 103 % 2° --> lunghezza del vettore casuale che otteniamo in uscita 104 for i = num_seqs :- 1:1 % vogliamo riordinare le sequenze secondo il nuovo ordine 105 bootseq (i). Header = primates (i). Header ; % bootseq allineamento fittizio 106 bootseq (i). Sequence = strrep (primates (i). Sequence (reorder_index ), '-' , '' ); % riordiniamo le posizioni secondo il nuovo indice e rimuoviamo i gap --> vogliamo riallineare le sequenze prima di calcolare di nuovo l'albero fittizio. 107 end 108 %salvo questa struct nell'array 109 boots {n} = bootseq ; 110 end 111 112 % Definisco un nuovo array di celle dove salvare gli alberi --> array di celle con 100 righe e una colonna 113 boot_trees = cell (num_boots ,1); 114 115 for n = 1:num_boots %vettore temporaneo distanze, calcolo tutto come prima 116 dist_tmp = seqpdist (boots {n}, 'Method' , 'Jukes-Cantor' , 'Alpha' , 'DNA' ); 117 boot_trees {n} = seqlinkage (dist_tmp , 'UPGMA' , boots {n}); 118 end 119 120 for i = num_seqs -1:- 1:1 121 branch_pointer = i + num_seqs ; %variabile che contiene il numero del nodo che voglio usare come radice 122 sub_tree = subtree (UPGMAtree ,branch_pointer ); %ricavo il sottoalbero partendo da branch_pointer come radice. Salvo i sottoalberi in un array di celle 123 orig_pointers {i} = getcanonical (sub_tree ); %getcanonical mi dà la forma canonica 124 orig_species {i} = sort (get (sub_tree , 'LeafNames' )); %salvo le specie contenute in questo sottoalbero in ordine alfabetico grazie alla funzione sort 125 end 126 127 for j = num_boots :- 1:1 %itero su tutti i 100 alberi 128 for i = num_seqs -1:- 1:1 129 branch_ptr = i + num_seqs ; 130 sub_tree = subtree (boot_trees {j}, branch_ptr ); 131 cluster_pointers {i,j} = getcanonical (sub_tree ); 132 cluster_species {i,j} = sort (get (sub_tree , 'LeafNames' )); 133 end 134 end 135 136 % Conto quante volte ogni ramo dell'albero originale ricorre nell'albero di bootstrap 137 count = zeros (num_seqs -1,1); 138 for i = 1:num_seqs -1 139 for j = 1:num_boots *( num_seqs -1) % numero totale sottoalberi 140 if isequal (orig_pointers {i}, cluster_pointers {j}) % il primo if controlla che la topologia sia uguale 141 if isequal (orig_species {i}, cluster_species {j}) 142 count (i) = count (i) + 1; 143 end 144 end 145 end 146 end 147 148 % RICAVO LA CONFIDENZA --> CALCOLO LE FREQUENZE 149 Pc = count ./ num_boots ; 150 151 % 2) SELEZIONARE RAMI CON CONFIDENZA > 0.9 152 [ptrs , dist , names ] = get (UPGMAtree ,'POINTERS' ,'DISTANCES' ,'NODENAMES' ); 153 % AGGIORNO IL VALORE DEI NODI AGGIUNGENDO LA CONFIDENZA 154 for i = 1:num_seqs -1 155 branch_ptr = i + num_seqs ; 156 names {branch_ptr } = [ names {branch_ptr } ', confidence: ' num2str (100 *Pc (i)) ' %' ]; 157 end 158 159 %RIDEFINISCO L'ALBERO USANDO QUESTI NUOVI NOMI 160 tr = phytree (ptrs ,dist ,names ); %ottengo un albero equivalente all'originale ma con names aggiornato con la confidenza del ramo rappresentato dal quel nodo. 161 162 %SELEZIONO I RAMI CON CONFIDENCE > 0.9 --> Per ottenere l'indice corretto devo aggiungere numseq 163 high_conf_branch_ptr = find (Pc > 0.9 ) + num_seqs ; 164 % Visualizzo il nuovo albero coi nodi in evidenza 165 view (tr ,high_conf_branch_ptr ) 166 167 ////////////////////////////////////////////////////////////////////// 168 169 %%%%%%%% LAB3 - information theory 170 171 %% ESERCIZIO 1 - identification of potential transcription factor binding site (TFBS) 172 173 % alignment.mat contains the alignment of 60 FIS DNA binding sites. 174 % 1) Compute the information content I for the 21 positions of the site. Since Fis binds to DNA on both strands, consider both the sequences and their reverse complement. Which positions have I>0.8 bits? 175 % NB: I = 2 - [H(l) + e(n)] ipotizziamo di avere un numero di sequenze tale da poter trascurare l'errore e(n) 176 % NB: entropia di Shannon H(l) = -sum(i=1->M)[Pi*log2(Pi)] 177 % 2) Reproduce the sequence logo by using seqlogo 178 % 3) Scan the E. coli sequence contained in K02672.fna to identify potential binding sites. Consider 400 bases upstream and 100 downstream in respect of the transcription start site (base #3395 in the file). If you use a threshold of 4.5 bits, how many positives do you find in this region? In which positions? 179 % NB: segnale S = sum(l=1->L)[I(l)*f(b,l)] 180 % NB: f(b,l) --> frequenza della base che troviamo nella sequenza in posizione l. 181 182 % [useful matlab functions: seqrcomplement, strfind, fastaread] 183 184 load ('alignment.mat' ) %matrice di 60 sequenze ciascuna con 21 posizioni 185 186 % 1) Per calcolare la weight matrix e il contenuto informativo vogliamo considerare sia queste 60 sequenze che le loro sequenze complementari. 187 A2 = A; % definiamo una nuova matrice che contenga sequenze + sequenze complementari. NB: Cerchiamo di non sovrascrivere una matrice, perciò usiamo A2. 188 189 for j = 1:60 %scorro le sequenze 190 A2 = [ A2 ; seqrcomplement (A2 (j,:))]; %aggiunge la sequenze complementari, trovate con la funzione seqrcomplement (riga j tutte le colonne). 191 end 192 193 % Calcolo la weight matrix a partire da questa matrice 194 wm = zeros (4,21 ); % 21=posizioni, 4=basi --> inizializzo la matrice a 0 195 196 for i = 1:21 % per ognuna delle colonne leggo quante volte occore una base, e così riempio la wm con le occorrenze. Iteriamo sulle posizioni, non sulle sequenze (21 posizioni) 197 wm (1,i) = length (strfind (A2 (:, i)', 'a' )); 198 % A2(:,i) selezioni la colonna i della matrice A2, tutte le righe 199 % A2(:,i)' faccio la trasposta poichè mi serve un vettore riga 200 % strfind cerca il secondo parametro all'interno del primo, che diamo come vettore riga poichè dobbiamo scriverlo come una stringa = vettore orizzontale di caratteri. Restituisce tutte le posizioni in cui la trova, ma a noi interesa solo quante volte viene trovata, quindi calcoliamo la lunghezza di questo vettore output (length()). 201 wm (2,i) = length (strfind (A2 (:, i)', 'c' )); 202 wm (3,i) = length (strfind (A2 (:, i)', 'g' )); 203 wm (4,i) = length (strfind (A2 (:, i)', 't' )); 204 end 205 206 wm = wm /120 ; % determino le frequenze di occorrenza dividendo per il numero totale di sequenze. 207 wm (wm == 0) = 1; % sostituisco 0 con 1 nella matrice per non avere problemi col log --> tanto comunque log1=0 quindi il risultato non cambia. 208 H = - wm .* log2 (wm ); % entropia di Shannon 209 H = sum (H); % sommo i valori delle 4 basi 210 211 % Calcolo il contenuto informativo I. Non considero l'errore 212 I = 2 - H 213 sol = find (I>0.8 ); % Mi da le posizioni in cui I è maggiore di 0.8 214 stem (I); % plottiamo il contenuto informativo (stem = steli di fiori) 215 title ('stem graphic' ); 216 217 % 2) Riprodurre il sequence logo della weight matrix. 218 seqlogo (wm ); % ogni colonna per posizione (21 posizioni). Per ogni colonna l'altezza complessiva rappresenta il contenuto informativo di quella posizione. La dimensione relativa delle lettere rappresenta la frequenza delle basi all'interno di quella posizione. 219 220 % 3) Determiniamo il segnale 221 promoter = fastaread ('K02672.fna' ); % promoter è adesso la sequenza tutta intera in formato fasta 222 % Voglio 500 basi totali, 400 all'inizio e 100 alla fine rispetto alla base 3395 nel file 223 l = length (promoter .Sequence ); 224 % sovrascrivo per ottenere solo l'intervallo di interesse nel campo Sequence (Header non ci interessa) 225 promoter = promoter .Sequence (3395 -400 :3395 +99 ); 226 227 S = zeros (1,500 -21 +1); % creo il vettore segnale e lo inizializzo. Ha queste dimensioni perchè non possiamo calcolare il segnale fino all'ultima base, perchè sto procedendo a finestre di 21 caratteri. 228 229 for j = 1:500 -21 +1 % sposto la finestra 230 for p = 1:21 % scorro i caratteri della finestra 231 switch promoter (j+p-1) 232 case 'A' %è come scrivere if promoter(j+p-1)=='A' 233 S(j) = S(j) + I(p) * wm (1,p); 234 case 'C' 235 S(j) = S(j) + I(p) * wm (2,p); 236 case 'G' 237 S(j) = S(j) + I(p) * wm (3,p); 238 case 'T' 239 S(j) = S(j) + I(p) * wm (4,p); 240 end 241 end 242 end 243 244 % Richiesta esercizio: posizioni con un segnale > 4.5 245 find (S>4.5 ); % in queste posizioni la sequenza è maggiormente simile al sito di legame. 246 % Posso plottare il segnale. 247 X = 1:500 -21 +1; 248 plot (X,S); 249 hold on 250 plot (X,4.5 *ones (size (X)), 'r' ); % linea orizzontale per vedere la soglia 251 252 % Ora in effetti vedo 4 picchi che superano la linea di soglia 253 254 ////////////////////////////////////////////////////////////////////// 255 256 %%%%%%%% LAB4 257 258 % TEORIA: 259 260 % La C in una coppia CpG viene spesso modificata dalla metilazione (ovvero, un Hatom viene sostituito da un gruppo CH3). C'è una probabilità relativamente alta che il metilC muti in una T. Quindi, le coppie CpG sono sottorappresentate nel genoma umano. I dinucleotidi CpG del DNA si verificano circa cinque volte meno frequentemente del previsto. 261 262 % Le isole CpG sono regioni non metilate del genoma associate alle estremità 5' della maggior parte dei geni di pulizia e di molti geni regolati. L'assenza di metilazione rallenta il decadimento di CpG, quindi le isole CpG possono essere rilevate nella sequenza del DNA come regioni in cui le coppie di CpG si verificano vicino alla frequenza prevista (vediamo più coppie di CpG in esse che altrove). 263 % - Circa l'80% delle isole CpG sono comuni all'uomo e al topo. 264 % - Circa il 56% dei geni umani e il 47% dei geni murini sono associati a CpG. 265 266 % Definizione classica: sequenza di DNA di lunghezza 200 con un contenuto di C + G del 50% e un rapporto tra il numero di CpG osservato e atteso superiore a 0,6. 267 268 %% ESERCIZIO 1 - an HMM for detecting CpG islands 269 270 % The model contains eight states that can emit the four symbols of the DNA alphabet --> Σ={A,C,G,T} 271 % Given a (short) sequence, how to decide whether this segment is from a CpG-island or not? 272 % Load the sequence CpGtest contained in ex1.mat and calculate the probability that the entire sequence is predicted to be in a CpG island and the probability that the entire sequence is not predicted to be in a CpG island. Which is the highest? 273 274 %aprire ex1 e state transition matrix 275 %split transition matrix 276 load ('ex1' ) 277 load ('state_transition_matrix' ) 278 279 % vogliamo separare il mondo delle GpC island da quello normale 280 trans_matrix_IN_CpG = a(1:4,1:4); % quadrante alto sinistra (mondo GpC) 281 trans_matrix_OUT_CpG = a(5:8,5:8); % quadrante basso destra (mondo normale) 282 % preparo la variabile in cui metto la probabilità cumulativa, ovvero la produttoria del modello markoviano. La inizializzo a 1 perchè è un prodotto 283 P_sequence_IN_CpG = 1; 284 P_sequence_OUT_CpG = 1; 285 286 L = length (CpGtest ); 287 sequence = upper (CpGtest ); % rendiamo la sequenza tutta maiuscolo 288 % setto lo stato di partenza per tutte e due 289 switch (sequence (1)) % k tiene conto dello stato di partenza della produttoria 290 case 'A' 291 k = 1; 292 case 'C' 293 k = 2; 294 case 'G' 295 k = 3; 296 case 'T' 297 k = 4; 298 end 299 % con lo switch guardo il primo simbolo della sequenza poi itero su tutta la sequenza 300 % read sequence position by position 301 for pos = 2:L % l tiene conto dello stato di arrivo della produttoria 302 symbol = sequence (pos ); 303 switch (symbol ) 304 case 'A' 305 l = 1; 306 case 'C' 307 l = 2; 308 case 'G' 309 l = 3; 310 case 'T' 311 l = 4; 312 end 313 314 % verifico a(k,l) --> k stato precedente, l stato nuovo 315 a_IN_CpG = trans_matrix_IN_CpG (k,l); % 316 a_OUT_CpG = trans_matrix_OUT_CpG (k,l); 317 318 % update produttorie 319 P_sequence_IN_CpG = P_sequence_IN_CpG * a_IN_CpG ; 320 P_sequence_OUT_CpG = P_sequence_OUT_CpG * a_OUT_CpG ; 321 322 % passo allo stato successivo (prossima iterazione del ciclo) 323 k = l; 324 end 325 326 % Quale delle due probabilit è più elevata? 327 if P_sequence_IN_CpG > P_sequence_OUT_CpG 328 disp ('la sequenza è in una CpG island' ) 329 else 330 disp ('no' ) 331 end 332 333 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 334 335 %% ESERCIZIO 2 336 337 % Given a (long) sequence, how to find CpG-islands contained in it? 338 % Load the sequence contained in the ex2.mat and try to identify a potential CpG island. Compare the results with UCSC Genome browser annotations. 339 340 % ALGORITMO DI VITERBI: 341 for i = 2:L 342 for l = 1:m 343 for k = 1:m 344 temp (k) = V(k,i-1) * a(k,l)); 345 end 346 [V(l,i), tr (l,i)] - max (temp ); 347 V(l,i) = V(l,i) * e(l,x) 348 end 349 end 350 % NB: passando a matlab non usiamo la produttoria perchè la moltiplicazione di molti numeri molto piccoli genera un risultato che potrebbe essere inferiore alla soglia di matlab. Per questo motivo trasformiamo la produttoria come logaritmo poichè il logaritmo del prodotto equivale alla somma di logaritmi --> log(a*b) = log(a) + lob(b). 351 352 % Ho una sequenza di 2000 nucleotidi e voglio sapere se ci sono delle CpG island. 353 354 load ('ex2' ) 355 load ('state_transition_matrix' ) 356 load ('emission_probabilities_matrix' ) 357 358 SYMBOLS ='ACGT' ; % Definisco l'alfabeto 359 E = max (eps ,E); % eps è un numero piccolo generico, dobbiamo sostituire tutti gli zeri con eps perchè log(0) non esiste. Sostituendo gli zeri con eps (molto piccolo) abbiamo comunque numeri molto piccoli che verranno scartati. 360 E = log2 (E); % sostituisco E con il suo logaritmo 361 362 [m,n] = size (E); 363 L = length (x); 364 % m = numero di stati 365 % n = dimensione alfabeto 366 % L = lunghezza sequenza 367 368 V = zeros (m,L); 369 tr = zeros (m,L); 370 temp = zeros (1,m); 371 path = zeros (1,L); 372 373 % Inizializzazione della prima colonna di V 374 for k = 1:m 375 v(k,1) = log2 (1/m) + E(k,strfind (SYMBOLS ,x(1))); % Questa strfind mi ritorna l'indice di A in ACGT 376 % log(a) = log(1/m) poichè la probabilità di trovare quella base in una posizione è esattamente 1 / numero di posizioni 377 end 378 379 % Iterate over sequence 380 for i = 2:L 381 for l = 1:m 382 for k = 1:m 383 temp (k) = V(k,i-1) + log2 (a(k,l)); 384 end 385 [V(l,i), tr (l,i)] = max (temp ); 386 V(l,i) = V(l,i) + E(l,strfind (SYMBOLS ,x(i))); 387 end 388 end 389 390 % Trovo lo score vincente 391 [maxscore ,path (L)] = max (V(:, L)); 392 393 % Traceback 394 for i = L:- 1:2 395 path (i-1) = tr (path (i), i); 396 end 397 398 figure 399 plot (path ) 400 xlabel ('bp' ); 401 ylabel ('state' ); 402 %sono CpG islands le posizioni tra circa 400 e 1000 perchè sono nello stato 403 %1-4 (quelle in basso) 404 405 ////////////////////////////////////////////////////////////////////// 406 407 %%%%%%%% LAB5 - low level analysis for microarray data 408 409 %% ESERCIZIO 1 410 411 % Import data corresponding to the two channels of a cDNA microarray experiment (microdata.mat) and visualize them. 412 % 1) derive signal-background intensities for the two channels 413 % 2) remove bad quality spots, deleting elements for which the flag value is -50 414 % 3) remove negative values 415 % 4) draw a log2(R) vs log2(G) plot 416 % 5) draw an MA plot 417 % 6) draw an MA plot using the Matlab function mairplot 418 % 7) apply a lowess normalization 419 420 load ('microdata' ) 421 % importo il file microdata.mat 422 % all'interno ci sono 5 variabili 423 % segnali: Gsig e Rsig 424 % flags: vettore che contiene un valore di qualità per ogni spot del microarray 425 % rumore di fondo: Rbkg e Gbkg 426 427 % 1) vado a rimuovere i rumori di fondo per ottenere i miei segnali puliti: 428 G = Gsig - Gbkg ; 429 R = Rsig - Rbkg ; % --> adesso contengono solo i segnali dovuti all'effettiva attività biologica. Però adesso qualche spot può essere diventato negativo (se il background era maggiore del segnale), quindi voglio rimuoverli. 430 431 % 3) Creo un vettore che contiene gli spot negativi e li unisco in un unico vettore: 432 neg_valueR = find (R uniforma le distribuzione di intensità per renderli confrontabili tra di loro. 500 % mi salvo in un vettore i nomi dei campioni, che mi serviranno dopo per il %plot: 501 samples = probeData .CELNames ; 502 503 % 3) voglio plottare 3 diversi grafici all'interno della stessa figura; lo faccio tramite la funzione subplot, che mi crea una griglia. Passo a subplot righe, colonne, numero della casella in cui metto il plot 504 figure 505 subplot (3,1,1) % 3 righe, 1 colonna, casella della griglia da cui vogliamo partire (1°) 506 507 % in input alla funzione boxplot passo i log2 dei valori di intensità grezzi e plotto soltanto i primi 4 campioni: 508 maboxplot (log2 (probeData .PMIntensities (:, 1:4)), samples (1:4), 'title' ,'Raw intensities' ,'orientation' ,'horizontal' ) 509 % input: valori, nomi, titolo, orientamento del boxplot 510 511 % vado avanti con la seconda casella 512 subplot (3,1,2) % voglio plottare le intensità dopo la rimozione del segnale di background. Prendo ancora solo i primi 4 campioni: 513 maboxplot (log2 (pms_bg (:, 1:4)), samples (1:4), 'title' ,'After background adjustment' , 'orientation' ,'horizontal' ) 514 515 % aggiungo la terza casella, in cui voglio plottare i dati finali, ovvero quelli ottenuti dopo la normalizzazione quantile: 516 subplot (3,1,3) 517 maboxplot (log2 (pms_bgnorm (:, 1:4)), samples (1:4), 'title' ,'After quantile normalization' ,'orientation' ,'horizontal' ) 518 519 % OSSERVAZIONE: all'inizio ho valori estremi, poi le distribuzioni diventano tutte più simili tra loro, e questo significa che posso confrontare i campioni tra loro, cosa che prima non sapevo fare. La normalizzazione quantile rende uguali le mediane e lo spazio interquartile. 520 521 % 4) ULTIMO STEP 522 cns_rma_exp = rmasummary (probeData .ProbeIndices , pms_bgnorm ); 523 % Le probe dello stesso probeset sono state accorpate quindi ora ho una riga per ogni probeset e non per ogni probe. 524 % Ora devo plottare i valori di intensità per i probset 3500 e 4221: 525 figure 526 subplot (2,1,1) %posizione 1 527 plot (cns_rma_exp (3500 ,:)) 528 subplot (2,1,2); %posizione 2 529 plot (cns_rma_exp (4221 ,:)) 530 531 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 532 533 %% ESERCIZIO 3 - differential analysis 534 535 % Consider a gene expression profiling experiment on breast cancer samples. Identify the upregulated genes (Bonferroni corrected pvalue 1.3 (concentrazione soglia). 732 % 2) Disegnare la curva ROC al variare della concentrazione soglia 733 734 load ('ex1.mat' ); 735 n_samples = length (disease ); 736 737 % cerco quali pazienti sono realmente malati/sani 738 POS = find (disease == 1); 739 NEG = find (disease == 0); 740 741 % Il classificatore è la concentrazione e la soglia è 1.3 --> per cercare i veri positivi, devo cercare all'interno di POSITIVI quelli che hanno concentrazione > 1.3 742 TP = length (find (concentration (POS ) > 1.3 )); 743 TN = length (find (concentration (NEG )