2008 avril « Xtof78000’s Weblog

Archives pour avril 2008

29 avril 2008 – interface sganarrayl

avril 29, 2008

11h24 :
Pour faire l’interface d’ajout du fichier gal, j’ai voulu reprendre celle d’ajout du fichier gff.
celle-ci étant assez moyenne, j’ai donc décidé de la reprendre.
pour que ce soit cohérent avec l’interface finale de sganarrayl, il faut deja pensé a l’intégralité de l’interface.
Depuis la page d’accueil pour ajouter des puce, on pourra donc voir la liste des puces disponibles et aussi du fichier de définition des features utilisé.
on pourra depuis cette page aussi ajouter de nouveaux fichiers de puce et de feature.
l’ajout de feature va renvoyer sur la page add_gff.php
cette page propose un champ browse pour parcourir le disque local et ajouter un fichier.
l’envoie du formulaire va renvoyer sur la meme page.
cette fois ci, comme il y a un fichier uploader, ca va déclencher le script php.
celui-ci copie le fichier uploadé dans un répertoire spécifique.
il va aussi modifier la page pour faire apparaitre le nom du fichier traité et faire disparaitre le champs browse.
la page va aussi avoir un script javascript (add_gff.js).
celui-ci va se déclencher s’il y a un nom de fichier a traiter.
dans ce cas il va faire appelle au script add_gff.py pour traiter le fichier gff et renvoyer la réponse qui correspond au traitement et l’afficher.
add_gff.py parse le fichier gff à l’aide des du parser gff de parser_tools puis utilise la classe garden_state du module sganarrayl_interface pour se connecter à la base luce_la_puce et ajouter les feature par la méthode add_feature_data.

add_gff.py renvoie un fichier xml avec les balises ou selon le déroulement du code.
le script javascript va interpréter ce xml et afficher les messages correspondants dans une balise prévue a cette effet et ingénieusement intitulée gff
voici le flux de données

L’interface pour le fichier gff fonctionne bien maintenant et le code est propre.
passons au fichier gal

déroulement de la séquence.
l’utilisateur doit ajouter une puce dans la base.
il va alors aller sur la page add_gal.php.
il n’y a que 2 champs pour le moment.
le premier pour donner un nom à la puce.
le 2eme pour mettre le fichier gal qui contiendra les info sur les spots.
le nom de la puce sera celui utilisé par la suite pour faire référence à celle-ci.
le script va ensuite ajouter le nom de la puce dans la table array_chip_types s’il n’existe pas encore. s’il existe deja, l’ajout de la puce est refusé. par la suite il devra permettre la mise a jour de la puce.
ensuite il faut parcourir chaque spot, déterminer son type, vérifier dans la table array_spot_type si le type existe et sinon l’ajouter. on recherche ensuite dans les features de la table bank_orf qui ont la meme annotation que le spot.
si on trouve on garde la ref de la feature, sinon, la ref est nulle. on ajoute ensuite les info du spot dans la table array_spots, en ajoutant les liaisons vers les tables array_chip_types array_spot_type et bank_orf.

todo -> modifier sganarrayl_interface pour pouvoir sélectionner les id des orf ayant l’annotation du spot

Publié dans SganArrayL, ajax, gal, gff, pygresql, python | Laisser un commentaire »

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

avril 17, 2008

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

Publié dans Uncategorized | Laisser un commentaire »

15 avril 2008 – Parser pour le fichier gff

avril 16, 2008

le parser pour fichier GFF est assez simple
description du format gff –> ici

un header de taille variable marqué par “#”
puis des données en colonne séparées par des tabulations.
il y a 9 colonnes.
il n’y a pas de ligne de titre
exemple :

less /home/gim2/bank/saccharomyces_cerevisiae.gff

##gff-version 3#date Wed Jul 18 19:35:09 2007## Saccharomyces cerevisiae S288C genome## Features from the 16 nuclear chromosomes labeled chrI to chrXVI,# plus the mitochondrial genome labeled chrMito and the 2-micron plasmid.## Created by Saccharomyces Genome Database (http://www.yeastgenome.org/)## Weekly updates of this file are available via Anonymous FTP from:# ftp://ftp.yeastgenome.org/yeast/data_download/chromosomal_feature/saccharomyces_cerevisiae.gff## Please send comments and suggestions to yeast-curator@yeastgenome.org## SGD is funded as a National Human Genome Research Institute Biomedical Informatics Resource from# the U. S. National Institutes of Health to Stanford University.  The staff of SGD is listed at:# http://www.yeastgenome.org/SGD-staff.html#chrI    SGD     chromosome      1       230208  .       .       .       ID=chrI;dbxref=NCBI:NC_001133chrI    SGD     repeat_family   1       62      .       -       .       ID=TEL01L-TR;Name=TEL01L-TR;Note=Terminal%20stretch%20of%20telomeric%20repeats%20on%20the%20left%20arm%20of%20Chromosome%20I;dbxref=SGD:S000028864

le contenu des colonnes est (dans l’ordre) :

seqname

Le nom de la séquences, dans le cas de la levure, le chromosome

source

La source de la feature. soit le programme ayant réalisé la prédiction, ou la base de donnée publique… dans le cas de la levure, c’est en général SGD.

feature

Le nom de la feature si possible un nom standardisé.

start, end

debut et fin de la feature, fin doit etre supérieur a debut

score

score de la feature ou ‘.’ s’il n’y en a pas

strand

‘+’ ou ‘-’ selon le brin ou alors ‘.’ si ca ne s’applique pas

frame

‘0′, ‘1′, ‘2′, ou ‘.’ indique la position du premier codon

Attribute

Cette colonne est différente car elle contient nom pas une valeur mais une liste d’attribut. ceux-ci sont séparé par des “;” et

Le fichier peut contenir des meta données identifées par un “#”

le parser doit donc
parser le header par la méthode par défaut
utiliser une liste de nom de colonne prédéfinie
parser les lignes
puis pour chaque item parser l’élément de la dernière colonne pour en extraire les attributs et les ajouter a item.

class Gff_parser(Items):    def __init__(self, **arg):        '''        x = Gff_parser(**arg)        '''

        arg['column_names']=['chr','origin','type','start','end','score','strand','frame', 'group']

        self.attr_group=set([])        #call the init method of Items        Items.__init__(self, **arg)        self.attr_list.extend(self.attr_group)

    def parse_group(self, item):        group = item.get_attribute('group')

        groups = split(group, ';')

        group_items=(self.parse_group_item(item) for item in groups)

        item.set_attribute(group_items)

        self.attr_group.update(item.get_attr_list())

        return item

    def special_data_parsing_fn(self, item):        return self.parse_group(item)

    def parse_group_item(self, group_item):        return split(group_item, '=')

j’ai modifié la méthode d’assignation des item.
j’ai rajouté une fonction qui s’applique sur chaque item juste après sa création.
cette méthode sera a surclasser dans les parser si besoin

    def list2item(self, data_list):        '''        data_item = x.list2item(data_list)        data_list -> list of list        data_item -> list        create instances of Item for all element of data_list        return all instances in a list        '''        try :            #create instance of Item for all element of data_list            return (self.special_data_parsing_fn(Item(x, self.attr_list)) for x in data_list)        except TypeError :            #if data_list can't be map, return an ampty list            return []

    def special_data_parsing_fn(self, item):        return item

