20 mai 2008 – analyse des résultat de sage – restriction WT et double mutant (suite)

10h40
après avoir obtenu les segment, je vais donc chercher a les localiser par rapport aux orf environnantes.
Tout d’abord j’exporte la liste des segment sur douanierousseau
scp berthe:sage/total/stats/ratio2_5/WT&mutant/segments_WT_CUTs ~/sage/total/stats/ratio2_5/WT&mutant/
ensuite j’utilise segment.py (dans le workspace genepy) pour transformer ce fichier du format matlab au format csv puis pour localiser les orf
j’ai modifier légèrement segment.py pour que les nom des colonnes correspondent bien a ce qui est dedans.

j’obtiens le fichier segments_WT_CUTs_density_max.csv que j’exporte sur bambino pour appliquer la macro excel de search_features_A_and_B.xls dessus
j’enregistre les résultats dans segments_WT_CUTs_density_max.xls et j’exporte dans le fichier segments_WT_CUTs_dist_to_feature.txt
j’élimine les colonne de texte et j’enregistre en csv
puis j’exporte sur berthemorisot et je l’importe dans matlab

j’obtiens la matrice seg_WT_CUTs_dist2feature et la matrice colheader_seg_WT_CUTs_dist2feature qui contient les intitulés des colonnes.
je représentes les nb de segment de chaque type en fontion de leur distance respectives avec soit le debut ou la fin des orf
j’ai écrit un script qui traces les courbes de façon automatique en
mettant la plupart des bon parametres

function [distance_hist]=plot_dist2feature(data, n, smooth_window)

extreme=cell(1,2);
extreme{1,1}='start';
extreme{1,2}='stop';

sens=cell(1,2);
sens{1,1}='sens';
sens{1,2}='antisens';

classe=cell(2,2);
classe{1,1}='Oligo-dT';
classe{2,1}=8:10;
classe{1,2}='CUTs';
classe{2,2}=1:2;

cols=cell(4,3);

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

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

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

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

max_plot=0;
min_plot=0;
distance_hist=cell(8,5);
i=1;
for id_classe=1:2
for id_extreme=1:2
for id_sens=1:2
for j=1:4
feature_extreme=extreme{1,id_extreme};
feature_sens=sens{1,id_sens};
if strcmp(cols{j,1}, feature_extreme) && strcmp(cols{j,2}, feature_sens)
columns=cols{j,3};
end
end

distance_hist{i,1}=mk_distance_hist(data, classe{2,id_classe}, columns, -1000, 1000, n, 20);
distance_hist{i,2}=feature_extreme;
distance_hist{i,3}=feature_sens;
distance_hist{i,4}=classe{1,id_classe};
distance_hist{i,5}=n;
h=zeros(2001,1);
tmp = distance_hist{i,1};
h(min(tmp)+1001:max(tmp)+1001)=hist(tmp,max(tmp)-min(tmp)+1)';
distance_hist{i,6}=h;
distance_hist{i,7}=smooth(h, smooth_window, 0);
i=i+1;
end
end
end
zero_line=[0 0];
figure
set(gcf,'Color','k');
plot_ids = 1:2:7;
for i=1:4
if i max_plot
max_plot= limite(2);
end
end

for i=1:4
subplot(2,2,i)
ylim([min_plot max_plot]);
line(xlim, zero_line, 'Color', 'w');
line(zero_line, ylim, 'Color', 'w');
end


j’ai rajouté qq éléments sur le graphique.
ce qui me donne ca pour n=4 :

quelques chiffres sur ces données
on a sélectioné uniquement les tag provennant de la souche sauvage (oligo-dT) et du double mutant Trf4/Rrp6 (CUTs)
on a :

oligo-dT -> 49284 tags
Cuts -> 66884 tags

environ 1/3 des tags des Cuts a été éliminé de façon aléatoire afin d’avoir le meme nombre pour les 2 catégories.

on a alors

oligo-dT -> 49284 tags
Cuts -> 49077 tags

il y a au total 98361 tags

la détection des segments sur ces tags
permet de trouver 12021 clusters
ces clusters contiennent au moins 1 tag
si on augmente le nb minimal de tag on a :

1-> 12021
2-> 6001
3-> 4047
4-> 3069
5-> 2482
6-> 2072
7-> 1791
8-> 1572
9-> 1410
10-> 1277

le calcul du ratio Oligo-dT / (CUTs + Oligo-dT) puis le fractionnement des clusters dans 10 classes en fonction de ce ratio donne

