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 commandedensity = 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
