15 mai 2008 – analyse des résultat de sage – restriction WT et double mutant

10h25
Je vais reprendre l’analyse des sage a la base en ne prennant en compte que les tag obtenu sur les manip wild type et double mutant Trf4/Rrp6

Je dois refaire le clustering afin de n’avoir que des segment correspondant a des groupes de tag de ces manip.

1- ouverture dans matlab du fichier
~/sage/total/stats/tag_tot_genome.mat

j’ai dans ce fichier les matrice correspondant aux 12 manip de sage (header01 à header12)
chaque matrice à une taille de 24 million ce qui corresponde au nombre de nucleotide de l’ensemble du génome mis bout a bout en séparant les brin C et W
les matrice chr_len et chr_ref donnent les positions de chaque chromosome dans cette matrice

2- la recherche des segments (cluster) va se faire uniquement sur les manip 1 et 8, correspondant au Wild Type (WT) et au double mutant Trf4/Rrp6

3- j’utilise le programme segment.m avec la commande
header = header01 + header08;
header(:,2) = header01;
[header_region] = segment(header, dt, nt, chr_ref);
je récupère la matrice ratio2_1_100_1 qui a été crée
et je représente les histogramme des classe de ratio

figurefor i=1:9subplot(3,3,i)hist(ratio2_1_100_1(find(ratio2_1_100_1(:,4)>i-1),1))xlabel('Ratio Oligo-dT / (CUT + Oligo-dT)', 'color',[1 1 1]);ylabel('Number of clusters', 'color',[1 1 1]);title(strcat('n = ', num2str(i)), 'color',[1 1 1])set(gca, 'Xcolor', [1 1 1], 'Ycolor', [1 1 1])end

et j’obtiens la figure suivante

on voit qu’il y a un nb plus important de cluster de la classe 1 que ceux de la classe 10
si on regarde le nb de tag, on constate aussi qu’il y a environ 20 % de plus de CUTSque d’Oligo-dT

On va donc chercher a réduire artificiellement le nb de CUT en en éliminant de façon aléatoire jusqu’à avoir le meme nombre
script : delete_rand_tag

h8 = header08;tag = find(h8);s = size(tag,1);while sum(h8) > 49284p = floor(rand(1000,1) * s);ptag = tag(p);h8(ptag) = 0;end

header_rand = header01 + h8;
header_rand(:,2) = header01;
je fais ensuite [header_region_rand] = segment(header_rand, dt, nt, chr_ref);
je renome ratio2_1_100_1 en ratio2_1_100_1_ rand
puis je représente les memes histogrammes

figurefor i=1:9subplot(3,3,i)hist(ratio2_1_100_1_rand(find(ratio2_1_100_1_rand(:,4)>i-1),1))xlabel('Ratio Oligo-dT / (CUT + Oligo-dT)', 'color',[1 1 1]);ylabel('Number of clusters', 'color',[1 1 1]);title(strcat('n = ', num2str(i)), 'color',[1 1 1])set(gca, 'Xcolor', [1 1 1], 'Ycolor', [1 1 1])end

et j’obtiens

on voit que c’est plus homogène entre les classes extrèmes.
pour vérifier cela, on cherche a voir si la répartition du nombre de tag dans les segment est la meme pour toute les classes.