Nb tag minimum	Classes -> 1	2	3	4	5	6	7	8	9	10
1	6168	179	94	178	330	60	161	216	200	4435
2	2573	179	94	178	330	60	161	216	200	2010
3	1568	179	94	178	97	60	161	216	200	1294
4	1169	179	94	62	97	60	59	216	200	933
5	935	179	46	62	58	60	59	148	200	735
6	787	140	46	42	58	30	59	107	200	603
7	691	117	46	27	40	30	46	107	176	511
8	621	100	37	27	32	20	46	97	152	440
9	572	85	28	24	29	20	41	85	137	389
10	538	73	19	18	24	16	37	74	121	357

on va considerer que les CUTs sont les segments des classes 1 et 2 et les oligo-dT ceux des classes 8,9 et 10. de plus on ne va s’intéresser qu’aux segments ayant au moins 4 tags
ce qui fait
CUTs -> 1348 segments
oligo-dT -> 1349 segments

dans ces segments on a
38 % (508 segments) des oligo-dT sont entre 600 et 0 avant le start des orf en sens
88 % (1198 segments) des oligo-dT sont entre 0 et 400 après le stop des orf en sens
48 % (646 segments) des oligo-dT sont entre -100 et 400 par rapport au stop des orf en antisens
45 % (605 segments) des CUTs sont entre -600 et 400 par rapport au start des orf en sens
45 % (613 segments) des CUTs sont entre -1000 et -200 par rapport au start des orf en antisens
40 % (538 segments) des CUTs sont entre -600 et 100 par rapport au stop des orf en antisens

15 mai 2008 – analyse des résultat de sage – restriction WT et double mutant

10h25
Je vais reprendre l’analyse des sage a la base en ne prennant en compte que les tag obtenu sur les manip wild type et double mutant Trf4/Rrp6

Je dois refaire le clustering afin de n’avoir que des segment correspondant a des groupes de tag de ces manip.

1- ouverture dans matlab du fichier
~/sage/total/stats/tag_tot_genome.mat

j’ai dans ce fichier les matrice correspondant aux 12 manip de sage (header01 à header12)
chaque matrice à une taille de 24 million ce qui corresponde au nombre de nucleotide de l’ensemble du génome mis bout a bout en séparant les brin C et W
les matrice chr_len et chr_ref donnent les positions de chaque chromosome dans cette matrice

2- la recherche des segments (cluster) va se faire uniquement sur les manip 1 et 8, correspondant au Wild Type (WT) et au double mutant Trf4/Rrp6

3- j’utilise le programme segment.m avec la commande
header = header01 + header08;
header(:,2) = header01;
[header_region] = segment(header, dt, nt, chr_ref);
je récupère la matrice ratio2_1_100_1 qui a été crée
et je représente les histogramme des classe de ratio

figurefor i=1:9subplot(3,3,i)hist(ratio2_1_100_1(find(ratio2_1_100_1(:,4)>i-1),1))xlabel('Ratio Oligo-dT / (CUT + Oligo-dT)', 'color',[1 1 1]);ylabel('Number of clusters', 'color',[1 1 1]);title(strcat('n = ', num2str(i)), 'color',[1 1 1])set(gca, 'Xcolor', [1 1 1], 'Ycolor', [1 1 1])end

et j’obtiens la figure suivante

on voit qu’il y a un nb plus important de cluster de la classe 1 que ceux de la classe 10
si on regarde le nb de tag, on constate aussi qu’il y a environ 20 % de plus de CUTSque d’Oligo-dT

On va donc chercher a réduire artificiellement le nb de CUT en en éliminant de façon aléatoire jusqu’à avoir le meme nombre
script : delete_rand_tag

h8 = header08;tag = find(h8);s = size(tag,1);while sum(h8) > 49284p = floor(rand(1000,1) * s);ptag = tag(p);h8(ptag) = 0;end

header_rand = header01 + h8;
header_rand(:,2) = header01;
je fais ensuite [header_region_rand] = segment(header_rand, dt, nt, chr_ref);
je renome ratio2_1_100_1 en ratio2_1_100_1_ rand
puis je représente les memes histogrammes

figurefor i=1:9subplot(3,3,i)hist(ratio2_1_100_1_rand(find(ratio2_1_100_1_rand(:,4)>i-1),1))xlabel('Ratio Oligo-dT / (CUT + Oligo-dT)', 'color',[1 1 1]);ylabel('Number of clusters', 'color',[1 1 1]);title(strcat('n = ', num2str(i)), 'color',[1 1 1])set(gca, 'Xcolor', [1 1 1], 'Ycolor', [1 1 1])end