j’ai ajouté un module de test gff qui reste a completer.

    def test_gff_parser(self):        test_gff = Gff_parser(file_arg = '/home/gim2/bank//saccharomyces_cerevisiae.gff')        self.assert_(isinstance(test_gff,Gff_parser))        print test_gff[0].keys()

j’ai en sortie

tmp_file format is not correct..col_0 is not a valide attribute....the argument is not a regular file objectfile toto doesn't existpermission denied for file tmp_parse.txt..----------------------------------------------------------------------Ran 8 tests in 4.603s

OK

Publié dans gff, parser | Laisser un commentaire »

15 avril 2008 – Parser pour le fichier gal

avril 15, 2008

11h04
maintenant que le parser global est construit, je vais faire un parser pour le fichier gal dérivé de celui-ci
Format d’un fichier GAL -> GenePix Array List
exemple du fichier /home/gim2/bioinfo/microarray/barcode12k_v2final.gal
la commande
less /home/gim2/bioinfo/microarray/barcode12k_v2final.gal
donne
ATF 1 51 6 Type=GenePix ArrayList V1.0 BlockCount=48 BlockType=0 "Block1= 3530, 8000, 120, 24, 180, 24, 180" "Block2= 8030, 8000, 120, 24, 180, 24, 180" "Block3= 12530, 8000, 120, 24, 180, 24, 180" "Block4= 17030, 8000, 120, 24, 180, 24, 180" "Block5= 3530, 12500, 120, 24, 180, 24, 180" "Block6= 8030, 12500, 120, 24, 180, 24, 180" "Block7= 12530, 12500, 120, 24, 180, 24, 180" "Block8= 17030, 12500, 120, 24, 180, 24, 180" "Block9= 3530, 17000, 120, 24, 180, 24, 180" "Block10= 8030, 17000, 120, 24, 180, 24, 180" "Block11= 12530, 17000, 120, 24, 180, 24, 180" "Block12= 17030, 17000, 120, 24, 180, 24, 180" "Block13= 3530, 21500, 120, 24, 180, 24, 180" "Block14= 8030, 21500, 120, 24, 180, 24, 180" "Block15= 12530, 21500, 120, 24, 180, 24, 180" "Block16= 17030, 21500, 120, 24, 180, 24, 180" "Block17= 3530, 26000, 120, 24, 180, 24, 180" "Block18= 8030, 26000, 120, 24, 180, 24, 180" "Block19= 12530, 26000, 120, 24, 180, 24, 180" "Block20= 17030, 26000, 120, 24, 180, 24, 180" "Block21= 3530, 30500, 120, 24, 180, 24, 180" "Block22= 8030, 30500, 120, 24, 180, 24, 180" "Block23= 12530, 30500, 120, 24, 180, 24, 180" "Block24= 17030, 30500, 120, 24, 180, 24, 180" "Block25= 3530, 35000, 120, 24, 180, 24, 180" "Block26= 8030, 35000, 120, 24, 180, 24, 180" "Block27= 12530, 35000, 120, 24, 180, 24, 180" "Block28= 17030, 35000, 120, 24, 180, 24, 180" "Block29= 3530, 39500, 120, 24, 180, 24, 180" "Block30= 8030, 39500, 120, 24, 180, 24, 180" "Block31= 12530, 39500, 120, 24, 180, 24, 180" "Block32= 17030, 39500, 120, 24, 180, 24, 180" "Block33= 3530, 44000, 120, 24, 180, 24, 180" "Block34= 8030, 44000, 120, 24, 180, 24, 180" "Block35= 12530, 44000, 120, 24, 180, 24, 180" "Block36= 17030, 44000, 120, 24, 180, 24, 180" "Block37= 3530, 48500, 120, 24, 180, 24, 180" "Block38= 8030, 48500, 120, 24, 180, 24, 180" "Block39= 12530, 48500, 120, 24, 180, 24, 180" "Block40= 17030, 48500, 120, 24, 180, 24, 180" "Block41= 3530, 53000, 120, 24, 180, 24, 180" "Block42= 8030, 53000, 120, 24, 180, 24, 180" "Block43= 12530, 53000, 120, 24, 180, 24, 180" "Block44= 17030, 53000, 120, 24, 180, 24, 180" "Block45= 3530, 57500, 120, 24, 180, 24, 180" "Block46= 8030, 57500, 120, 24, 180, 24, 180" "Block47= 12530, 57500, 120, 24, 180, 24, 180" "Block48= 17030, 57500, 120, 24, 180, 24, 180" Block Column Row ID Name Sequence 1 1 1 YHR185C-D-R PFS1 AGAGATCCACTTCCCATAATT 1 2 1 YHR185C-D-R PFS1 AGAGATCCACTTCCCATAATT 1 3 1 YIR043C-D CTGCACATTCGATTACAGCG 1 4 1 YIR043C-D CTGCACATTCGATTACAGCG 1 5 1 YHR067W-D RMD12 TTATGCACCGTGACGAGGCT 1 6 1 YHR067W-D RMD12 TTATGCACCGTGACGAGGCT 1 7 1 YBR296C-U PHO89 CACGACCGACAATATGGTGA 1 8 1 YBR296C-U PHO89 CACGACCGACAATATGGTGA 1 9 1 YOL076W-U MDM2 ACGGATGGATCAGTTGCTAT 1 10 1 YOL076W-U MDM2 ACGGATGGATCAGTTGCTAT 1 11 1 YNL204C-U SPS18 GCGCTGGCACAAGAATACCA 1 12 1 YNL204C-U SPS18 GCGCTGGCACAAGAATACCA 1 13 1 YLR205C-U HMX1 AACTGAACATACCCGGTGAC 1 14 1 YLR205C-U HMX1 AACTGAACATACCCGGTGAC 1 15 1 YJL039C-D NUP192 TTGAGCCCGATCAGTCGATG 1 16 1 YJL039C-D NUP192 TTGAGCCCGATCAGTCGATG 1 17 1 YGR029W-D ERV1 TTACGTCCCGGATGCCGTTT 1 18 1 YGR029W-D ERV1 TTACGTCCCGGATGCCGTTT 1 19 1 YEL072W-D RMD6 TGCAGCACGCAAGACCATGA 1 20 1 YEL072W-D RMD6 TGCAGCACGCAAGACCATGA 1 21 1 YDR022C-D CIS1 GGTCGATAATAACACGCCAC 1 22 1 YDR022C-D CIS1 GGTCGATAATAACACGCCAC 1 23 1 YBL086C-U ACCACTCACTAAGGAGGATC 1 24 1 YBL086C-U ACCACTCACTAAGGAGGATC 1 1 2 YLR327C-D-R GATAACGACTCAGTGC 1 2 2 YLR327C-D-R GATAACGACTCAGTGC 1 3 2 YIL033C-D BCY1 GGATATTAGCCATCTACGTG

