20 mai 2008 – montage d’un systeme de fichier distant (berthemorisot sur douanierousseau)

11h13
J’ai 2 machines, douanierousseau fonctionnant avec ubuntu feisty 7.04 et berthemorisot avec Centos 4
berthemorisot est le serveur de l’unité.
on y stock bcp de choses.
c’est aussi la que se trouvent les script créer pour le labo et fonctionnant avec une interface web
Le but : avoir un répertoire sur ma machine qui reprend le système de fichier de berthemorisot.
Comment faire : en utilisant nfs (Network File System)
il faut configuere les deux coté, le serveur qui va fournir les fichier et le client qui va les lire.
donc ici le serveur est berthemorisot et le client douanierousseau.

1- configuration de berthemorisot
je m’aide de cette page
http://www.centos.org/docs/5/html/Deployment_Guide-en-US/s1-nfs-server-export.html
je configure le fichier /etc/exports pour que douanierousseau soit accepté comme client et que ce soit le repertoire /home/gim qui soit partagé
j’ecris dans ce fichier :
/home/gim/ douanierousseau.gim.pasteur.fr(rw,all_squash)

puis je lance le serveur nfs
sudo /sbin/service nfs start

pour le serveur, c’est fini, passons au client

2- douanierousseau
je m’aide de cette page
http://doc.ubuntu-fr.org/nfs
j’installe le paquet nfs-common
le point de montage sera /home/gim2/berthe/

et je rajoute la ligne
berthemorisot.gim.pasteur.fr:/home/gim /home/gim2/berthe nfs user 0 0
dans /etc/fstab

si le point de montage est tjs vide alors je fais
sudo mount -t nfs berthemorisot.gim.pasteur.fr:/home/gim /home/gim2/berthe

ca marche bien mais je n’ai pas les droit de modification sur les fichiers, c’est chiant

je modifie alors plusieurs choses
d’une par je fais en sorte que l’utilisateur distant soit identifié comme l’utilisateur gim en modifiant /etc/exports
je tape la la ligne
/home/gim/ douanierousseau.gim.pasteur.fr(rw,all_squash,anonuid=501,anongid=100,)

501 étant le code de l’utilisateur gim et 100 celui du group users

et aussi sur douanier rousseau je remplace l’id de l’utilisateur gim2 qui était 1000 par 501 en modifiant le fichier /etc/passwd
j’ai donc maintenant accès a berthemorisot sans probleme.
j’ai modifié /usr/sbin/ntsysv, pour que nfs démarre au démarrage du serveur

17 avril 2008 – génération des tables pour le fichier gal

17h34
les table reprennant les infos de la puce sont
array_chip_type
array_spot
array_spot_type

array_chip_type pourrait a mon avis contenir plus de choses. en gros toutes les infos globales sur une puces.
array_spot rassemble les infos liées aux spot de la puce
et array_spot_type regroupe les différents types de spots, pour l’instant je pense feature, control ou empty.

le type de la puce est SLRI-Yeast-Barcode-13k
les infos sont dans barcode12k_v2final.gal
comment doit se dérouler le programme :
1- récupérer le nom de la puce (et plus tard toute sa description)
2- ajouter ce nom dans le fichier array_chip_type
3- récupérer le fichier gal
4- le parser et pour chaque ligne
5- extraire le type (feature, spotting buffer ou empty)
6- ajouter le type dans la table s’il n’existe pas
7- extraite l’ID
8- eliminer les element en trop de l’ID
9- rechercher dans la base des feature si la feature existe et récupérer sont identifiant
10-ajouter ensuite la ligne dans array_spot en la liant au type de la puce, au type de spot qui va bien et à la feature qui va bien.

en reprennant le code pour introduire les données de features (ici) je me rend compte que le code est assez moyen, et surtout pas adapté a l’utilisation dans un cgi.
de plus il n’utilise pas le nouveau parser que j’ai créé.
donc je reprends tout, enfin en partie.