et j’obtiens

on voit que c’est plus homogène entre les classes extrèmes.
pour vérifier cela, on cherche a voir si la répartition du nombre de tag dans les segment est la meme pour toute les classes.

j’extrait de la matrice ratio2_1_100_1_ rand les lignes pour lesquelles le ratio est faible (inf a 0.2 CUTs) et fort (sup à 0.8 WT) et je fait l’histogramme du nombre de tag (colonne 4). je choisi pour ces histogramme d’avaoir comme nombre de catégorie la valeur maximale de tag pour chaque type (3318 pour le WT et 1522 pour les CUTs.
puis je trace les courbe correspondantes en échelle log

h_0=hist(ratio2_1_100_1_rand(find(ratio2_1_100_1_rand(:,1)0.8),4), max(ratio2_1_100_1_rand(find(ratio2_1_100_1_rand(:,1)>0.8),4)));
figure;
plot(h_0,'r');
hold on
plot(h_1,'g');

on voit que les segments avec un faible nb de tag sont plus nombreux pour les CUTs. cependant a partir de 4 tags par segment on observe la meme répartition.
on peut donc dire que a partir de 4 tags par segment on a un nombre similaire de segments avec le meme nombre de tag dans ceux-ci.

Je me suis rendu compte que les position des tags sont fausses. en fait elles sont décalée par rapport a ce qu’on veux car on s’interresse a la fin des tag et non au début
j’ai fait un programme qui corrige ca (tag_position_correction)
je corrige donc les 2 matrice contenant les résultats des WT et CUTs

chr_WT=tag_position_correction(header01,chr_len);chr_CUTs=tag_position_correction(h8,chr_len);chr_WT_CUTs = chr_WT + chr_CUTs;chr_WT_CUTs(:,2)=chr_WT;[regionWT_CUTs] = segment(chr_WT_CUTs, dt, nt, chr_ref);

cpt =

       12021

et je récupère la matrice segments_WT_CUTs
je vais maintenant rechercher le max de densité de chaque segment
je calcul d’abord la matrice de densité pour les tag WT et CUts
density_WT_CUTs = densite_sage(chr_WT_CUTs(:,1), 100);
puis je cherche le max de densité
segments_WT_CUTs=find_max_segment(density_WT_CUTs, segments_WT_CUTs);
j’exporte avec save(strcat('segments_WT_CUTs','.mat.csv'),'segments_WT_CUTs', '-ascii','-tabs');

TODO -> export sur douanierousseau, localiser les segment par rapport aux orf et faire les graph

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

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

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

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.

26 mars 2008 – répartition aleatoire de sage sur le génome

17h18 :
Je veux savoir si une répartition des sage aléatoire sur le génome produirait les memes résultats en terme de segments.
l’idée est de répartir sur les 24 million de position possible les 267000 tags de façon aleatoire, puis de rechercher les segments de la meme manière qu’avec les vraies données et ensuite de regarder la tailles des segments produites.
1- je pars des valeurs suivante
genome : 24 313 356 bases
tag, après filtration (juste tag unique et correct) 267023 tag
2- dans matlab je fait une répartition random de ces tags.
pour cela je genere 267023 valeurs aléatoire que je multiplie par la taille du génome
j’ai alors 267023 valeurs comprise entre 0 et 24313356

rand_tag=zeros(267023,1);for i=1:267023rand_tag(i,1)=gen*rand();end

on reparti dans une matrice de la taille du genome ces différents tag

g=zeros(24313356,1);for i=1:267023g(ceil(rand_tag(i,1)))=g(ceil(rand_tag(i,1))) + 1;end

je fait appelle ensuite à la fonction region_hit

cpt = region_hit(g, 1,1,1,1);

je recupère en sortie une matrice
–> cette méthode n’est pas la bonne car je vais etudier la répartition des segments sur une ditribution aléatoire. hors ce n’est pas exactement ce que je veux montrer.
je veux finalement montrer que c’est la distrivution des sage (ou le tirage) qui est ou pas aléatoire.
je suis parti sur l’idée d’un tirage aléatoire.
je ne le fait que sur les tag des manip 1 et 2.
ca en fait 77000
je fait un randperm sur une matrice de la taille du génome et je sélection les 77000 premie résultats.
ce seront les 77000 position aléatoire de mes tag

>> taille_gen=24313356;>> gen=[1:taille_gen];>> x=randperm(taille_gen);>> rand_tag_12=gen(x(1:77212));>> rand_tag_12=sort(rand_tag_12');>> dist_rand_tag12=diff(rand_tag_12);>> hist(dist_rand_tag12,100)

et je récupère en sortie

je fait la meme chose sur les tag des manip 1 & 2
j’utilise le fichier matlab ~/total/stats/tag_tot_genome.mat

n12 = sum(header01) + sum(header02);tag12=zeros(n12,1);j=1;for i=1:size(header01,1)t=header01(i,1);if (t>0)  for k=1:t     tag12(j,1)=i;     j=j+1;  endendendfor i=1:size(header02,1)t=header02(i,1);if (t>0)  for k=1:t     tag12(j,1)=i;     j=j+1;  endendendtag12=sort(tag12);dist_tag12=diff(tag12);hist(dist_tag12(find(dist_tag12<3500)),100)

et je récupère

on voit bien que les distributions sont différentes mais il faudrait une valeur numérique pour le prouver.
on a pour la distribution des tag 1 et 2 :

>>mean(dist_tag12)ans =306.8112>> median(dist_tag12)ans =0>> std(dist_tag12)ans =1.4071e+03

et pour la dsitribution aléatoire

>> mean(dist_rand_tag12)ans =314.8891>> median(dist_rand_tag12)ans =218>> std(dist_rand_tag12)ans =316.0702

on voit bien que la mediane est très différente.
mais ce n’est peut-etre pas un critère suffisant
le test du chi2 permet de voir si une distribution est aléatoire ou non
prennons la taille du génome totale -> 24313356
prennons le nb de tag 1 et 2 -> 77212
s’il sont répartis de façon homgène sur le génome ca fait 1 tag tout les 314 nucleotide environ.
imaginons qu’il soient réparti selon une distribution normal centrée autour de cette distance moyenne.

>> g2=randn(77212,1);>> g2=g2+abs(min(g2));>> mean(g2)ans =4.3959>> coeff = 314/4.3954coeff =71.4383>> g2=g2*coeff;>> hist(g2,100)>> [h,p]=chi2gof(g2)h =0p =0.7848

la distribution de ces distances est la suivante


avec une p-value sur le test de chi2 de 0.7848,
donc une proba très forte que cette distribution soit aléatoire.
Si je représente la distribution des distance entre tag des manip 1 et 2
(voir au-dessus), on voit que celle-ci est centrée sur 0, donc la plupart
des tag sont très proches entre eux

le test du chi2 sur cette distribution donne :

>> [h,p]=chi2gof(dist_tag12, 'nbins',32)h =1p =0

une p-value de 0, soit une proba nulle que ce soit une distribution aléatoire.
on calcule pour la distrib aléatoire

>> median(g2)ans =  314.8353>> mean(g2)ans =  314.8909>> std(g2)ans =   71.6086

et pour les manip

>>mean(dist_tag12)ans =306.8112>> median(dist_tag12)ans =0>> std(dist_tag12)ans =1.4071e+03

je fais ensuite la meme chose pour les tag des manip 7 et 8

n78 = sum(header07) + sum(header08);tag78=zeros(n78,1);j=1;for i=1:size(header07,1) t=header07(i,1); if (t>0)     for k=1:t         tag78(j,1)=i;         j=j+1;     end endendfor i=1:size(header08,1) t=header08(i,1); if (t>0)     for k=1:t         tag78(j,1)=i;         j=j+1;     end endendtag78=sort(tag78);dist_tag78=diff(tag78);

puis je fais une distribution aléatoire de tag équivalente

g=randn(69549,1);>> g=g+abs(min(g));>> mean(g)ans = 4.3684>> 24313356/69549ans =349.5860>> coeff = 349/4.3684coeff =79.8920>> g=g*coeff;>> hist(g,100)

j’obtiens pour la distribution aléatoire

>> [h,p]=chi2gof(g)h =   0p =  0.7109

soit une proba de 0.7 que ce soit une distrib aléatoire
avec la distrib suivante

et pour les tag

>> [h,p]=chi2gof(dist_tag78, 'nbins',25)h =   1p =   0

soit une proba nulle que ce soit aléatoire
avec la distrib suivante

pour cette distribution on a

>> mean(dist_tag78)ans =  348.9803>> median(dist_tag78)ans =     0>> std(dist_tag78)ans =   1.5912e+03

et pour la distrib aléatoire

>> median(g)ans =  348.9698>> mean(g)ans =  348.9981>> std(g)ans =   80.0903

les images sont dans ~/Desktop/screenshot/tags