description des différentes parties
le fichier est spéraré en 3 parties
les 2 première ligne donne des info sur la structure du fichier
il y a ensuite le header du fichier
puis les données
1ere ligne -> ATF : conforme au format Axon Text File
1 : version du format
2eme ligne -> 51 : nombre de ligne dans le header
6 : nombre de colonne dans les données
3éme à 53 ème ligne ->header
Type=GenePix ArrayList V1.0 -> type du fichier
BlockCount=48 -> nombre de block
BlockType=0 ->Type de block decrit :
0 = rectangulaire.
1 = orange-packing #1.
2 = orange-packing #2.
“Block1= 3530, 8000, 120, 24, 180, 24, 180″ ->position et dimension de chaque block
xOrigin : position x du centre de la feature la plus en haut a gauche
yOrigin : position y du centre de la feature la plus en haut a gauche
FeatureDiameter : Diamètre des features du block
xFeature : Nombre de colonnes de features dans le block
xSpacing : espacement des colonnes
yFeature : Nombre de lignes de features dans le block
ySpacing : espacement des lignes
54eme ligne -> ligne de titre des colonnes
55eme ligne a la fin -> information pour chaque features

Le parser doit donc spécifiquement
lire la première ligne et stocker le format de fichier et sa valeur
lire la 2eme ligne et stocker le nb de ligne de header et le nb de colonne (pas vraimenet utile)
lire toute les ligne du header, les parser pour séparer la clé des données
(possible de rajouter un mini parser pour les données de block)
puis lire les données.
je vais donc faire un parser héritant de la classe Items.
je vais surclasser parse_header pour traiter les 1ere ligne et le header
puis parser le rester par les méthodes de Items.

le code du parser


class Gal_parser(Items):'''Gal_parserparser for *.gal filesx = Gal_parser(**arg)options :file_arg -> file path or file handleitem_list -> list of item instances

'''

def __init__(self, **arg):'''x = Gal_parser(**arg)'''#call the init method of ItemsItems.__init__(self, **arg)

def parse_header(self, fhin):'''x.parse_header(fhin)fhin -> fileparse the header of the file.

'''#get the format of the file (1st line of the file)file_format = fhin.readline()#check if it's a ATF fileif not "ATF" in file_format :    raise FormatError("%s format is not correct" % fhin.name)#store the file_format dataself.file_format=split(file_format, '\t')

#get the file structure data (2nd line of the file)try :    self.nb_header_rows, self.nb_cols=split(fhin.readline(), '\t')except ValueError :    self.nb_header_rows = fhin.readline()

self.nb_header_rows = int(self.nb_header_rows)

#create the structure for the headerself.header={}#get all the header lines and parse themfor i in range(self.nb_header_rows) :    try :        k,v = self.parse_header_line(fhin.readline())    except AttributeError :        continue    except ValueError :        continue    self.header[k] = v#set the parse_header flag to true, so it won't be parsed again laterself.flag['parse_header'] = True

def parse_header_line(self, line):'''k,v = x.parse_header_line(line)clean the line and split the key and the value'''#cleaning of the linetry :    #remove backspace    line = remove_end_line(line)    #remove "    line=replace(line,'"', '')    #remove tab    line=replace(line,'\t', '')

except AttributeError :    print "no replace for %s" % line    raise#split the line with the field separator =#and return the resulting tuple ok key/valuetry :    return split(line,"=")except ValueError :    print "no = in %s" % line    raise

def get_header(self, *arg):'''header = x.get_header(*arg)return the header dictionnaryarg ->    keys : return keys of header    values : return values of header    items : return tuple ok key/value of header'''#check the arg optionif 'keys' in arg : return self.header.keys()if 'values' in arg : return self.header.values()if 'items' in arg : return self.header.items()#if no optionreturn self.header

class FormatError(exceptions.Exception):def __init__(self, msg=''):self.errmsg = msg

def __str__(self):return self.errmsg

j’ai ajouté un classe FormatError pour gérer les erreur de format de fichier.
j’ai modifié la fonction de test

def test_gal_parser(self):'''test si l'heritage fonctionne bien'''file_format="ATF"file_format_version = 1header_row_nb = 2data_col_nb = 3header_key_1 = "key1"header_val_1 = "val1"header_key_2 = "key2"header_val_2 = "val2"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"data_l2_c3 = "data_c3"data_l1_c3 = "data_c3"test_file='''%s\t%s%s\t%s%s=%s"%s=%s"%s\t%s\t%s%s\t%s\t%s%s\t%s\t%s''' % (file_format, file_format_version,header_row_nb, data_col_nb,header_key_1, header_val_1,header_key_2, header_val_2,col_1, col_2, col_3,data_l1_c1, data_l1_c2, data_l1_c3,data_l2_c1, data_l2_c2, data_l2_c3)self.write_tmp_file(filetext = test_file)test_gal = Gal_parser(file_arg = 'tmp_file')self.assert_(isinstance(test_gal,Gal_parser))self.assert_(len(test_gal.get_header())==2)

self.assert_(test_gal.get_header('keys')==[header_key_2, header_key_1])self.assert_(test_gal.get_header('values')==[header_val_2, header_val_1])self.assert_(len(test_gal)==2)self.assert_(test_gal[0].get_attribute(col_1)==data_l1_c1)

file_format="ATG"test_file='''%s\t%s%s\t%s%s=%s"%s=%s"%s\t%s\t%s%s\t%s\t%s%s\t%s\t%s''' % (file_format, file_format_version,header_row_nb, data_col_nb,header_key_1, header_val_1,header_key_2, header_val_2,col_1, col_2, col_3,data_l1_c1, data_l1_c2, data_l1_c3,data_l2_c1, data_l2_c2, data_l2_c3)self.write_tmp_file(filetext = test_file)self.assertRaises(FormatError, Gal_parser, filename = 'tmp_file')

remove('tmp_file')

test_gal = Gal_parser(file_arg = '/home/gim2/bioinfo/microarray/barcode12k_v2final.gal');self.assert_(isinstance(test_gal,Gal_parser))self.assert_(len(test_gal.get_header())==51)

et le test me renvoie

tmp_file format is not correct.col_0 is not a valide attribute....the argument is not a regular file objectfile toto doesn't existpermission denied for file tmp_parse.txt..----------------------------------------------------------------------Ran 7 tests in 2.844s

OK