imaginons la structure de l’interface.
le plus simple, une page web avec un champ parcourir pour aller chercher le gff et un bouton valider
ca envoie la requete sur une page qui va traiter la requete.
cette page va récuperer le fichier et l’enregistrer dans un répertoire temporaire.
puis un objet loader va le tester pour savoir s’il existe bien, et ouvrir le filehandle
puis il va etre transmis au parser qui va vérifier le format (rajouter la fonction check dans le parser) et puis traiter toutes les données.
il faut ensuite une fonction pour ajouter ces données dans la base
il me faut alors une classe pour effectuer cette action.
ca va etre finalement une classe d’interface avec luce_la_puce.
cette classe pourrait s’appeler Garden_state qui traitera par une méthode appropriée l’objet de feature.
il me faut donc un loader.
c’est assez simple, le loader doit récupérer un nom de fichier
vérifier qu’il existye et l’ouvrir
récupérer le parser qui va le traiter
appliquer ce parser et renvoyer l’objet créér

class file_loader():

    def __init__(self, **arg):        filename = arg.get('filename', '')        try :            self.fhin = self._open_file(filename)        except IOError, e :            raise

    def _open_file(self, filename):        '''        fhin = x._open_file(filename)        filename -> str        fhin -> file        open a file and return the filehandle        '''        #test id the file exists and is a file        if not isfile(filename) :            raise IOError, "file %s doesn't exist" % filename        #test if reading access is granted        if not access(filename, R_OK) :            raise IOError, "permission denied for file %s" % filename        return open(filename, 'r')

    def get_filehandle(self):        try :            return self.fhin        except AttributeError :            return None

    def apply_parser(self, parser=None):        try :            if not parser:                parser = self.parser        except AttributeError :            raise

        return parser(file=self.fhin)

    def set_parser(self, parser=None):        self.parser = parser

peut-etre aussi faut’il rajouter le moyen de réouvrir le file handle et puis par la suite une méthode pour copier le fichier dans un autre emplacement.
ensuite il faut donc que je créé l’interface pour luce_la_puce Garden_state

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.

4 avril 2008 – Classification des segments, améliorée

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 donc
distance_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 donc
distance_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 par
hist_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


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

la commande est distance = plot_segment_distance(data_using_max, 0);
et j’obtiens


———————fin edit———————————-

3 avril 2008 – Détail de la classification des segments

11h04 les segments après avoir été répertorié sont classifié selon leur prosition relative au 2 orf les plus proches.
on recense 4 grande classe qui dérive des position relative des orf entre elle.
les orf peuvent etre soit :

• en tandem

• divergentes

• convergentes

les segments peuvent alors se positionner en fonction de ces orf.
ils vont etre soit entre des orf
l’orf an jaune est l’orf upstream du segment
l’orf en rouge est l’orf downstrean du segment

• convergentes

• divergentes

• en tandem sur le meme brin (en sens)

• en tandem sur le brin opposé (en antisens)

On a donc alors 4 classes principales notées C (convergente), D (divergente), Ts (tandem sens) et Ta (tandem antisens).
Cependant, dans chacun des ces cas il y a des sous-catégories selon la position du segment. s’il est bien entre les orf, ou légèrement chevauchant ou totalement inclut dans une orf.
on a pu ainsi definir 16 catégories, 4 pour chaque classe :

• de 1 à 4 pour la classe Ts

• catégorie 1

• catégorie 2

• catégorie 3

• catégorie 4

• de 5 à 8 pour la classe Ta

• catégorie 5

• catégorie 6

• catégorie 7

• catégorie 8

• de 9 à 12 pour la classe D

• catégorie 9

• catégorie 10

• catégorie 11

• catégorie 12

• de 13 à 16 pour la classe C

• catégorie 13

• catégorie 14

• catégorie 15

• catégorie 16

