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