j’ai donc un parser de fichier gal qui fonctionne bien.
quel est le but de tout ca ?
l’idée est de parser le fichier de description de la puce et de mettre les infos dans la base de donnée. je cherche également a liéer les spot aux orf du genome de la levure.
j’ai précédemment examiné le fichier gal et d’autre fichiers décrivant les puces (12 mars, 13 mars (1) et 13 mars (2)
l’analyse de ces fichier a permis de constater que l’annotation de certain spot ne coorespond plus a une orf existante. mais que cette ancienne annotation est généralement toujours présente dans les alias des orf.
l’analyse plus fine de cela a meme permis de déterminer que seule 7 annotations avaient réellement totalement disparues.
YAR037W YAR040C YAR043C YCL006C YCL013W YCL026C YCL053C
ce qui correspond a 30 spot au total
grep YCL053C ~/bioinfo/microarray/barcode12k_v2final.gal
1 23 19 YCL053C-U CTATTGTTGAAATGCCGGGA 1 24 19 YCL053C-U CTATTGTTGAAATGCCGGGA 39 23 20 YCL053C-D CATAGTCGAGAACCGGAGAC 39 24 20 YCL053C-D CATAGTCGAGAACCGGAGAC
grep YAR043C ~/bioinfo/microarray/barcode12k_v2final.gal
34 3 8 YAR043C-U ATTCTAGCGGCAGATCCGTG 34 4 8 YAR043C-U ATTCTAGCGGCAGATCCGTG
grep YCL013W ~/bioinfo/microarray/barcode12k_v2final.gal
13 23 22 YCL013W-D CGCTCGAACATAATTGGGTA 13 24 22 YCL013W-D CGCTCGAACATAATTGGGTA 25 23 22 YCL013W-U CCTGTCAGTAAACCGAGAGA 25 24 22 YCL013W-U CCTGTCAGTAAACCGAGAGA
grep YCL026C ~/bioinfo/microarray/barcode12k_v2final.gal
5 23 21 YCL026C-D CCTCCGAACAGAGAGTCTTA 5 24 21 YCL026C-D CCTCCGAACAGAGAGTCTTA 13 3 4 YCL026C-A-U FRM2 GCTCACCGAACATCAGATTA 13 4 4 YCL026C-A-U FRM2 GCTCACCGAACATCAGATTA 17 3 21 YCL026C-U-R CCTCTGCTAAGTAGTAGA 17 4 21 YCL026C-U-R CCTCTGCTAAGTAGTAGA 17 23 21 YCL026C-U CCCTCTGCTAAAGTAGTAGA 17 24 21 YCL026C-U CCCTCTGCTAAAGTAGTAGA 42 5 19 YCL026C-A-D FRM2 GGCGGACTACAACACATTCA 42 6 19 YCL026C-A-D FRM2 GGCGGACTACAACACATTCA
grep YCL006C ~/bioinfo/microarray/barcode12k_v2final.gal
23 23 23 YCL006C-D CCCGCTAGTCAATAATCGTA 23 24 23 YCL006C-D CCCGCTAGTCAATAATCGTA 29 5 15 YCL006C-D CCCGCTAGTCAATAATCGTA 29 6 15 YCL006C-D CCCGCTAGTCAATAATCGTA 35 23 23 YCL006C-U CCTGAAGATAAATCCCGTCA 35 24 23 YCL006C-U CCTGAAGATAAATCCCGTCA 41 5 15 YCL006C-U CCTGAAGATAAATCCCGTCA 41 6 15 YCL006C-U CCTGAAGATAAATCCCGTCA
grep YAR040C ~/bioinfo/microarray/barcode12k_v2final.gal
7 3 9 YAR040C-U CCATCTCAGTGGGTGCAATG 7 4 9 YAR040C-U CCATCTCAGTGGGTGCAATG
grep YAR037W ~/bioinfo/microarray/barcode12k_v2final.gal
23 3 9 YAR037W-U AGCTAGACTATCGCCCAATG 23 4 9 YAR037W-U AGCTAGACTATCGCCCAATG
a
l’ideal serait de faire cette identification lorsque j’ajoute le fichier gal dans la base.
si on suppose que la liste des orf est dans la base de donnée.
il faut donc extraire le nom des annotation des spot, retirer les caractère qui ont été ajouté pour les puc (U, D et R) en faire une liste d’occurences uniques et les rechercher dans la base de données.
je vais faire une premiere version mais en utilisant un fichier gff au lieu de de la base de donnée.
il faut deja récupéré les ID des spots.
pour ca j’utilise la méthode get_attr_list qui me renvoie la liste des valeur d’un attribut pour l’ensemble des items.
j’utilise ensuite un set pour stocker ces valeur, ainsi je n’ai que celles qui sont distincte.
j’utilise ensuite une expression régulière pour eliminer les -D, -U, et -R

def extract_annotation_id(annot_id):return re.sub(r'-[RUD]+[1-9]*','',annot_id)

def analyse_gal_file():gal_file = '/home/gim2/bioinfo/microarray/barcode12k_v2final.gal'gal = Gal_parser(file_arg = gal_file)id = gal.get_attr_list('ID')print "il y a %i ID de spot dans %s" % (len(id), gal_file)id_set = set(id)print "il y a %i ID distintes dans %s" % (len(id_set), gal_file)id=map(extract_annotation_id, id)id_set = set(id)print "il y a %i annotations d'orf dans %s" % (len(id_set), gal_file)

j’obtiens

il y a 27648 ID de spot dans /home/gim2/bioinfo/microarray/barcode12k_v2final.galil y a 12677 ID distintes dans /home/gim2/bioinfo/microarray/barcode12k_v2final.galil y a 5920 annotations d'orf dans /home/gim2/bioinfo/microarray/barcode12k_v2final.gal

j’arrive donc a récupérer les 5920 orf utilisée dans ces puces
il faut ensuite aller les comparer avec celles du fichiers gff
je refais donc un parser gff rapide (ici)

puis je sélectionne les id des orf
gff_file = '/home/gim2/bank//saccharomyces_cerevisiae.gff' gff = Gff_parser(file_arg = gff_file) id2 = gff.get_attr_list('ID')
il y a 16316 ID dans /home/gim2/bank//saccharomyces_cerevisiae.gff
si je prends les tag name et gene en plus j’ai :
id2 = gff.get_attr_list('ID') + gff.get_attr_list('Name') + gff.get_attr_list('gene')
il y a 48948 ID, Name et gene dans /home/gim2/bank//saccharomyces_cerevisiae.gff
si j’utilise les set pour connaitre n’avoir que des valeurs distinctes j’ai :
il y a 16316 ID dans /home/gim2/bank//saccharomyces_cerevisiae.gff il y a 8098 ID distinctes dans /home/gim2/bank//saccharomyces_cerevisiae.gff il y a 48948 ID, Name et gene dans /home/gim2/bank//saccharomyces_cerevisiae.gff il y a 12879 ID, Name et gene distincts dans /home/gim2/bank//saccharomyces_cerevisiae.gff
j’ai besoin aussi des alias des orf car ce sont parfois les anciens noms de ces orfs

alias=gff.get_attr_list('Alias')  for x in alias :      if not x : continue      x=split(x,',')      id2 += x  print len(id2)  print 'il y a %i ID, Name, gene et alias dans %s' % (len(id2), gff_file)  id_set2 = set(id2)  print 'il y a %i ID, Name, gene et alias distincts dans %s' % (len(id_set2), gff_file)

et j’obtiens

il y a 65786 ID, Name, gene et alias dans /home/gim2/bank//saccharomyces_cerevisiae.gffil y a 15774 ID, Name, gene et alias distincts dans /home/gim2/bank//saccharomyces_cerevisiae.gff

je vais donc comparer les id du fichier gal et ceux du fichier gff en faisant une simple différence de set
print "il y a %i ID presente dans %s qui sont absentes de %s" % (len(id_set2), gal_file, gff_file) print "ce sont les ID : ", id_set.difference(id_set2)

il y a 15774 ID presentes dans /home/gim2/bioinfo/microarray/barcode12k_v2final.gal qui sont absentes de /home/gim2/bank/saccharomyces_cerevisiae.gffce sont les ID :  set(['YCL053C', 'YAR043C', 'Spotting Buffer', 'YCL013W', 'Empty', 'YCL026C', 'YCL006C', 'YAR040C', 'YAR037W'])

Ca correspond a celles trouvées précédemment

Publié dans SganArrayL, gal, parser | Laisser un commentaire »

14 avril 2008 – ajout du module de conversion des données

avril 14, 2008

10h48
suite de la construction du parser global.
pour mieux utiliser les données dans le parser, il est utile de les convertir au bon format.
elles sont toutes reconnues comme texte.
les seuls autres types pouvant exister pour des données sont integer ou float.
il pourrait y avoir des booléens mais on va partir du principe qu’il n’y en a pas.
les types listes ou autre propres a python ne sont pas présents dans des fichiers de données.
la conversion des données est intéressante mais ralenti le parsing.
il me semble que malgré tout il y a un gain de temps par la suite
La conversion utilise la fonction try. qui permet e s’affranchir de pas mal de test qui sont réalisé par les fonction de conversion elle-meme. ces fonctions sont int et float.
code du convertisseur

 def set_type(self, value):     '''     val = x.set_type(value)     value -> str     val -> str, int or float     convert the data to the correct data type     '''     #store the data into the return variable     val = value     #try converting     try :         #try first to convert into a float         val = float(value)         #if succeed, try to convert into an integer         val = int(value)     except ValueError:         #catch a converting error, either from float or int         pass     #return the stored value     return val

et le code du test correspondant

 def test_type_parsed_data(self):     '''     teste si le type des donnees parsees est correcte     '''     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 = "texte"     data_l1_c2 = 1     data_l2_c1 = 1.5     data_l2_c2 = "1e+5"     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_(type(test_item[0].get_attribute(col_1)) == type(data_l1_c1))     self.assert_(type(test_item[0].get_attribute(col_2)) == type(data_l1_c2))     self.assert_(type(test_item[1].get_attribute(col_1)) == type(data_l2_c1))     self.assert_(type(test_item[1].get_attribute(col_2)) == type(float(data_l2_c2)))

     remove('tmp_file')

ce qui donne en sortie de la fonction test

..the argument is not a regular file objectfile toto doesn't existpermission denied for file tmp_parse.txt..----------------------------------------------------------------------Ran 4 tests in 0.065s

OK

11h34
le parser est bien avancé. il lui manque quelques fonctions comme l’impression

c’est facile a ajouter. l’impression sera par defaut basique en utilisant les marker et séparateur par défaut.
touts les attributs seront imprimé par défaut.
il sera cependant possible de spécifier ces paramètres

DANS Items def __str__(self, text_sep='', attr_list=[]):     '''     print x  x.__str__()     print the data in a text format, using default options

     '''     #if no text separator was specified, use the default one     if not text_sep : text_sep = self.text_sep     #if no attribute list was specified, use the default one if exist or nothing     if not attr_list :         try :             attr_list = self.attr_list         except AttributeError :             pass     #if header exist add the header_mark at the top of the each line     #and add all lines into the text to print     try :         text_list = map(lambda x: self.header_mark + x , self.header)     except AttributeError :         text_list = []     #if exist, add the attribute list, join by the text separator     if attr_list : text_list += [join(attr_list,text_sep)]     #call the print method of the item using the text separator and the attribute list as options     #and add the resulting text into the text to print     text_list += [item.__str__(text_sep=text_sep, attr_list=attr_list) for item in self]     #join element of the text to print by a \n and return it     return join(text_list, '\n')

DANS Item def __str__(self, text_sep = '\t', attr_list = []):     '''     print x  x.__str__()     print all element of x using tabulation as default text separator     options :     x.__str__([text_sep, [attr_list]])     print only the values of the attributes specified in attr_list     return the values as text using text_sep as text separator

     '''     #if no attribute list, use all the attributes of self         if not attr_list : attr_list = self.keys()     #get the value of the attributes in attr_list or an empty string if attribute does'nt exist     text_val = [self.get(attr, '') for attr in attr_list]     #return the values in a text format join by text_sep     return join(map(str,text_val), text_sep)

et j’ai donc ajouté un module de test pour cette fonction

 def test_printed_data(self):     '''     teste si les donnees parsees sont imprimees correctement     '''     header_info_1 = "header_01"     header_info_2 = "header_02"     col_1 = "col_1"     col_2 = "col_2"     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 with a complete file     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.__str__() == test_file)

     #test using a specific attribute list     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')     attr_list=[col_1]     test_file='''#%s#%s%s%s%s''' % (header_info_1, header_info_2, col_1, data_l1_c1, data_l2_c1)     self.assert_(test_item.__str__(attr_list=attr_list) == test_file)

     #test using a specific text separator     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')     attr_list=[col_1]     test_file='''#%s#%s%s;%s%s;%s%s;%s''' % (header_info_1, header_info_2, col_1, col_2, data_l1_c1, data_l1_c2, data_l2_c1, data_l2_c2)     self.assert_(test_item.__str__(text_sep=";") == test_file)

     #test with no header     test_file='''%s\t%s%s\t%s%s\t%s''' % (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.__str__() == test_file)

     #test with no attribute list     test_file='''%s\t%s%s\t%s''' % (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.__str__() == test_file)

     #test with no data     test_file='''%s\t%s''' % (col_1, col_2)     self.write_tmp_file(filetext = test_file)     test_item = Items(file_arg = 'tmp_file')     self.assert_(test_item.__str__() == test_file)

     #test with only header     test_file='''#%s#%s''' % (header_info_1, header_info_2)     self.write_tmp_file(filetext = test_file)     test_item = Items(file_arg = 'tmp_file')     self.assert_(test_item.__str__() == test_file)

     #test with nothing     test_file=''     self.write_tmp_file(filetext = test_file)     test_item = Items(file_arg = 'tmp_file')     self.assert_(test_item.__str__() == test_file)

     remove('tmp_file')     def valid_object(self, nb_header_line, nb_col, nb_line):     '''     x.valid_object(nb_header_line, nb_col, nb_line)     create a tmp file with the given parameters     parse the file using Items     test if all data have been correctly parsed     '''     #create the name of the file using the parameters (it could have been an arbitrary filename)     filename = 'file_%s_%s_%s' % (nb_header_line, nb_col, nb_line)     #call the method for creating the file     self.create_tmp_file(filename, nb_header_line, nb_col, nb_line)

     self.valid_item(filename=filename, nb_header_line=nb_header_line, nb_col=nb_col, nb_line=nb_line)

     #delete te tmp file     self.delete_tmp_file(filename)

