Biomedical Engineering - Bioinformatica e Genomica Funzionale

Matlab's codes for the exam

Other

%% LAB 2 - esercizio 1 load alignment.mat A2=A; for j=1:60 A2=[A2; seqrcomplement(A2(j,:))]; %aggiunge la sequenze complementari end %inizializzo la weight matrix wm=zeros(4,21); %21=posizioni, 4=basi %per ognuna delle colonne leggo quante volte occore una base, %e così riempio la wm con le occorrenze 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 %adesso voglio trovare le frequenze: wm=wm/120; %sostituisco 0 con 1 nella matrice per non avere problemi col log wm(wm==0)=1; %% calcolo entropia di shannon H=-wm.*log2(wm); %sommo i valori delle 4 basi H=sum(H); %calcolo il contenuto informativo I %non considero l'errore I=2-H sol=find(I>0.8); %mi da le posizioni in cui I è maggiore di 0.8 stem(I); title('stem graphic'); seqlogo(wm); %% LAB 2 - esercizio 2 %per l'esecuzione di questo esercizio bisogna tenere valide tutte le %variabili dell'esercizio 1, è come fosse un unico script promoter=fastaread('K02672.fna'); %promoter è adesso la sequenza tutta intera %voglio 500 basi totali, 400 all'inizio e 100 alla fine rispetto alla ba\ se %3395 nel file l=length(promoter.Sequence); %% sovrascrivo promoter=promoter.Sequence(3395-400:3395+99); %creo il vettore segnale e lo inizializzo S=zeros(1,500-21+1); %ha queste dimensioni perchè non possiamo calcolare il segnale fino %all'ultima base, perchè sto procedendo a finestre di 21 caratteri. %% for j=1:500-21+1 %sposto la finestra for p=1:21 %scorro i caratteri della finestra switch promoter(j+p-1) case 'A' %è come scrivere if promoter(j+p-1)=='A' S(j)=S(j)+I(p)*wm(1,p); case 'C' S(j)=S(j)+I(p)*wm(2,p); case 'G' S(j)=S(j)+I(p)*wm(3,p); case 'T' S(j)=S(j)+I(p)*wm(4,p); end end end %richiesta esercizio: posizioni con un segnale >4.5 length(find(S>4.5)); find(S>4.5); %in queste posizioni la sequenza è maggiormente simile al sito di lega\ me %posso plottare il segnale X=1:500-21+1; plot(X,S); hold on plot(X,4.5*ones(size(X))); %per vedere la soglia %ora in effetti vedo 4 picchi che superano la linea di soglia %% LAB 3 - ESERCIZIO 1 %aprire ex1 e state transition matrix %split transition matrix load('ex1.mat') load('state_transition_matrix.mat') Trans_matrix_IN_CpG = a(1:4,1:4); %quadrante alto sinistra Trans_matrix_OUT_CpG = a(5:8,5:8); %quadrante basso destra %preparo la variabile in cui metto la probabilità cumulativa %ovvero la produttoria del modello markoviano %la inizializzo a 1 perchè è un prodotto P_sequence_IN_CpG = 1; P_sequence_OUT_CpG = 1; L=length(CpGtest); sequence=upper(CpGtest); %setto lo stato di partenza per tutte e due %set starting state switch (sequence(1)) case 'A' k=1; case 'C' k=2; case 'G' k=3; case 'T' k=4; end %con lo switch guardo il primo simbolo della sequenza %poi itero su tutta la sequenza %read sequence position by position for pos=2:L current_symbol=sequence(pos); switch(current_symbol) case 'A' l=1; case 'C' l=2; case 'G' l=3; case 'T' l=4; end trans_prob_IN_CpG = Trans_matrix_IN_CpG(k,l); %k stato precedente, l stato nuovo trans_prob_OUT_CpG = Trans_matrix_OUT_CpG(k,l); P_sequence_IN_CpG = P_sequence_IN_CpG*trans_prob_IN_CpG; P_sequence_OUT_CpG = P_sequence_OUT_CpG*trans_prob_OUT_CpG; k=l; end %result check if P_sequence_IN_CpG > P_sequence_OUT_CpG IN=1; else IN=0; end %mi esce IN=1 %% LAB 3 - ESERCIZIO 2 %ho una sequenza di 2000 nucleotidi e voglio sapere se ci sono delle %CpG island. %definisco l'alfabeto SYMBOLS='ACGT'; E = max(eps,E);%eps è un numero piccolo generico E = log2(E); % m = number of states [m,n] = size(E); L=length(x); V = zeros(m,L); tr = zeros(m,L); temp = zeros(1,m); path = zeros(1,L); %inizializzazione della prima colonna di V for k=1:m v(k,1) = log(1/m) + E(k,strfind(SYMBOLS,x(1))); %questa strfind mi ritorna l'indice di A in ACGT end %iterate over sequence for i=2:L for l=1:m for k=1:m temp(k)=V(k,i-1) + log2(a(k,l)); end [V(l,i),tr(l,i)]=max(temp); V(l,i)=V(l,i) + E(l,strfind(SYMBOLS,x(i))); end end temp=V(:,L); [maxscore,path(L)]=max(temp); %traceback for i=L:-1:2 path(i-1) = tr(path(i),i); end figure plot(path) xlabel('bp'); ylabel('state'); %sono cpg islands le posizioni tra circa 400 e 1000 perchè sono nello\ stato %1-4 (quelle in basso) %% LAB 4 - ESERCIZIO 1 load primates distances = seqpdist(primates, 'method', 'Jukes-Cantor', 'Alpha', 'DNA'); %% 1. costruisco l'albero UPGMAtree = seqlinkage(distances, 'UPGMA', primates); %eseguo e vedo che nel workspace appare UPGMAtree %% posso plottare l'albero h=plot(UPGMAtree, 'orient','top'); title('UPGMA Distance Tree of Primates using Jukes-Cantor model'); ylabel('Evolutionary distance'); %% 2. FILTRARLO ? impongo una soglia di distanza %voglio solo le foglie vicine a european human %salvo i nomi delle foglie dell'albero nella variabile "names" names = get(UPGMAtree, 'LeafNames'); %selezione delle foglie da tenere, salvo 2 vettori %3 è la foglia relativa a european human %come criterio per la selezione delle foglie deve utilizzare la distanza\ %rispetto alla foglia di reference %la soglia (treshold) di distanza è 0.3 [h_all,h_leaves] = select(UPGMAtree, 'reference' ,3, 'criteria', 'distance',... 'threshold', 0.3); subtree_names = names (h_leaves); %h_leaves contiene valori 0 e 1; gli 1 sono quelli che seleziono. %% adesso devo fare un vettore con le foglie che voglio RIMUOVERE DALL'A\ LBERO %utilizzo la funzione prune %per fare l'inverso uso la tilde ~ leaves_to_prune = ~h_leaves %definisco il nuovo albero pruned_tree = prune(UPGMAtree,leaves_to_prune); %plotto questo nuovo albero ancora con orientamento verticale h = plot(pruned_tree, 'orient', 'top'); title('albero filtrato') view(UPGMAtree, h_leaves); %vedo in rosso queste foglie sull'albero totale %% LAB 4 - ESERCIZIO 2 %DEFINISCO LA VARIABILE CHE CONTIENE IL NUMERO DI ITERAZIONI DA FARE num_boots = 100; %DEFINISCO ALTRI VALORI UTILI: %numero di sequenze: num_seqs = length(primates); %lunghezza dell'allineamento: seq_len = length(primates(1).Sequence); %ho preso la prima sequenza tanto sono tutte uguali %% 1. DEVO CREARE I NUOVI ALLINEAMENTI %creo l'arrey che li contenga boots = cell(num_boots,1); %ITERO SULLE 100 ITERAZIONI %ITERO SULLE 12 SEQUENCE for n = 1:num_boots %definisco il nuovo indice posizione reorder_index = randsample(seq_len, seq_len, true); %metto true perchè voglio che ci sia rimpiazzo for i=num_seqs:-1:1 bootseq(i).Header = primates(i).Header; bootseq(i).Sequence = strrep(primates(i).Sequence(reorder_index),\ ... '-', ''); end %salvo questa struct nell'array boots{n} = bootseq; end %% DEFINISCO UN NUOVO ARRAY %è un array di celle con 100 righe e una colonna boot_trees = cell(num_boots,1); for n = 1:num_boots %vettore temporaneo distanze dist_tmp = seqpdist(boots{n},'Method', 'Jukes-Cantor', 'Alpha', 'DNA'); boot_trees{n} = seqlinkage(dist_tmp, 'UPGMA', boots{n}); end for i = num_seqs-1:-1:1 %variabile che contiene il numero del nodo che voglio usare come radice branch_pointer = i + num_seqs; %ricavo il sottoalbero partendo da branch pointer come radice sub_tree = subtree(UPGMAtree,branch_pointer); %salvo i sottoalberi in un array di celle orig_pointers{i} = getcanonical(sub_tree); %getcanonical mi dà la forma canonica %salvo le specie contenute in questo sottoalbero in ordine alfabetico orig_species{i} = sort(get(sub_tree, 'LeafNames')); %sort ordina in ordine alfabetico end %itero su tutti i 100 alberi for j = num_boots:-1:1 for i = num_seqs-1:-1:1 branch_ptr = i + num_seqs; sub_tree = subtree(boot_trees{j}, branch_ptr); cluster_pointers{i,j} = getcanonical(sub_tree); cluster_species{i,j} = sort(get(sub_tree, 'LeafNames')); end end %% contare quante volte ogni ramo dell'albero originale ricorre nell'alb\ ero di bootstrap count = zeros(num_seqs-1,1); for i = 1:num_seqs-1 for j = 1:num_boots*(num_seqs-1) %numero totale sottoalberi %il primo if controlla che la topologia sia uguale if isequal(orig_pointers{i},cluster_pointers{j}) if isequal(orig_species{i},cluster_species{j}) count(i) = count(i) + 1; end end end end %% RICAVO LA CONFIDENZA ? CALCOLO LE FREQUENZE Pc = count./num_boots; %% 2. SELEZIONARE RAMI CON CONFIDENZA >0.9 [ptrs, dist, names] = get(UPGMAtree, 'POINTERS', 'DISTANCES', 'NODENAMES'); % AGGIORNO IL VALORE DEI NODI AGGIUNGENDO LA CONFIDENZA for i = 1:num_seqs-1 branch_ptr = i + num_seqs; names{branch_ptr} = [names{branch_ptr} ', confidence: ' num2str(100*Pc(i)) ' %']; end %RIDEFINISCO L'ALBERO USANDO QUESTI NUOVI NOMI tr = phytree(ptrs,dist,names); %ottengo un albero equivalente all'originale ma con names aggiornato con\ la %confidenza del ramo rappresentato dal quel nodo. %SELEZIONO I RAMI CON CONFIDENCE>0.9 %per ottenere l'indice corretto devo aggiungere numseq high_conf_branch_ptr = find(Pc>0.9 + num_seqs); %visualizzo il nuovo albero coi nodi in evidenza view(tr, high_conf_branch_ptr) %% LAB5 - ESERCIZIO 1 - ESPRESSIONE GENICA (MICROARRAY) % importo il file microdata.mat % all'interno ci sono 5 variabili % segnali: Gsig e Rsig % flags: vettore che contiene un valore di qualità per ogni spot del % microarray % rumore di fondo: Rbkg e Gbkg %vado a rimuovere i rumori di fondo per ottenere i miei segnali puliti: G = Gsig - Gbkg; R = Rsig - Rbkg; %adesso contengono solo i segnali dovuti all'effettiva attività biolo\ gica. %però adesso qualche spot può essere diventato negativo (se il background era maggiore del segnale %quindi voglio rimuoverli. Creo un vettore che contiene gli spot negativ\ i: neg_valueR = find(R=t1(length(t1))))]; nt=N*ones(1,length(t1)); for j=2:length(t1) nt(j)=nt(j-1)-eventi1(j-1)-c1(j-1); end %adesso ho tutto per calcolare la probabilità %devo cancolare innanzitutto la probabilità di ogni tempo f=eventi1./nt %s è la prob per ogni singolo intervallo, non la prob condizionata s=1-f; preKM=s(1); %inizializzo la produttoria KM=1; %KM valore della curva per ogni istante; lo inizializzo a 1 for k=1:length(t1)-1 %aggiorno la produttoria preKM=preKM*s(k+1); KM=[KM preKM]; end %% faccio il grafico figure stairs(t1, KM) if max(EVENT)>=max(censored) t1=[t1 max(EVENT)+1]; else %devo allungare l'asse x per includere anche l'esclusione avvenuta fdopo\ l'ultimo evebto di morte t1=[t1 max(censored)+1] %+1 solo perchè viene più ordinato il grafico end %% riplotto la curva col nuovo asse figure stairs(t1,[KM KM(end)]) %x e y devono avere la stessa dimensione quindi ho aggiungo un elemento anche a KM hold on %aggiungo i simboli che corrispondono alle esclusioni for i=1:length(censored) %se l'evento cade nel mezzo di un intervallo metto la croce T=find(t1