9 avril 2008 – recherche du maximum de densité des segments (2)

10h43
suite de recherche du maximum de densité des segments
je dois donc modifier mon code pour tenir compte du décalage de 18 nucléotides.
2 possibilités:
- je modifie les positions des tags sur chr1278
- je rajoute une étape dans le script de find_max_segment qui repositionne le segment comme par rapport au start des sages.

la 2eme solution est simple (la première aussi cela dit) mais elle implique que le génome utilisé sera toujours décalé de 18.
c’est finalement idiot, il serait mieux de partir avec un génome bien recalé.
quand j’étudie en détail les position je constate que:
- les positions sur les brins watson sont celles du premier nucléotide du sage (coté 5′)
- les positions sur les brins crick sont celles du dernier nucléotide avant le sage (en 5′)
ce qui nous intéresse c’est finalement la position du premier nucléotide après le sage, en 3′
je dois donc faire :
-> pos + 18 pour les tag sur les brin watson
-> pos – 19 pour les tag sur le brin crick

la matrice chr_len donne les positions de chaque chromosome et brin sur le génome linéarisé. a partir de cette matrice je peux effectuer les corrections nécessaires.
j’ai écrit un script pour ca

function res = tag_position_correction(genome, chr_len)%res = tag_position_correction(genome, chr_len)%fonction permettant de modifier la position des tag sur un génome linearisé%de +18 s'ils sont sur le brin watson et -19 sur le brin crick%la position des brin est donnée dans chr_len

%recupere la taille du genomenb_chr = size(chr_len,1);%alloue l'espace pour le resultatres = zeros(size(genome));%pour chaque birn de chaque chromosomefor i = 1:nb_chr %recupère l'information du brin strand = chr_len{i,2}; if strcmp('c', strand)     %si c'est crick alors decalage de -19     decalage = -19; else     %si c'est watson alors decalage de +18     decalage = 18; end %teste de la valeur de i dans la boucle if i>1     %si i>1 alors on recupère la position du début du chromosome sur la     %ligne précédente     %on ajoute +1 car la ligne précédente donne le dernier nucléotide     %du brin précédent     %on rajoute + 19 pour anticiper le décalage     deb = chr_len{i-1,4} + 1 + 19; else     %si i=1 alors il n'y a pas de ligne précédente, on demarre au début     %du chromosome

     deb=1 + 19; end %on récupère la position de fin et on retire 18 pour anticiper le %décalage fin = chr_len{i,4} - 18; %on copie dans res la partie du génome entre deb et fin %on les affecte au position deb + decalage à fin + decalage res(deb + decalage:fin + decalage) = genome(deb:fin);end

je l’appelle avec la commande chr1278_2 = tag_position_correction(chr1278, chr_len);

je peux maintenant réappliquer le script density sur cette matrice
density = densite_sage(chr1278_2, 20);
et j’applique le script find_max_segment dessus
res=find_max_segment(density, ratio2_5_100_1);
dans ratio2_5_100_1 j’avais également une erreur de position pour les segment des brin crick qui sont decalés de +1
j’ai créé segment_list a partir de ratio2_5_100_1 en corrigeant ce décalage.
es brins ont tous une référence (de 1 à 34)
les brins crick sont tous impaires
je cherche pour chaque segment si la référence du brin sur lequel il se trouve
est paire ou non

l=mod(res(:,6),2);

comme j’ai 1 pour les brin crick et 0 pour les brin watson, et que le décalage est de +1
je retire la matrice l aux colonnes contenant les positions des brins

segment_list(:,[2 3 7 8]) = segment_list(:,[2 3 7 8]) - l*ones(1,4);

puis je recherche le max de densité des segments

res=find_max_segment(density, segment_list);

le résultat semble correct.
toutefois quand les tags d’un meme segments sont éloignés d’une distance supérieur à la fenetre utilisée pour calculer la densité, le max de densité est mal placé.
il faut peut-etre utiliser 100 comme fenetre pour etre sur de couvrir les segments les plus grand.
density_100 = densite_sage(chr1278_2, 100);
res_100=find_max_segment(density_100, segment_list);
cette fois il ne semble plus y avoir de probleme