qui renvoie

...the argument is not a regular file objectfile toto doesn't existpermission denied for file tmp_parse.txt..----------------------------------------------------------------------Ran 5 tests in 0.065s

OK

je vais aussi ajouter une fonction assez pratique au parser.
il s’agit d’une sorte d’index des infos qui sont dans les attributs.
je ne suis pas sur que ce soit pertinenet car ca risque d’etre lourd d’un point de vue traitement.
il faut probablement que ce soit fait a la demande.
comment procéder ?
il faut une fonction de requete.
cette fonction va créer l’index nécessaire s’il n’existe pas et rechercher dedans par la suite.
elle va donc passer par une fonction de création d’index. cette fonction var créer l’index souhaité pour le ou les attributs passé en argument.
en laissant la possibilité d’avoir plusieurs attributs en argument, on s’assure que la fonction pourra etre utilisée dans d’autre cas, comme la création de ces index de manière explicite.
si aucun argument est passé a cette fonction alors l’ensemble des index pour chaque attribut sera créé (risque d’etre long)

def create_indexes(self, attr_list=[]):       '''       x.create_indexes([attr_list])       attr_list -> list of string or str       create indexes for a list of attribute       '''       #test if there is an argument       if not attr_list :           #if not try to get the default attribute list           try :               attr_list = self.attr_list           except AttributeError :               #catch the error if it does'nt exist               print "no attribute for these data, indexes cannot be created"               #and stop function by raising error               raise       #build the index for each item       map(lambda x : self.build_index(x, attr_list), self)

   def build_index(self, item, attr_list):       '''       x.build_index(item, attr_list)       item -> instance of item       attr_list -> list of str or str       '''       #if attr_list is a single string, convert, put it in a list       if type(attr_list) == str : attr_list=[attr_list]       #for each attr       for attr in attr_list :           #test if it's a valid attribute           if attr not in self.attr_list :               print "%s is not a valide attribute" % attr               continue           #try to access to the value associated with the key attr           #if it doesn't exist, create the key as an empty dict and return it           try :               attr_index_list = self.indexes.setdefault(attr,{})           except AttributeError :               #catch the error if the attribute indexes does'nt exist,               #and create it               self.indexes = {}               attr_index_list = self.indexes.setdefault(attr,{})           #get the list associated with the key of the attribute of the item           #or an empty list if the key does'nt already exist           attr_index_item_list = attr_index_list.setdefault(item.get_attribute(attr), [])           #add item into this list           attr_index_item_list.append(item)

   def get_items_by_value(self, attr, value):       '''       item_list = x.get_items_by_value(attr, value)       attr -> str       value -> str, int or float       item_list -> instance of Items       return a list of items having a specified value for the attribute attr       '''

       try :           #check if the attr is in the attribute list           if attr not in self.attr_list :               #if not raise a keyError               print "%s is not a valide attribute" % attr               raise IndexError, "%s is not a valide attribute" % attr       except AttributeError :           #catch the AttributeError if ther is no attribute list           print "no attribute list available"           raise       #try to access to the index of the given attribute       try :           attr_index_list = self.indexes[attr]       except AttributeError :           #catch the error if the attribute indexes doesn't exist           #or if the key attr does'nt exist           #the index for the given attribuite is created           self.create_indexes(attr)           #and put into attr_index_list           attr_index_list = self.indexes[attr]       except KeyError :           #catch the error if the attribute indexes doesn't exist           #or if the key attr does'nt exist           #the index for the given attribuite is created           self.create_indexes(attr)           #and put into attr_index_list           attr_index_list = self.indexes[attr]       #return an instance of items with the list of the items matching with the given value       return Items(item_list = attr_index_list.get(value, []))