la classification des segments se fait a l’aide d’un programme simple réalisé en VBA sur excel.
il se trouve sur bambino dans F:/users/christophe/search_features_A_and_B.xls a

Sub classification()

Range("AH2").Selectligne = 0Valeur = Range("A2").Offset(ligne, 0).Value

Do Until Valeur = ""

 S = ActiveCell.Offset(ligne, -29) E = ActiveCell.Offset(ligne, -28) m = Int((S + E) / 2) sm = m - S em = E - m If Range("C2").Offset(ligne, 0) = "c" Then     signe = -1 Else     signe = 1 End If up_val = test_prox(ActiveCell.Offset(ligne, -21).Value, ActiveCell.Offset(ligne, -14).Value)

 If up_val Then     col_orf_A = -24     signe2 = 1 Else     col_orf_A = -16     signe2 = -1 End If ORF_A_name = ActiveCell.Offset(ligne, col_orf_A).Value ORF_A_gene = ActiveCell.Offset(ligne, col_orf_A + 1).Value ORF_A_start = m - (ActiveCell.Offset(ligne, col_orf_A + 2).Value * signe * signe2) ORF_A_end = m - (ActiveCell.Offset(ligne, col_orf_A + 3).Value * signe * signe2) If Not IsEmpty(ActiveCell.Offset(ligne, -6)) Then     down_val = True     col_orf_B = -8     m2S = ActiveCell.Offset(ligne, -6).Value     m2E = ActiveCell.Offset(ligne, -5).Value     signe2 = 1 ElseIf Not IsEmpty(ActiveCell.Offset(ligne, -2)) Then     down_val = False     col_orf_B = -4     m2E = -ActiveCell.Offset(ligne, -2).Value     m2S = -ActiveCell.Offset(ligne, -1).Value     signe2 = -1 Else     down_val = test_prox(ActiveCell.Offset(ligne, -18).Value, ActiveCell.Offset(ligne, -9).Value)     If down_val Then         col_orf_B = -20         m2S = ActiveCell.Offset(ligne, -18).Value         m2E = ActiveCell.Offset(ligne, -17).Value         signe2 = 1     Else         col_orf_B = -12         m2E = -ActiveCell.Offset(ligne, -10).Value         m2S = -ActiveCell.Offset(ligne, -9).Value         signe2 = -1     End If End If ORF_B_name = ActiveCell.Offset(ligne, col_orf_B).Value ORF_B_gene = ActiveCell.Offset(ligne, col_orf_B + 1).Value ORF_B_start = m - (ActiveCell.Offset(ligne, col_orf_B + 2).Value * signe * signe2) ORF_B_end = m - (ActiveCell.Offset(ligne, col_orf_B + 3).Value * signe * signe2)

 If (up_val And down_val) Then     ActiveCell.Offset(ligne, 0).Value = "Ts"     base_classe = 0 End If If (Not up_val And Not down_val) Then     ActiveCell.Offset(ligne, 0).Value = "Ta"     base_classe = 4 End If If (Not up_val And down_val) Then     ActiveCell.Offset(ligne, 0).Value = "D"     base_classe = 8 End If If (up_val And Not down_val) Then     ActiveCell.Offset(ligne, 0).Value = "C"     base_classe = 12 End If classe = test_case(sm, em, m2S, m2E) + base_classe ActiveCell.Offset(ligne, 1).Value = ORF_A_name ActiveCell.Offset(ligne, 2).Value = ORF_A_gene ActiveCell.Offset(ligne, 3).Value = ORF_A_start ActiveCell.Offset(ligne, 4).Value = ORF_A_end

 ActiveCell.Offset(ligne, 5).Value = ORF_B_name ActiveCell.Offset(ligne, 6).Value = ORF_B_gene ActiveCell.Offset(ligne, 7).Value = ORF_B_start ActiveCell.Offset(ligne, 8).Value = ORF_B_end ActiveCell.Offset(ligne, 9).Value = classe ligne = ligne + 1 Application.StatusBar = ligne Valeur = Range("A1").Offset(ligne, 0).ValueLoop