j’extrait de la matrice ratio2_1_100_1_ rand les lignes pour lesquelles le ratio est faible (inf a 0.2 CUTs) et fort (sup à 0.8 WT) et je fait l’histogramme du nombre de tag (colonne 4). je choisi pour ces histogramme d’avaoir comme nombre de catégorie la valeur maximale de tag pour chaque type (3318 pour le WT et 1522 pour les CUTs.
puis je trace les courbe correspondantes en échelle log

h_0=hist(ratio2_1_100_1_rand(find(ratio2_1_100_1_rand(:,1)0.8),4), max(ratio2_1_100_1_rand(find(ratio2_1_100_1_rand(:,1)>0.8),4)));
figure;
plot(h_0,'r');
hold on
plot(h_1,'g');

on voit que les segments avec un faible nb de tag sont plus nombreux pour les CUTs. cependant a partir de 4 tags par segment on observe la meme répartition.
on peut donc dire que a partir de 4 tags par segment on a un nombre similaire de segments avec le meme nombre de tag dans ceux-ci.

Je me suis rendu compte que les position des tags sont fausses. en fait elles sont décalée par rapport a ce qu’on veux car on s’interresse a la fin des tag et non au début
j’ai fait un programme qui corrige ca (tag_position_correction)
je corrige donc les 2 matrice contenant les résultats des WT et CUTs

chr_WT=tag_position_correction(header01,chr_len);chr_CUTs=tag_position_correction(h8,chr_len);chr_WT_CUTs = chr_WT + chr_CUTs;chr_WT_CUTs(:,2)=chr_WT;[regionWT_CUTs] = segment(chr_WT_CUTs, dt, nt, chr_ref);

cpt =

       12021

et je récupère la matrice segments_WT_CUTs
je vais maintenant rechercher le max de densité de chaque segment
je calcul d’abord la matrice de densité pour les tag WT et CUts
density_WT_CUTs = densite_sage(chr_WT_CUTs(:,1), 100);
puis je cherche le max de densité
segments_WT_CUTs=find_max_segment(density_WT_CUTs, segments_WT_CUTs);
j’exporte avec save(strcat('segments_WT_CUTs','.mat.csv'),'segments_WT_CUTs', '-ascii','-tabs');

TODO -> export sur douanierousseau, localiser les segment par rapport aux orf et faire les graph

8 avril 2008 – recherche du maximum de densité des segments

17h52 :
au lieu de localiser les segments par rapport à la position du maximum de tag dans un segment, je vais le faire par rapport à la position du maximum de densité.
pour evaluer la densité, je fais avancer une fenetre dont la taille L est a déterminer.
je compter alors le nb de sage dans la fenetre et cette valeur sera alors la densité pour la position correspondant au milieu de la fenetre.
ainsi si plusieurs tag se suivent mais sans se chevaucher, la fenetre de densité va augmenter qd meme en s’approchant du milieu du segment car c’est a cet endroit qu’il y aura un max de segment environnant.



le schema ci-dessus représente un comptage théorique pour un segments dont aucun tag ne se chevauchent. les points bleus représentent la courbe de densité qui serait obtenu. on voit que le max de densité serait environ au 2/3 du segment.
j’ai écrit une fonction pour le calcul de la densité

function density = densite_sage(chr, window)%density = densite_sage(chr, window)%-input%  chr -> entire genome%  window -> size of the window to determine the densitty%-output%  density -> result, matrice with the size of chr

%recupere la taille du genometaille_genome=length(chr);%alloue l'espace mémoire necessaire pour le resultatdensity=zeros(taille_genome,1);%calcul la taille de demi-fenetre%je ne garde que la partie entièrew=floor(window/2);%calcul la somme de tout les tag pour une position i sur une fenetre%encandrant cette position%la fenetre est de taille window si window est impaire et window +1 si%window est pairefor i=1+w:taille_genome-w    density(i)=sum(chr(i-w:i+w));end

le calcul est très simple et renvoie une matrice de densité
j’ouvre le fichier /data/users_data/helen/454/sageData/total/stats/ratio2_5/max_segments/matlab.mat
puis je créé la matrice chr1278 qui contient tous les tags des manip 1, 2 7 et 8 réparti sur l’ensemble du génome.
chr1278=header(:,5)
puis j’execute la commande
density = densite_sage(chr1278, 20);
j’ai donc ma matrice densité.
il faut maintenant associer a chaque segment son max de densité.
j’ai une matrice de segment ratio2_2_100_1
chaque ligne correspond a un segment dont j’ai les coordonnées dans les 2eme et 3eme colonnes.
je vais donc parcourir cette matrice
récupérer les coordonnées du segment sur le génome linéarisé (chr1278)
puis regarder dans cette section sur density ou se trouve le max
s’il y a plusieurs positions max, ca sera la plus centrale qui sera prise.
la position trouvée est inscrite dans la matrice des segment dans la dernière colonne

j’ai écrit un script pour ca

function segment_data = find_max_segment(chr_data, segment_data)%fn pour trouver dans chr_data le max pour chaque segment de segment_data%si chr_data est la liste des tag réparti sur le génome alors la fonction%renvoie la position du nb max de tag%si chr_data est la densité de sage sur tout le génome, alors la fonction%renvoie la position de la densité max de tag dans le segment   [l,c] = size(segment_data);   for i=1:l      start_seg = segment_data(i,2);      stop_seg = segment_data(i,3);      segment=chr_data(start_seg:stop_seg,1);      max_seg=find(segment==max(segment));      pos = floor(mean(length(max_seg)));      segment_data(i,c+1)=segment_data(i, 7) + max_seg(pos) - 1;   end

script très simple.
je le teste par la commande res=find_max_segment(density, ratio2_5_100_1);
le programme marche bien mais en regardant les résultats j’ai constaté que les valeurs étaient assez étrange et souvent mal placée.
en regardant de plus pret j’ai compris l’erreur que j’avais faite.
les positions des tag sur le génnome correspondent au debut des sage sur les chromosome alors que les position des segments dans la liste des segment correspondent à la première position après la fin du sage, soit 18 nucleotides plus loin. la recherche du max se fait donc avec un decalage de 18, mais pas toujours dans le meme sens selon que la partie du génome considérée soit sur le brin watson ou crick.
je vais devoir corriger cela