j’ai du coup du modifier un peu le __init__ pour qu’il accepte des liste de item pour créer l’objet. cela permet dans une fonction de recherche de renvoyer un sous ensemble de l’objet de la meme forme que l’objet (une instance de items) plutot que une simple liste.
j’ai ajouté le code de test correspondant

def test_indexed_data(self):       '''       teste si les index creer sont correct       '''       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"       data_l2_c3 = "data_c3"       data_l1_c3 = "data_c3"       test_file='''#%s#%s%s\t%s\t%s%s\t%s\t%s%s\t%s\t%s''' % (header_info_1, header_info_2, col_1, col_2, col_3, data_l1_c1, data_l1_c2, data_l1_c3, data_l2_c1, data_l2_c2, data_l2_c3)       self.write_tmp_file(filetext = test_file)       test_item = Items(file_arg = 'tmp_file')       items_index_1 = test_item.get_items_by_value(col_1, data_l1_c1)       self.assert_(isinstance(items_index_1,Items))       self.assert_(len(items_index_1)==1)       self.assert_(isinstance(items_index_1[0], Item))       self.assert_(items_index_1[0].get_attribute(col_1) == data_l1_c1)

       items_index_2 = test_item.get_items_by_value(col_3, data_l1_c3)       self.assert_(isinstance(items_index_2,Items))       self.assert_(len(items_index_2)==2)       self.assert_(isinstance(items_index_2[0], Item))       self.assert_(isinstance(items_index_2[1], Item))       self.assert_(items_index_2[0].get_attribute(col_3) == data_l1_c3)       self.assert_(items_index_2[1].get_attribute(col_3) == data_l2_c3)

       items_index_3 = test_item.get_items_by_value(col_1, 'toto')       self.assert_(isinstance(items_index_3,Items))       self.assert_(len(items_index_3)==0)       self.assertFalse(items_index_3)

       self.assertRaises(IndexError, test_item.get_items_by_value, attr = 'col_0', value = 'toto')     

       remove('tmp_file')

le test renvoie alors

col_0 is not a valide attribute....the argument is not a regular file objectfile toto doesn't existpermission denied for file tmp_parse.txt..----------------------------------------------------------------------Ran 6 tests in 0.170s

OK

Le parser a l’air a peu pret complet, a voir a l’usage.
–> todo list tester le parser sur un plus gros fichier, type gal

Publié dans SganArrayL, conversion, parser, python, type | Laisser un commentaire »

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

avril 9, 2008

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.

Publié dans Uncategorized | Laisser un commentaire »

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

avril 8, 2008

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

Publié dans densité, maximum, sage, segment | Laisser un commentaire »

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

avril 4, 2008

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———————————-

Publié dans Uncategorized | Laisser un commentaire »

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

avril 3, 2008

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.

Publié dans Uncategorized | Laisser un commentaire »

1er avril 2008 – détermination du maximum des segments de sage

avril 1, 2008

les résultats de la manip de sage se présentes sous la formes de séquences dont une partie a été identifiée et localisée sur le génomes.
ces séquences ou tag, sont répartis de façon non aléatoire sur le génome et foirmes des groupes ou cluster.
l’identification de ces cluster a été effectuée par la méthode CNC (Closest Neighbour Clustering). cette méthode met bout a bout chaque brins de tous les chromosomes de la levure. on parcours ensuite l’ensemble du génome ainsi formé jusqu’a ce qu’un sage soit rencontré. on enregistre cette position (S0) puis on parcours a nouveau le génome. si on rencontre a nouveau un sage entre la position sur une distance N depuis la position du sage précédent, on enregistre cette nouvelle position pui on parcours a nouveau le génome sur une distance N et on note si un sage est présent sur cette distance. on repète cette opération tant qu’un sage est trouvé. si aucun sage est trouvé on note la position du dernier sage trouvé (E). on calcule alors entre les position S et E combien de sage il y a et s’il y en a plus que le seuil limite n (a déterminé au préalable) alors on considère que la zone de S à E est un cluster.
puis on recommence la recherche.
cette recherche a donc 2 paramètres, N et n, qui sont très important et font varier les résultats de façon importante.
en prennant N = 100 et n = 1, on obtient 14258 segments.
les résultats sont stockés dans le fichier ~/sage/total/stats/ratio2_5.mat
ils se présente sous la forme d’une matrice avec 11 colonnes (ratio2_5_100_1)
colonnes :
1 -> ratio Wt oligo d(T) / (Wt oligo d(T) + Cbp20-TAP)
2 -> position du start du segment sur le genome reconstitué
3 -> position du end du segment sur le genome reconstitué
4 -> nb de tag dans le segment
5 -> nb de tag dans le segment
6 -> reference du chromosome
7 -> start du segment sur lme chromosome
8 -> end du segment sur le chromosome
9 -> taille du segment
10 -> classe du segment (de 1 à 10 en fn du ratio)