End Sub

Function test_prox(val_sens, val_anti)If Abs(val_sens) < Abs(val_anti) Then    test_prox = TrueEnd IfEnd Function

Function test_case(sm, em, m2S, m2E)    test_sS = (m2S - sm < 0)    test_eS = (m2S + em < 0)    test_sE = (m2E - sm < 0)    test_eE = (m2E + em < 0)    If test_sS And test_eS And test_sE And test_eE Then        test_case = 1    ElseIf test_sS And Not test_eS And test_sE And test_eE Then        test_case = 2    ElseIf Not test_sS And Not test_eS And test_sE And test_eE Then        test_case = 3    ElseIf Not test_sS And Not test_eS And test_sE And Not test_eE Then        test_case = 4    Else        test_case = 100    End If

End Function

Je reprends ce code pour voir pourquoi j’ai tant de différence antre les 2 méthodes et aussi pour le commenter un peu.
en le reprennant je constate que j’ai fait une faute et que je n’ai pas pris les bonnes valeur, j’ai oublié de tenir compte de la présence d’une colonne en plus (maxi)
j’ai aussi oublié que cette classification se base sur la position du milieu du segment.
je dois modifier ca aussi.

31 3ars 2008 – parser global, la suite

12h29 : le parser actuellment ne gère pas la possibilité qu’il n’y ai pas de ligne de titre pour les colonnes, je vais le rajouter.
l’option est assez simple. si rien n’est mentionné le parser suppose qu’il y a une ligne de titre a récupéré. si une liste nom de colonne est passée en argument (genre column_names=['colonne A', 'colonne B', 'colonne C'] elle est utilisée comme liste par défaut. si l’argument est column_names=None, alors la liste de nom de colonne sera choisie arbitrairement (par ex col_0, col_1, col_2).
j’en ai profité pour modifié un peu l’__init__.
on peut lui passer une liste d’argument en entrée.
les mots clés reconnus pour l’instant sont filename, file, file_arg, tous les 3 pour le nom du fichier a parser (c’est juste pour laisser plus de liberté) et column_names.

def __init__(self, **arg):        '''        x = Items(**arg)        open and read the file to parse.        parse it according to the specified options        arg keys :        file, filename or file_arg -> str, identifier for the file to parse        column_names -> optional. must be NOT specified if the column names are in the file to parse         must be specified if there is no column names in the file.            - column names are given as a list or tuple            - if value None affected, column names will be assigned automatically.        '''        #initialize a flag dictionary        self.flag = {}        #get the column_names argument        self.flag['column_names'] = arg.get('column_names', False)        #get the filename argument (3 possibilities)        file_arg = arg.get('filename', '')        if not file_arg : file_arg = arg.get('file', '')        if not file_arg : file_arg = arg.get('file_arg', '')        try :            #call get_fhin to open file_arg            fhin = self.get_fhin(file_arg)

        except TypeError, e :            #if error of type of the argument file_arg returned, raise the error            print e            raise        except IOError :            #if problem for opening file_arg, raise the error            raise        #if the file is open and readable start to parse the file        self.parse_file(fhin)        #do some after init actions like closing files        self.__end_init__(fhin)

le module de test tel qu’il est fait vérifie bien que toutes les infos ont été parsées mais ne vérifie pas que ce sont bien les bonnes données qui sont dans les bons attributs.
pour cela il me faut des accesseurs pour les items.
J’ai ajouté la méthode get_attribute() a item. cette méthode renvoie la valeur affectée à un attribut passé en argument et none si l’attribut n’existe pas.

    def get_attribute(self, attr):        '''        val = x.get_attribute(attr)        attr -> str        val -> str, int, float or boolean        return the value of an attribute if exists, None if not        '''

        return self.get(attr, None)

j’ai également ajouté get_attr_list() à items qui renvoie la liste des valeurs affectées a un attribut passé en argument pour l’ensemble des item stockés dans la collection

    def get_attr_list(self, attr):        '''        values = x.get_attr_list(attr)        attr -> str        values -> list of str, int, float, boolean or None        return a list containing values of the given attr for all items in x        '''        return [item.get_attribute(attr) for item in self]

puis j’ai ajouté un nouveau module de test pour verifier que toutes les info sont bien récupérées.

    def test_parsed_data(self):        '''        teste si les donnees parsees sont correctes        '''        header_info_1 = "header_01"        header_info_2 = "header_02"        col_1 = "col_1"        col_2 = "col_2"        col_3 = "col_3"        data_l1_c1 = "data_l1_c1"        data_l1_c2 = "data_l1_c2"        data_l2_c1 = "data_l2_c1"        data_l2_c2 = "data_l2_c2"        test_file='''#%s#%s%s\t%s%s\t%s%s\t%s''' % (header_info_1, header_info_2, col_1, col_2, data_l1_c1, data_l1_c2, data_l2_c1, data_l2_c2)        self.write_tmp_file(filetext = test_file)        test_item = Items(file_arg = 'tmp_file')        self.assert_(test_item.get_header() == [header_info_1,header_info_2])        self.assert_(test_item[0].get_attribute(col_1) == data_l1_c1)        self.assert_(test_item[0].get_attribute(col_2) == data_l1_c2)        self.assert_(test_item[1].get_attribute(col_1) == data_l2_c1)        self.assert_(test_item[1].get_attribute(col_2) == data_l2_c2)        self.assert_(test_item[0].get_attribute(col_3) == None)        self.assert_(test_item[1].get_attribute(col_3) == None)               self.assert_(test_item.get_attr_list(col_1) == [data_l1_c1, data_l2_c1])        self.assert_(test_item.get_attr_list(col_2) == [data_l1_c2, data_l2_c2])        self.assert_(test_item.get_attr_list(col_3) == [None, None])

        remove('tmp_file')

a faire demain -> les données sont pour l’instant parsées indépendamment de leur contenu. elles sont toutes stockées comme texte. il faut que je rajoute le module de conversion…

26 mars 2008 – parser global, browsing du fichier

10h44 : je vais donc écrire les méthodes pour parcourir les fichiers et effectuer les actions adaptées selon la ligne.
il faut prévoir tous les cas. voila ce que je pense faire.
une méthode de pré-traitement, vide par défault pour analyser le header de façon spécifique
une méthode de browsing avec un filtre qui met de coté les lignes du header (ca peut etre optionnel aussi)
les ligne de data sont récupéré dans une liste
une méthode d’analyse de la première ligne pour déterminer les nom de colonne, optionnel aussi
une méthode de traitement de toutes les autres lignes.

j’ai créé ces méthodes.
le programme suit donc la logique suivante :
parsing du header par une méthode spécifique. -> celle-ci est vide par défaut mais elle devra etre surclasser si besoin dans les parser qui dériveront ce celui-ci
browsing de tout le fichier et stockage des lignes dans une liste. au moment de ce stockage, les lignes sont testée pour savoir s’il elle sont valides.
je teste si ce n’est pas une ligne vide et si ce n’est pas du commentaire. a noter que ce que j’appelle commentaire est le header du fichier. je stocke ce header dans une liste mais un objet header serait plus malin. je ferai ca
une fois que j’ai récupéré les lignes valides, je cherche dans la première ligne les nom des colonne. il faut que je puisse m’affranchir de ca. je vais ajouter une option a l’instanciation de item pour préciser s’il y a une ligne de titre de colonne.
et sinon, il faut pouvoir fournir des noms de colonne ou les générer s’ils ne sont pas fournis
ensuite je parcours la liste des données et je parse chaque ligne en splitant sur le séparateur de champ et je créé une instance de Item en lui passant la liste des valeurs et des attributs.
Item va alors faire des couple de valeurs nom de colonne/valeur et va rentrer l’ensemble de ces couple dans son dictionnaire interne (car item hérite de dict).
je récupère chaque item dans une liste
et cette liste est ajoutée a la liste interne de Items (qui hérite de list)
et voial, c’est parsé. c’est beau non ?
je suis pas certain que ce soit l’algo le plus efficace car je parcours plusieurs fois les données. le problème est surtout que je stocke le fichier dans une liste temporaire. ce n’est pas très bien si on ne connait pas la taille du fichier avant.
mais bon, deja je peux faire une vérif de la taille du fichier avant et aussi la capacité des micro permet facilement de gerer ca je pense. les liste ne devraient pas faire plus de qq centaine de millier de ligne, a mon avis moins d’un million. en tou cas pour le moment. ce n’est pas très dur de modifier ce code pour mieux traiter les gros fichiers. il suffit de parser les ligne lors du browsing.

j’ai aussi mis en place la partie test de ces méthodes.
je créer des fichiers en faisant varier le nb de ligne, de colonnes et aussi le nb de ligne de header et je test si le fichier est bien parsé, si j’ai le bon nombre de lignes (taille de la liste items), le bon nombre de colonne (taille des objet item contenus dans items) et le bon nombre de ligne dans headers.
je teste aussi les fichiers vide, les fichier sans header, les fichiers sans valeurs, juste des noms de colonnes. il faut aussi que je teste les fichiers qui n’ont pas de ligne de titre de colonne. mais pour ca je dois modifier le programme pour gérer cette option.
il me faut aussi prévoir des fichiers qui auraient été mal fait. comme par exemple des lignes vides, des colonnes vides, des ligne sans le meme nombre de valeurs

en testant ces fichiers, j’ai vu que
-si une ligne est vide mais qu’elle contient juste une tabulation ca génère une erreur. a voir si c’est important. ca venait peut-etre aussi du fait que c’était la première ligne parsée, avant celle des attributs
-si j’ai moins d’élément dans la ligne que le nb de colonne, alors certain attributs de colonnes vont manquer dans l’objet item. ca peut-etre génant, mais c pas sur non plus
sinon ca marche bien.
lundi il faudra ajouter la gestion de la ligne de titre de colonne absente et gérer ces petites erreurs eventuellement.

13 mars 2008 – identification des tags du fichier gal

18h28 : après avoir vu que le fichier de délétion ne contient pas tous les tags du fichier gal (ici) je recherche dans uin autre fichier si je peux trouver ces séquences.
ce fichier ( /home/gim2/bioinfo/microarray/original_vs_repair.txt) contient les sequence des tag qui ont été identifiés sur la puce comme étant mal séquencé. il donne la coorespondance avec la bonne séquence.

ORF Batch Tag Original_sequence Repaired_sequence Avgerage_signal_repaired Avgerage_signal_original Defect_score Sequencing_data Recommended_version
YAL016W chr00_12 down GTAGACGGAGGATTATTCAC GTAGACGGAGGATTATCAC 1879.5 590.45 4 data_not_checked repaired
YBL022C chr2_1 down CTGCGAGCAATCAGCCGATA CTGTGAGCAATCAGCCGATA 79.6 74.65 4 data_ok repaired
YBL025W chr2_1 down CTGCTGCGAAGTTCCGAGAA CTGCTGCGAAGTTCCCGAGAA 3607.55 2536.5 8 data_not_checked repaired
YBL046W chr2_1 down CTTGCGAAGTGTATTCACCA CTTGGCGAAGTGTATTCACCA 1230.55 648.4 4 data_not_checked repaired
YBL054W chr2_1 down GAAGTGCGGCTAATATGCTA GAAGTGCGGCTAACATGCTA 4143.4 2695.45 7 data_not_checked repaired
YBL060W chr2_1 down GACCCAATTCTACAGCGTAA GACCCATTCTACAGCGTAA 3901.4 2245.7 7 data_not_checked repaired
YBL077W chr2_1 down GAGAGACCATGCAGCCGATA GAGAGACCATGCAGCGATA 3599.85 2338.7 5 data_not_checked repaired
YBL080C chr2_1 down GAGAGTGGAATCGCTCATAA GAGAGTGGATCGCTCATAA 2906.7 1226 9 data_not_checked repaired
YBL085W chr2_1 down GAGGGTCTAATCCTGAGTAA GAGGGTCAATCCTGAGTAA 2805.55 783.9 8 data_not_checked repaired


je vais donc vérifié que les nouvelles séquences ainsi définie correspondent aux séquences absente du fichier de délétion.
pour ca je vais deja parser le fichier.
j’ai fait le parser reparatio_parser.py
je le fait tourner et je fais matcher les Repaired_sequence avec celle du fichier gal

ouverture de barcode12k_v2final.gal en 2.398 sec
barcode12k_v2final.gal contient 27648 lignes
ouverture de original_vs_repair.txt en 0.092 sec
original_vs_repair.txt contient 818 lignes
('Repaired_sequence',)
original_vs_repair.txt contient 818 tags distincts
('Sequence',)
barcode12k_v2final.gal contient 12683 tags distincts
818 elements de original_vs_repair.txt sont absents de barcode12k_v2final.gal
12683 elements de barcode12k_v2final.gal sont absents de original_vs_repair.txt

aucune des sequences n’est dans le fichier gal. j’ai le meme résultat avec la colonne original_sequence
si je recherche ces séquence dans le fichier de délétion j’obtiens

ouverture de Deletion_primers_PCR_sizes.txt en 1.325 sec
Deletion_primers_PCR_sizes.txt contient 6363 lignes
ouverture de original_vs_repair.txt en 0.063 sec
original_vs_repair.txt contient 818 lignes
('Original_sequence',)
original_vs_repair.txt contient 818 tags distincts
('UPTAG', 'DNTAG')
Deletion_primers_PCR_sizes.txt contient 12481 tags distincts
818 elements de original_vs_repair.txt sont absents de Deletion_primers_PCR_sizes.txt
12481 elements de Deletion_primers_PCR_sizes.txt sont absents de original_vs_repair.txt

j’ai le meme résultat avec la colonne repaired_sequence
il semble donc qu’il n’y ai aucune correspondance entre ce fichier et les séquences existante

en regardant de plus pret on constate que les séquence du fichier de reparation semble a l’envers
je modifie mon parser pour inverser les séquences. j’obtiens en comparant les séquences d’origine avec le fichier de délétion.

ouverture de Deletion_primers_PCR_sizes.txt en 1.504 sec
Deletion_primers_PCR_sizes.txt contient 6363 lignes
ouverture de original_vs_repair.txt en 0.080 sec
original_vs_repair.txt contient 818 lignes
('Original_sequence',)
original_vs_repair.txt contient 818 tags distincts
('UPTAG', 'DNTAG')
Deletion_primers_PCR_sizes.txt contient 12481 tags distincts
1 elements de original_vs_repair.txt sont absents de Deletion_primers_PCR_sizes.txt
sequences absentes set(['GCCTCTTCGCAAAGTAACAA'])
11664 elements de Deletion_primers_PCR_sizes.txt sont absents de original_vs_repair.txt

je retrouve bien toute les séquence originelles
cependant si je compare les séquences réparée avec le fichier gal j’obtiens

ouverture de barcode12k_v2final.gal en 2.186 sec
barcode12k_v2final.gal contient 27648 lignes
ouverture de Deletion_primers_PCR_sizes.txt en 1.369 sec
Deletion_primers_PCR_sizes.txt contient 6363 lignes
ouverture de original_vs_repair.txt en 0.083 sec
original_vs_repair.txt contient 818 lignes
('Repaired_sequence',)
original_vs_repair.txt contient 818 tags distincts
('Sequence',)
barcode12k_v2final.gal contient 12683 tags distincts
595 elements de original_vs_repair.txt sont absents de barcode12k_v2final.gal
12460 elements de barcode12k_v2final.gal sont absents de original_vs_repair.txt

je retrouve donc environ 200 séquences mais presque 600 sont encore absentes.

13 mars 2008 – association des tag des spot des microarray avec les ORF

10h45 : hier j’ai pu constater que environ une centaine d’anotations des spots des puces n’était pas présente dans le fichier gff. voir ici
mais en prenant au hasard quelque orf j’ai juste constaté qu’elle était bien présente dans le fichier GFF mais sous le nom orf-A ou orf-B. du coup la correspondance ne marche plus en python.
je vais essayer de voir si je peux faire matcher ces orf quand meme…
pour cela j’ai une première idée, c’est tout simplement de ne garder que le nom d’orf sans le -A quand j’extrais les noms des orf avec le parser gff_parser et de voir si j’ai des meilleurs match.
par cette méthode ce sont juste 29 élément que je ne retrouve plus dans le fichier gff

YKL199C YAR044W YCL060C YCL006C YCL062W YBR100W YER108C YLR391W YKL158W YOR240W
YAR037W YKL200C YCL012W YCL053C YAR043C YGL046W YFL018W YFL006W YJL018W YOL053C
YML010C YDR474C YJL021C YML033W YOR088W YCL013W YFL035C YAR040C YCL003W

Cependant je n’aime pas trop cette méthode car je perd de l’info.
finalement je parse le fichier gal (identique a priori au barcode et je récupère les id des spot. avec une expression régulière je retire les -D, -U, -R présents a la fin des id
re.sub(r'-[RUD]+[1-9]*','',tag) j’ai mis des chiffres a la suite des D, U et R car parfois il y en a dans les id.
du coup j’obtiens 40 orf absente du gff

YFL035C-B YFL035C-A YKL199C YBR090C-A YHR079C-B YAR044W YCL060C YCL006C YCL062W YAR037W
YER108C YLR391W YKL158W YHR039C-B YOR240W Spotting Buffer Empty YKL200C YCL012W YFR024C
YCL053C YAR043C YGL046W YFL006W YJL018W YBR100W YDR474C YFL018W-A YCL026C YJL021C
YJL206C-A YML033W YOR088W YOL053C-A YCL013W YML010C-B YDL134C-A YFL035C YAR040C YCL003W

si je fais une recherche sur le premier, YFL035C-B, je trouve sur le chr VI l’orf YFL034-B qui a pour alias YFL035-B. je vais donc modifié le parser gff pour aussi tenir compte des alias
de cette manière je n’ai plus que 7 orf non présentes dans le gff

YCL053C YAR043C Spotting Buffer YCL013W Empty YCL026C YCL006C YAR040C YAR037W


10 mars 2008 – interface d’ajout des puces (fichier gal)

15h53 : après avoir ajouter les features, il faut créer une nouvelle puce, a partir du fichier gal correspondant (fichier de définition des puce d’afficmetrix) et peut-etre qq infos en plus.
deja j’ai plusieurs questions.
il est facile de créer une puce et de mettre ne relation les element du fichiers gal avec ceux de la table des feature.
cependant que faut il faire quand on modifie le fichier des features ? faut t’il créer une nouvelle puce liée a ces nouvelle feature qui seront utilisées a partir de ce moment, ou alors remplacer dans la liste les lien vers les nouvelles features.
idem si on importe un nouveau fichier gal, faut-il ajouter en plus ou mettre a jour. dans ce cas je pense qu’il faut ajouter car les précédentes manip auront été faite avec une ancienne puce et donc un ancien fichier gal.
il semble que finalement le plus interessant c’est de toujours tout mettre a jour.
par contre je pense qu’il faut coupler tout ca avec la base de délétion disponible ici
je vais surement créer une nouvelle table pour ca.