13h23 :
en reprenant le code VBA pour classer les segments, j’ai constaté quelques erreurs.
je l’ai donc repris a la base.
je l’ai finalement totalement changer.
après quelques tatonnement j’ai abouti a un code correct.
mon problème principal a été l’orientation des segments et des orf.
en effet, pour respecter le sens 5′->3′ de l’adn, selon le brin sur lequel se trouve une séquence, la valeur de la position de début d’une séquence sera supérieure a celle de la fin si on se trouve sur le brin crick et inférieure pour le brin watson.
cependant je n’ai pas respecté cette convention d’ecriture pour les segments mais uniquement pour les orf.
c’est un peu dommage et surtout ca pose des problème lors des calculs.
j’ai donc choisi d’inverser aussi les valeurs pour les orf lors du calcul, pour me simplifier la vie.
-description de l’algo de classification
pour chaque segments trouvé
je récupère le start et le end. attention les start est focement inférieur ou égal au end. le brin n’est pas pris en compte
je cherche d’abord l’orf A. par définition elle est forcement après le segment et ne le chevauche pas. je la cherche donc parmis les orf upstream en sens et antisens.
je cherche ensuite l’orf B. celle-ci peut etre chevauchante ou recouvrir en intégralité le segment. elle se trouve parmi les orf downstream et overlap en sens et antisens. la plus proche sera celle choisie. si il y a une orf overlap en sens et en antisens, celle en sens sera privilégiée.
on regarde ensuite la disposition du segment par rapport a ces 2 orfs
- si les 2 orf et le segment sont dans le meme sens c’est une configuration tandem sens
- si les 2 orf sont dans le meme sens et le segment en antisens c’est une configuration tandem antisens
- si l’orf A est dans le meme sens que le segment et B en antisens c’est une configuration convergente
- si l’orf A est en antisens du segment et que B est dans le meme sens, c’est une configuration divergente
on recherche ensuite ou se situe le segment par rapport a l’orf B, soit avant, soit chevauchant le début, soit inclu soit chevauchant la fin.
pour sa on caclul la différence entre les positions start et end du segment et de l’orf.
au préalable, si l’orf est sur le brin crick, on aura inversé les valeurs de start et de end.
on teste si ces différences sont positives ou négatives.
si le segment est devant l’orf, toute les différence seront positives par exemple.
dans le cas ou le segment est sur le brin crick, les résultat sont inversé. les différence sont alors toute négative si le segment est avant l’orf. pour prendre en compte cette configuration je multiplie par -1 les différence si le segment est sur crick et par 1 si il est sur watson.
si j’associe 1 à une valeur de différence positive et zero a une valeur négative, et que ensuite je fait la somme, j’obtiens alors un chiffre compris entre 0 et 4.
- 4 signifie que le segment est intégralement devant l’orf (cas 1, 5, 9, et 13)
- 3 signifie que le segment chevauche le début de l’orf (cas 2, 6, 10 et 14)
- 2 signifie que le segment est inclus intégralement dans l’orf (cas 3, 7, 11 et 15)ou que il recouvre l’orf intégralement. ce dernier cas n’avait pas été envisagé. je ne le traite pas différemment car il est rare, sno rna en général
- 1 signifie que le segment chevauche la fin de l’orf (cas 4, 8, 12 et 16)
- 0 signifie que le segment est intégralement derrière l’orf, ce qui n’est pas possible a priori selon la méthode de sélection des orf environnantes
le code vba
Attribute VB_Name = "search_features"Sub classification()
Range("AH2").Selectligne = 0Valeur = Range("A2").Offset(ligne, 0).Value
Do Until Valeur = ""'recupération des start et end du segmentsegment_start = ActiveCell.Offset(ligne, -30)segment_end = ActiveCell.Offset(ligne, -29)'calcul du milieu du segment'segment_mid = Int((segment_start + segment_end) / 2)'récupération de la position du max du segmentsegment_max = ActiveCell.Offset(ligne, -28)'calcul des distances entre le start et le milieu'segment_start2mid = segment_mid - segment_start'calcul des distances entre le end et le milieu'segment_end2mid = segment_end - segment_mid'calcul des distances entre le start et le max'segment_start2max = segment_max - segment_start'calcul des distances entre le end et le max'segment_end2max = segment_end - segment_max
'teste si le brin est crick ou watson pour initialiser la variable signe'cette variable servira par la suite pour des calculs de distanceIf Range("C2").Offset(ligne, 0) = "c" Then'si crick on affecte a signe la valeur -1signe = -1Else'sinon on lui donne la valeur 1signe = 1End If
'********************recherche de l'orf A*********************************'l'orf A est forcement l'orf upstream la plus proche, elle est definit comme ca'la recherche se fait donc entre les feature up en sens et antisens
'appel de text_prox. les valeurs en argument sont les distances des features up en sens et antisens'test_prox retourne vrai si la première (feature_up_sens) est la plus procheup_val = test_prox(ActiveCell.Offset(ligne, -21).Value, ActiveCell.Offset(ligne, -14).Value)'teste de up_valIf up_val Then's'il est vrai, alors l'orf A est l'orf up senscol_orf_A = -24'et une 2eme variable de signe est initialisée a 1signe2 = 1Else's'il est faux alors l'orf A est l'orf up antisenscol_orf_A = -16'et la variable signe2 est initialisée a -1signe2 = -1End If'je recupère les noms de l'orf AORF_A_name = ActiveCell.Offset(ligne, col_orf_A).ValueORF_A_gene = ActiveCell.Offset(ligne, col_orf_A + 1).Value'je calcul la position du start de l'orf a partir du max du segment et de la distance a celui ci'les varaible de signe permette de tenir compte du brin'et de la position relative de la feature par rapport au segmentORF_A_start = segment_max - (ActiveCell.Offset(ligne, col_orf_A + 2).Value * signe * signe2)ORF_A_end = segment_max - (ActiveCell.Offset(ligne, col_orf_A + 3).Value * signe * signe2)
'********************recherche de l'orf B*********************************'l'orf B va etre l'orf la plus proche parmis les orf restante (downstream et overlap)'l'algo cherche d'abord s'il existe des orf overlap et sinon recherche entre les orf down laquelle est la plus proche
'teste s'il y a une orf overlap en sensIf Not IsEmpty(ActiveCell.Offset(ligne, -6)) Then'si oui, la variable down_val est vrai, cette variable signifie que l'orf B est sur le meme brin que le segmentdown_val = True'l'orf B sera alors l'orf qui overlap en senscol_orf_B = -8'recupération des distances du segment a l'orfmax2orf_start = ActiveCell.Offset(ligne, -6).Valuemax2orf_end = ActiveCell.Offset(ligne, -5).Value'affectation du signe + a la variable signe2 qui signifie que B et le segment sont dans le meme sensesigne2 = 1'teste s'il y a une orf overlap en antisens -- s'il ya des orf overlap en sens et antisens, celle en sens sera la seule considéréeElseIf Not IsEmpty(ActiveCell.Offset(ligne, -2)) Then'la variable down_val est alors faussedown_val = False'l'orf B est l'orf overlap antisenscol_orf_B = -4'recupération des distances du segment a l'orf, attention il y a une inversion entre start et end pour tenir compte de l'antisensmax2orf_end = -ActiveCell.Offset(ligne, -2).Valuemax2orf_start = -ActiveCell.Offset(ligne, -1).Value'signe - pour tenir compte de l'antisenssigne2 = -1'si il n'y a aucune orf overlapElse'appel de test_prox pour savoir quelle orf entre down sens et down antisens est la plus prochedown_val = test_prox(ActiveCell.Offset(ligne, -18).Value, ActiveCell.Offset(ligne, -9).Value)'test si down_val est vraiIf down_val Then'si oui, alors l'orf down sens est la plus proche'ce sera l'orf Bcol_orf_B = -20'on récupère les distances du segment a l'orfmax2orf_start = ActiveCell.Offset(ligne, -18).Valuemax2orf_end = ActiveCell.Offset(ligne, -17).Value'affectation du signe +signe2 = 1Else'sinon c'est la doown antisens la plus proche, ce sera l'orf Bcol_orf_B = -12'on recupère les distance, attention a l'inversion pour respecter l'antisensmax2orf_end = -ActiveCell.Offset(ligne, -10).Valuemax2orf_start = -ActiveCell.Offset(ligne, -9).Value'affectation du signe - pour l'antisenssigne2 = -1End IfEnd If'recup des noms de l'orf BORF_B_name = ActiveCell.Offset(ligne, col_orf_B).ValueORF_B_gene = ActiveCell.Offset(ligne, col_orf_B + 1).Value'calcul des start et end de l'orf B en utilisant les signes qui vont bienORF_B_start = segment_max - (ActiveCell.Offset(ligne, col_orf_B + 2).Value * signe * signe2)ORF_B_end = segment_max - (ActiveCell.Offset(ligne, col_orf_B + 3).Value * signe * signe2)
'********************classification des segments*********************************'la classif se fait en 2 étapes'on cherche d'abord la classe entre Tandem sens(Ts), Tandem antisens (Ta), convergente (C) et divergente (D)
'si les orf A et B sont en sens, alors c'est le cas tandem sensIf (up_val And down_val) Then'on inscrit TsActiveCell.Offset(ligne, 0).Value = "Ts"'la variable base_classe est initialisée a 0'cette variablme servira pour le calcul de la catégoriebase_classe = 0End If'si A et B sont en antisens alors c'est tandem antisensIf (Not up_val And Not down_val) ThenActiveCell.Offset(ligne, 0).Value = "Ta"base_classe = 4End If'si A est en antisens et B en sens alors c'est divergentIf (Not up_val And down_val) ThenActiveCell.Offset(ligne, 0).Value = "D"base_classe = 8End If'si A est en sens et B en antisens alors c'est convergentIf (up_val And Not down_val) ThenActiveCell.Offset(ligne, 0).Value = "C"base_classe = 12End If'appel de test_case pour déterminer la catégorie'recupère une valeur entre 1 et 4 qui est ajoutée a la variable base_classeclasse = test_case(segment_start, segment_end, ORF_B_start, ORF_B_end, signe, signe * signe2) + base_classe'ecrit les valeurs trouvées pour l'orf AActiveCell.Offset(ligne, 1).Value = ORF_A_nameActiveCell.Offset(ligne, 2).Value = ORF_A_geneActiveCell.Offset(ligne, 3).Value = ORF_A_startActiveCell.Offset(ligne, 4).Value = ORF_A_end
'ecrit les valeurs trouvées pour l'orf BActiveCell.Offset(ligne, 5).Value = ORF_B_nameActiveCell.Offset(ligne, 6).Value = ORF_B_geneActiveCell.Offset(ligne, 7).Value = ORF_B_startActiveCell.Offset(ligne, 8).Value = ORF_B_end'ecrit la categorieActiveCell.Offset(ligne, 9).Value = classe'incremente la ligneligne = ligne + 1'modifie la valeur dans la status barApplication.StatusBar = ligne'recupere la valeur de la permière cellule de la ligne suivante et continue la boucleValeur = Range("A1").Offset(ligne, 0).ValueLoop
End SubFunction test_prox(val_sens, val_anti)If Abs(val_sens) < Abs(val_anti) Thentest_prox = TrueEnd IfEnd Function
Function bool_2_val(bool)If bool Thenbool_2_val = 1Elsebool_2_val = 0End IfEnd Function
Function test_case(segment_start, segment_end, ORF_B_start2, ORF_B_end2, signe, brin_orf)'fonction test_case'permet de déterminer la catégorie du segment en fn de sa position'il a y 4 cas, entre les orf A et B(1), chevauchant le debut de B (2)'inclu dans B (3) et chevauchant la fin de B (4)'il y a aussi le cas 100 si jamais le cas n'est pas répertorié
'teste si l'orf est sur le brin crick ou watsonIf brin_orf < 0 Then'si c'est crick alors on iverse le start et le endtmp = ORF_B_end2ORF_B_end2 = ORF_B_start2ORF_B_start2 = tmpEnd If'teste si le start du segment est avant le start de l'orf (1=oui, 0=non)start_segment_b4_start_orfB = bool_2_val((ORF_B_start2 - segment_start) * signe > 0)'teste si le end du segment est avant le start de l'orf (1=oui, 0=non)end_segment_b4_start_orfB = bool_2_val((ORF_B_start2 - segment_end) * signe > 0)'teste si le start du segment est avant le end de l'orf (1=oui, 0=non)start_segment_b4_end_orfB = bool_2_val((ORF_B_end2 - segment_start) * signe >= 0)'teste si le end du segment est avant le end de l'orf (1=oui, 0=non)end_segment_b4_end_orfB = bool_2_val((ORF_B_end2 - segment_end) * signe >= 0)
'fait la somme de ces testssom_bool = start_segment_b4_start_orfB + end_segment_b4_start_orfB + start_segment_b4_end_orfB + end_segment_b4_end_orfB'si la somme est de 4, tout le segment est avant, on est dans le cas 1If som_bool = 4 Thentest_case = 1'si la some vaut 3, un bout du segment chevauche le debut de l'orfElseIf som_bool = 3 Thentest_case = 2'si la somme vaut 2, le segment est dans l'orf ou alors l'orf est dans le segmentElseIf som_bool = 2 Thentest_case = 3'si la somme vaut 1 alors le segment chevauche la fin de l'orfElseIf som_bool = 1 Or som_bool = 0 Thentest_case = 4Elsetest_case = 100End If
End Function
le script a été enregistré sur bambino dans F:/users/christophe/search_features_A_and_B.xls
il a été appliqué sur les données de F:/users/christophe/segments_description2_5_100_1.txt
j’ai remarqué que tous les segments du brin crick étaient mal positionné. ils sont décalé de +1. j’ai fait une correction des positions sur le start et end de ces segments avant d’appliquer le script. le résultat a été enregistré dans F:/users/christophe/segments_distance_to_feature_using_max.xls a
j’obtiens comme classification
| catégorie | mid to feature | max to feature |
| 1 | 2417 | 2423 |
| 2 | 196 | 189 |
| 3 | 1334 | 1347 |
| 4 | 22 | 17 |
| 5 | 1103 | 1110 |
| 6 | 240 | 241 |
| 7 | 2012 | 2032 |
| 8 | 33 | 27 |
| 9 | 1535 | 1539 |
| 10 | 218 | 222 |
| 11 | 1494 | 1522 |
| 12 | 24 | 21 |
| 13 | 1917 | 1926 |
| 14 | 378 | 379 |
| 15 | 1225 | 1243 |
| 16 | 24 | 21 |
| 100 | 11 | 0 |
| 104 | 21 | 0 |
| 108 | 32 | 0 |
| 112 | 23 | 0 |
c’est très proche de ce que j’obtenais avant.
je vais maintenant voir l’influence de la localisation des features par rapport au max des segments sur les distances obtenues
je génère sur matlab les graphique représentant la distribution.
pour cela j’exporte le fichier xls en format txt tabulé en ayant avant enlevé toutes les colonnes contenant du texte.
je transfert ce fichier sur berthemorisot /data/users_data/helen/454/sageData/total/stats/ratio2_5/max_segments/segments_distance_to_feature_using_max.csv
et je l’ouvre dans matlab
je l’enregistre dans la matrice data_using_max . les titre des colonnes sont dans textdata
je fais pareil pour les données calculée avec le milieu des segment qui sont enregistrée dan sla matrice data_using_mid.
j’enregistre le tout dans le fichier distance_using_max.mat
pour représenter les résultats graphiquement j’ai créé un script qui permet de filtrer les valeurs de distance.
function distance = mk_distance_hist(data, groups, cols, distance_min, distance_max, nb_tag_min, smooth_window)indice_groups = [];%selectionne les lignes faisant parti des groupes a utiliserfor g=groupsid = find(data(:,7)==g);indice_groups(end+1:end+length(id)) = id;end%conserve uniquement les lignes des groupes choisidata = data(indice_groups,:);%recherche les ligne avec un nb de tag superieur ou egal a la limiteindice_nb = find(data(:,5)<=nb_tag_min);%conserve uniquement ces lignes ainsi que les colonnes choisiesdata = data(indice_nb, cols);%recherche la valeur la plus proche de zero pour chque lignedistance=near_zero(data');distance=distance';%recherche les valeurs comprises dans l'intervalle spécifiédistance_id = find(distance<=distance_min & distance>=distance_max);distance=distance(distance_id);%cherche les extreme de la matricedist_max = max(distance);dist_min = min(distance);%applique la fonction histogramme en choisisant une catégorie par%valeurh=hist(distance, dist_max-dist_min + 1);%dessine la répartition des valeurs entre les extreme%figure;%plot(dist_min:dist_max, smooth_dist(h', smooth_window));%affiche le nb de segments aisni filtréslength(distance)
———————edit du 8 avril 2008———————-
on cherche a représenter les distances en fonction des classes, et de la configuration par rappoprt au orf.
par exemple, je veux voir la distance des segment par rapport au start des orf en sens.
on veux voir les segments proches des orf donc on va regarder les distance pour les orf downstream en sens mais aussi pour les orf overlap. c’est a dire les colonnes 10 et 16.
on s’intéresse aux segment de la classe 1.
on cherche pour tous les segments (n > 0)
et on cherche dans une fenetre réduite de -1000 a 1000
la commande est doncdistance_c1_start_sens=mk_distance_hist(data_using_max, 1, [10 16], -1000, 1000, 0, 20);
je recupère les valeurs dans la matrice distance_c1_start_sens
je vais appliquée une fonction permettant de calculer la fréquence dans une fenetre donnée avec un nombre de classes prédéterminé.
function histogramme = frequence(data, mini, maxi, nb) step = (maxi-mini) / nb;histogramme = [];for i = mini:step:maxideb = i;fin = i + step;histogramme(end+1,1) = deb;histogramme(end,2) = sum(logical(data >= deb & data < fin));endif (size(histogramme, 1) > nb)histogramme(end-1,2)= sum(histogramme(end-1:end,2));histogramme(end,:)=[];end
la commande est hist_c1_start_sens=frequence(distance_c1_start_sens,-1000,1000,200);
et j’obtiens
cette courbe varie en fonction du nombre de sage choisi
si on fait varier ce nombre de 1 à 10 on a
figure;hold all;for i=1:10distance_c1_start_sens=mk_distance_hist(data_using_max, 1, [10 16], -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
et on obtient
On peut aussi représenter sur un meme graphique les distances des segments des
différentes classes en sens et antisens
on prend les classes 1 et 10
on choisi d’abord de chercher les distance par rapport aux start des orf.
en sens il faut regarder par rapport au start des orf down et overlap en sens (colonnes 10 et 16) et en antisens il faut regarder par rapport aux orf upstream (on inverse) et overlap antisens (colonnes 12 et 18)
les commandes sont doncdistance_c1_start_sens=mk_distance_hist(data_using_max, 1, [10 16], -1000, 1000, 0, 20);distance_c1_start_antisens=mk_distance_hist(data_using_max, 1, [12 18], -1000, 1000, 0, 20);distance_c10_start_sens=mk_distance_hist(data_using_max, 10, [10 16], -1000, 1000, 0, 20);distance_c10_start_antisens=mk_distance_hist(data_using_max, 10, [12 18], -1000, 1000, 0, 20);
les histogrammes sont alors calculer parhist_c1_start_sens=frequence(distance_c1_start_sens,-1000,1000,200);hist_c1_start_antisens=frequence(distance_c1_start_antisens,-1000,1000,200);hist_c10_start_sens=frequence(distance_c10_start_sens,-1000,1000,200);hist_c10_start_antisens=frequence(distance_c10_start_antisens,-1000,1000,200);
puis on les représente sur un meme graphique
figurehold allplot(hist_c1_start_sens(:,1),hist_c1_start_sens(:,2))plot(hist_c10_start_sens(:,1),hist_c10_start_sens(:,2))plot(hist_c1_start_antisens(:,1),-hist_c1_start_antisens(:,2))plot(hist_c10_start_antisens(:,1),-hist_c10_start_antisens(:,2))
et on obtient la commande est 
j’ai écrit un script permettant de représenter les graph pour les start et les stop et de faire varier le nb de tag minimum
rappel
pour les distances par rapport au start on utilise
- les colonnes 10 et 16 en sens
- les colonnes 12 et 18 en antisens
pour les distances par rapport aux stop, on utilise
- les colonnes 9 et 17 en sens
- les colonnes 15 et 19 en antisens
la fonction
function plot_segment_distance(data, n)
distance = cell(8,3);distance{1,2}='distance from segment of classe 1 to start of orf in sens';distance{1,1}=mk_distance_hist(data, 1, [10 16], -1000, 1000, n, 20);distance{2,2}='distance from segment of classe 1 to start of orf in antisens';distance{2,1}=mk_distance_hist(data, 1, [12 18], -1000, 1000, n, 20);distance{3,2}='distance from segment of classe 10 to start of orf in sens';distance{3,1}=mk_distance_hist(data, 10, [10 16], -1000, 1000, n, 20);distance{4,2}='distance from segment of classe 10 to start of orf in antisens';distance{4,1}=mk_distance_hist(data, 10, [12 18], -1000, 1000, n, 20);distance{5,2}='distance from segment of classe 1 to end of orf in sens';distance{5,1}=mk_distance_hist(data, 1, [9 17], -1000, 1000, n, 20);distance{6,2}='distance from segment of classe 1 to end of orf in antisens';distance{6,1}=mk_distance_hist(data, 1, [15 19], -1000, 1000, n, 20);distance{7,2}='distance from segment of classe 10 to end of orf in sens';distance{7,1}=mk_distance_hist(data, 10, [9 17], -1000, 1000, n, 20);distance{8,2}='distance from segment of classe 10 to end of orf in antisens';distance{8,1}=mk_distance_hist(data, 10, [15 19], -1000, 1000, n, 20);
titles=cell(2,1);titles{1,1}=strcat('Distance from segment to start of features for n = ',num2str(n));titles{2,1}=strcat('Distance from segment to end of features for n = ',num2str(n));for i=1:8 distance{i,3}=frequence(distance{i,1},-1000,1000,200);endfor j=0:1 figure; xlabel('Distance (bp)'); ylabel('Frequency'); title(titles{j+1,1}); hold all for i=1:4 if even(i) signe = -1; else signe = 1; end id=(4*j) + i; plot(distance{id,3}(:,1),signe * distance{id,3}(:,2)) legend(distance{id,2})
end legend(distance((4*j)+1:4*(j+1),2))end
distance = plot_segment_distance(data_using_max, 0);
et j’obtiens
———————fin edit———————————-