Je dois a partir de cette matrice retrouver la position dans chaque segment pour laquelle le nb de tag est maximale.
j’ai écrit un script simple.
je parcours la matrice ligne par ligne.
je récupére les positions start et end du segment considéré sur le génome linearisé (stocké dans la matriceheader cette matrice a 24millions de lignes et 7 colonnes. ces colonnes correspondent a : 1-> manip 1 à 10; 2-> manip 1 & 2; 3-> manip 1 à 8; 4-> manip 1,2,3,4,9 et 10; 5-> manip 1,2, 7 et 8; 6-> manip 7 & 8; 7-> manip 1 à 12; c’est la colonne 5 qui nous intéresse dans ce cas)
j’extrait le segment du genome
je recherche le max de ce segment et détermine sa position dansle segment
si j’en ai plusieurs, je garde que le premier
je calcule la position du max dans le chromosome en l’ajoutant a celle du strat
je rajoute une colonne dans la matrice ratio2_5_100_1 ou je stocke cette valeur
et voila

for i=1:size(ratio2_5_100_1,1)d= ratio2_5_100_1(i,2);f=ratio2_5_100_1(i,3);segment=header(d:f,5);max_seg=find(segment==max(segment));ratio2_5_100_1(i,11)=ratio2_5_100_1(i,7) + max_seg(1) - 1;end

cette matrice est ensuite expostée en format texte tabulé avec l’extension .csv.mat par la fonction ~/sage/total/stats/save_mat.m
le fichier créé a été sauvegardé dans ~/sage/total/stats/max_segments/ratio2_5_100_1.mat.csv
je l’ai ensuite envoyé sur duanierousseau
je l’ai traité alors avec le programme ~/workspace/genepy/src/.segment.py
tout d’abord j’ai converti le fichier .csv.mat en fichier .csv lisible par excel ou openOffice

def matlab2csv(segment_path):

##===============================================================================##           Conversion fichier matlab en fichier csv##===============================================================================print "loading %s" % segment_pathfhin_mat=open("%s.mat.csv" % segment_path,'r')matlab_txt = convert_mat_file(fhin_mat)fhin_mat.close()fhin_csv=open("%s.csv" % segment_path,'w')print >>fhin_csv,"ratio\tstart_gen\tend_gen\tnb1\tnb2\tgenome_pos\tstart\tend\twidth\tgroup\tmaxi"print >>fhin_csv, matlab_txtfhin_csv.close()

le fichier crée est /home/gim2/sage/total/stats/ratio2_5/max_segments/ratio2_5_100_1.csv
puis j’ai déterminer la position des features mais en me basant non plus comme je l’ai deja fait sur le milieu des segment mais sur son maximum

def fetch_features_around_segment(segment_path):#===============================================================================#  recherche des features environantes des segments#===============================================================================#    segment_path = '%ssegments_location%s' % (dir_region, suffixe)seg_list = Segment_list(file = "%s.csv" % segment_path)print "%s read" % segment_pathdb_path='%s/bank/saccharomyces_cerevisiae.gff' % homedb_features = DBGFF.DB(db_path)fhout_path = '%ssegments_description%s' % (dir_region, suffixe)fhout = open("%s.csv" % fhout_path, 'w')segment_info_list=('id','chr', 'strand','start','end','maxi','nb1', 'ratio', 'group')feature_info_list=('feature_up_sens_name', 'feature_up_sens_gene', 'feature_up_sens_start_to_segment', 'feature_up_sens_end_to_segment',             'feature_down_sens_name', 'feature_down_sens_gene', 'feature_down_sens_start_to_segment', 'feature_down_sens_end_to_segment',             'feature_up_antisens_name', 'feature_up_antisens_gene', 'feature_up_antisens_start_to_segment', 'feature_up_antisens_end_to_segment',             'feature_down_antisens_name', 'feature_down_antisens_gene', 'feature_down_antisens_start_to_segment', 'feature_down_antisens_end_to_segment',             'feature_overlap_sens_name', 'feature_overlap_sens_gene', 'feature_overlap_sens_start_to_segment', 'feature_overlap_sens_end_to_segment',             'feature_overlap_antisens_name', 'feature_overlap_antisens_gene', 'feature_overlap_antisens_start_to_segment', 'feature_overlap_antisens_end_to_segment')print >>fhout, "%s\t%s" % (join(segment_info_list, '\t'),join(feature_info_list, '\t') )               

current_chr = ''no_feature_txt = 'no feature found'for segment in seg_list :    mid = segment.get_attribute('mid')    mid = segment.get_attribute('maxi')    chr = segment.get_attribute('chr')    rev = segment.get_attribute('rev')    if chr != current_chr:        feature_chr = db_features.getDBChr2(chr)        current_chr = chr        feature_w = feature_chr.select('rev', False)        feature_c = feature_chr.select('rev', True)        f_start_w = feature_w.get_list_attr('start')        f_end_w = feature_w.get_list_attr('end')        f_coord_w = zip(f_start_w, f_end_w)        f_w = dict(zip(f_coord_w, feature_w))        f_start_c = feature_c.get_list_attr('start')        f_end_c = feature_c.get_list_attr('end')        f_coord_c = zip(f_start_c, f_end_c)        f_c = dict(zip(f_coord_c, feature_c))        key_w = f_w.keys()        key_w.sort(key=operator.itemgetter(0))

        key_c = f_c.keys()        key_c.sort(key=operator.itemgetter(0))        print chr

    dist_f_start_w = map(lambda x : x - mid, f_start_w)    dist_f_end_w = map(lambda x : x - mid, f_end_w)    dist_f_start_c = map(lambda x : x - mid, f_start_c)    dist_f_end_c = map(lambda x : x - mid, f_end_c)

    dist_f_start_w_up = filter(lambda x : x  0, dist_f_start_w)    dist_f_end_w_down = filter(lambda x : x > 0, dist_f_end_w)    dist_f_start_c_up = filter(lambda x : x  0, dist_f_start_c)    dist_f_end_c_down = filter(lambda x : x > 0, dist_f_end_c)

    f_start_w_up = extreme(dist_f_start_w_up, max)    f_end_w_up = extreme(dist_f_end_w_up, max)    f_start_w_down = extreme(dist_f_start_w_down, min)    f_end_w_down = extreme(dist_f_end_w_down, min)    f_start_c_up = extreme(dist_f_start_c_up, max)    f_end_c_up = extreme(dist_f_end_c_up, max)    f_start_c_down = extreme(dist_f_start_c_down, min)    f_end_c_down = extreme(dist_f_end_c_down, min)

    f_coord_w_up = (f_start_w_up + mid, f_end_w_up + mid)    f_coord_w_down = (f_start_w_down + mid, f_end_w_down + mid)    f_coord_c_up = (f_start_c_up + mid, f_end_c_up + mid)    f_coord_c_down = (f_start_c_down + mid, f_end_c_down + mid)