il faut donc exporter ce résultat.
le résultat est enregistré dans la matrice seg_list_dist_density_max
avant de l’exporter j’ajoute une colonne avec l’ancienne valeur de milieu de segment.
cela permettra de ne pas modifier les références des segments des brins crick.
je modifierai les références plus tard.
pour ajouter la valeur du milieu du segment je fais
seg_list_dist_density_max(:,end+1)=floor(mean(ratio2_5_100_1(:,7:8)'))';
je sauve ensuite la matrice en csv
save(strcat('seg_list_dist_density_max','.mat.csv'),'seg_list_dist_density_max', '-ascii','-tabs');
puis je l’exporte sur douanier rousseau
scp berthe:/home/gim/sage/total/stats/ratio2_5/max_segments/seg_list_dist_density_max.mat.csv /home/gim/sage/total/stats/ratio2_5/max_segments/

j’execute sur ce fichier les script de segment.py
en modifiant le nom du fichier d’entrée.
j’ai modifié aussi le script de matlab2csv pour qu’il ajoute dans la ligne de titre les colonnes density_max et mid

le fichier contenant les résultat brut en csv est ~/sage/total/stats/ratio2_5/max_segments/seg_list_dist_density_max.csv
le fichier contenant les distances avec les features environanantes est dans ~/sage/total/stats/ratio2_5/max_segments/segments_description_density_max.csv

je transfert ce fichier sur bambino
j’applique dessus la macro classification du fichier F:/users/christophe/search_features_A_and_B.xls a
je récupère le fichier F:/users/christophe/sements_description_density_max.xls
je l’exporte sur douanierousseau dans
~/sage/total/stats/ratio2_5/max_segments/segments_description_density_max.xls
j’élimine les colonnes avec du texte
je sauve dans ~/sage/total/stats/ratio2_5/max_segments/segments_distance_to_feature_using_density_max.csv
j’exporte sur berthemorisot
je l’ouvre dans matlab dans la matrice data_using_density_max
j’ouvre data_using_max et data_using_mid dans distances_using_max.mat
je compare les résultats obtenus avec les autres matrices.
j’ai écrit un script qui affiche les courbes pour les 3 matrices en fn de la classe, et de la config des features

function plot_with_n_variable(classe, feature_extreme, feature_sens, data_using_max, data_using_density_max, data_using_mid)  cols=cell(4,3);

  cols{1,1}='start';  cols{1,2}='sens';  cols{1,3}=[10 16];

  cols{2,1}='start';  cols{2,2}='antisens';  cols{2,3}=[12 18];

  cols{3,1}='end';  cols{3,2}='sens';  cols{3,3}=[9 17];

  cols{4,1}='end';  cols{4,2}='antisens';  cols{4,3}=[15 19];

  for i=1:4      if strcmp(cols{i,1}, feature_extreme) && strcmp(cols{i,2}, feature_sens)          columns=cols{i,3};      end  end

  figure;  title(strcat('distance from MAX of segment to ',feature_extreme, '  of features in  ',   feature_sens));  hold all;  for i=1:10  distance_c1_start_sens=mk_distance_hist(data_using_max, classe, columns, -1000, 1000, i, 20);  hist_c1_start_sens=frequence(distance_c1_start_sens,-1000,1000,200);  plot(hist_c1_start_sens(:,1),hist_c1_start_sens(:,2))  end

  figure;  title(strcat('distance from DENSITY MAX of segment to ',feature_extreme, ' of features in ', feature_sens));  hold all;  for i=1:10  distance_c1_start_sens=mk_distance_hist(data_using_density_max, classe, columns, -1000, 1000, i, 20);  hist_c1_start_sens=frequence(distance_c1_start_sens,-1000,1000,200);  plot(hist_c1_start_sens(:,1),hist_c1_start_sens(:,2))  end

  figure;  title(strcat('distance from MID of segment to ',feature_extreme, ' of features in ', feature_sens));  hold all;  for i=1:10  distance_c1_start_sens=mk_distance_hist(data_using_mid, classe, columns, -1000, 1000, i, 20);  hist_c1_start_sens=frequence(distance_c1_start_sens,-1000,1000,200);  plot(hist_c1_start_sens(:,1),hist_c1_start_sens(:,2))  end

je peux alors utiliser la commande
plot_with_n_variable(1, 'start', 'sens', data_using_max, data_using_density_max, data_using_mid)




on observe une dispersion légèrement plus faible pour les résultats obtenus en utilisant le maximum de densité des segments.