#                print mid#                print f_coord_w_up#                print f_coord_w_down#                print f_coord_c_up#                print f_coord_c_down

    f_start_w_over = 0    feature_w_up = f_w.get(f_coord_w_up, None)    while not feature_w_up :        if not f_start_w_up :            break        f_start_w_over = f_start_w_up        dist_f_start_w_up.remove(f_start_w_up)        f_start_w_up = extreme(dist_f_start_w_up, max)        f_coord_w_up = (f_start_w_up + mid, f_end_w_up + mid)

        feature_w_up = f_w.get(f_coord_w_up, None)

    f_end_w_over = 0    feature_w_down = f_w.get(f_coord_w_down, None)    while not feature_w_down :        if not f_end_w_down :            break        f_end_w_over = f_end_w_down        dist_f_end_w_down.remove(f_end_w_down)        f_end_w_down = extreme(dist_f_end_w_down, min)        f_coord_w_down = (f_start_w_down + mid, f_end_w_down + mid)        feature_w_down = f_w.get(f_coord_w_down, None)

    f_start_c_over = 0    feature_c_up = f_c.get(f_coord_c_up, None)    while not feature_c_up :        if not f_start_c_up :            break        f_start_c_over = f_start_c_up        dist_f_start_c_up.remove(f_start_c_up)        f_start_c_up = extreme(dist_f_start_c_up, max)        f_coord_c_up = (f_start_c_up + mid, f_end_c_up + mid)

        feature_c_up = f_c.get(f_coord_c_up, None)

    f_end_c_over = 0    feature_c_down = f_c.get(f_coord_c_down, None)    while not feature_c_down :        if not f_end_c_down :            break        f_end_c_over = f_end_c_down        dist_f_end_c_down.remove(f_end_c_down)        f_end_c_down = extreme(dist_f_end_c_down, min)        f_coord_c_down = (f_start_c_down + mid, f_end_c_down + mid)        feature_c_down = f_c.get(f_coord_c_down, None)

    f_coord_w_over = (f_start_w_over + mid, f_end_w_over + mid)    feature_w_over = f_w.get(f_coord_w_over, None)

    f_coord_c_over = (f_start_c_over + mid, f_end_c_over + mid)    feature_c_over = f_c.get(f_coord_c_over, None)

    if rev :        feature_sens_up = feature_c_down        f_dist_sens_up = (f_end_c_down,f_start_c_down)        feature_sens_down = feature_c_up        f_dist_sens_down = (f_end_c_up,f_start_c_up)        feature_anti_up = feature_w_down        f_dist_anti_up = (-f_start_w_down,-f_end_w_down)        feature_anti_down = feature_w_up        f_dist_anti_down = (-f_start_w_up,-f_end_w_up)        feature_sens_over = feature_c_over        f_dist_sens_over = (f_end_c_over,f_start_c_over)        feature_anti_over = feature_w_over        f_dist_anti_over = (-f_start_w_over,-f_end_w_over)

    else :        feature_sens_up = feature_w_up        f_dist_sens_up = (-f_start_w_up,-f_end_w_up)        feature_sens_down = feature_w_down        f_dist_sens_down = (-f_start_w_down,-f_end_w_down)        feature_anti_up = feature_c_up        f_dist_anti_up = (f_end_c_up,f_start_c_up)        feature_anti_down = feature_c_down        f_dist_anti_down = (f_end_c_down,f_start_c_down)        feature_sens_over = feature_w_over        f_dist_sens_over = (-f_start_w_over,-f_end_w_over)        feature_anti_over = feature_c_over        f_dist_anti_over = (f_end_c_over,f_start_c_over)

    fhout.write(segment.__str__(segment_info_list))    feature_info = ()    if feature_sens_up :        feature_sens_up_info = (feature_sens_up.getattribute('name'),feature_sens_up.getattribute('gene')) + f_dist_sens_up    else :        feature_sens_up_info = ('', '', '', '')    feature_info += feature_sens_up_info

    if feature_sens_down :        feature_sens_down_info = (feature_sens_down.getattribute('name'),feature_sens_down.getattribute('gene')) + f_dist_sens_down    else :        feature_sens_down_info = ('', '', '', '')    feature_info += feature_sens_down_info

    if feature_anti_up :        feature_anti_up_info = (feature_anti_up.getattribute('name'),feature_anti_up.getattribute('gene')) + f_dist_anti_up    else :        feature_anti_up_info = ('', '', '', '')    feature_info += feature_anti_up_info

    if feature_anti_down :        feature_anti_down_info = (feature_anti_down.getattribute('name'),feature_anti_down.getattribute('gene')) + f_dist_anti_down    else :        feature_anti_down_info = ('', '', '', '')    feature_info += feature_anti_down_info

    if feature_sens_over :        feature_sens_over_info = (feature_sens_over.getattribute('name'),feature_sens_over.getattribute('gene')) + f_dist_sens_over    else :        feature_sens_over_info = ('', '', '', '')    feature_info += feature_sens_over_info

    if feature_anti_over :        feature_anti_over_info = (feature_anti_over.getattribute('name'),feature_anti_over.getattribute('gene')) + f_dist_anti_over    else :        feature_anti_over_info = ('', '', '', '')    feature_info += feature_anti_over_info

    if feature_sens_over and feature_anti_over :        print mid#                    feature_over_info = zip(feature_sens_over_info, feature_anti_over_info)#                    feature_over_info = tuple(map(lambda x : join(x, ' / ') , feature_over_info))    print >>fhout, join(map(str,feature_info), '\t')fhout.close()

le fichier créé est /home/gim2/sage/total/stats/ratio2_5/max_segments/ segments_description2_5_100_1.csv
le fichier est ouvert avec excel sur bambino et la macro du fichier F:/users/christophe/search_features_A_and_B.xlsest utilisée pour trouver les features A & B et pour classifier les segment
j’obtiens le fichier F:/users/christophe/segment_distance_to_features_using_max.xls
j’exporte ce fichier sur douanierousseau.
je vais comparer la classification obtenue evec celle obtenue pour le calcul des distance avec le milieu des segments.
j’ai donc 2 fichiers
~sage/total/stats/ratio2_5/max_segments/segment_distance_to_features_using_max.xls
~sage/total/stats/ratio2_5/max_segments/segment_distance_to_features_using_mid.xls
j’utilise la fonction grep pour trouver toute les catégories

cut -f42 segments_distance_to_feature_using_mid.txt | sort | uniq > categories

et j’ai dans categories

11010010410811112121314151623456789

j’utilise grep pour compter chacune des lignes appartenant a ces catégories

while read line   do       n=$(cut -f42 segments_distance_to_feature_using_mid.txt | sort | grep -w $line -c)       echo categorie $line    $n segments   done &lt categories

et je récupère pour les deux fichiers

 






catégorie
mid to feature
max to feature


1
2417
2451


2
196
0


3
1334
1388


4
22
0


5
1103
1133


6
240
0


7
2012
2083


8
33
0


9
1535
1571


10
218
0


11
1494
1571


12
24
0


13
1917
2024


14
378
0


15
1225
1271


16
24
0


100
11
137


104
21
193


108
32
162


112
23
274

les résultats sont assez proche mais on vois que certaine catégories ont disparues au profit des catégories non définies àla fin.
je referai le point sur ces catégories pour savoir exactement ce qu’il se passe.

Publié dans Matlab, classification, cluster, localisation, python, sage, segment | Laisser un commentaire »

Xtof78000’s Weblog

Archives pour avril 2008

29 avril 2008 – interface sganarrayl

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

15 avril 2008 – Parser pour le fichier gff

15 avril 2008 – Parser pour le fichier gal

14 avril 2008 – ajout du module de conversion des données

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

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

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

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

1er avril 2008 – détermination du maximum des segments de sage

Pages

Archives

Catégories

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