0
|
1 #!/usr/bin/perl
|
|
2
|
|
3 use strict;
|
|
4 use warnings;
|
|
5
|
|
6 my $num_args = $#ARGV + 1;
|
6
|
7 if ($num_args != 9) {
|
0
|
8 print "Il n'y a pas le bon nombre d'arguments !\n";
|
|
9 exit;
|
|
10 }
|
|
11
|
|
12 # INPUT_
|
|
13 my $matrix_file = $ARGV[0]; # fichier tabulé : une liste d'orthogroupes qui se retrouvent ou non dans les différentes souches
|
|
14 my $species_file = $ARGV[1]; # association de chaque souche à son espèce (fichier tabulé également)
|
|
15 my $annotation = $ARGV[2]; # collection de fichiers tabulés qui contiennent pour chaque gène la ou les catégories de COG associée(s)
|
|
16 my $order = $ARGV[3]; # cette entrée correspond simplement au nom des souches qui sont rentrées dans le même ordre que les fichiers d'annotation : cela permet de savoir pour un fichier COG à quelle souche et donc plus tard à quelle espèce il correspond
|
|
17 my $annotation_GFF = $ARGV[4]; # fichiers avec les GFF
|
6
|
18 # my $order_GFF = $ARGV[5];
|
0
|
19
|
|
20 # OUTPUT_
|
6
|
21 my $output = $ARGV[5]; # liste des espèces avec leurs orthogroupes (présence-absence)
|
|
22 my $output2 = $ARGV[6]; # fichier des moyennes
|
|
23 my $output3 = $ARGV[7]; # fichier de la liste des valeurs pour chaque catégorie de COG et pour chaque espèce
|
|
24 my $output4 = $ARGV[8]; # fichier avec les catégories de COG pour core-génome / génome accessoire / gènes spé
|
0
|
25
|
|
26
|
|
27 # print "ok\n";
|
|
28 # exit;
|
|
29
|
2
|
30 my @list_gff = split(',', $annotation_GFF); # liste des différents fichiers GFF (qui se retrouvent dans le dossier Annotation Maker)
|
0
|
31 my %hSpecies = (); # HASH -> key: N_Id (ex NF_AR12) ; val: nom de l'esp (ex Naegleria Fowleri)
|
|
32
|
|
33 ######################## LE SPECIES_FILE ###########################
|
|
34 open (S, $species_file);
|
|
35 while (my $line = <S>){
|
|
36
|
|
37 $line =~s/\n//g; $line =~s/\r//g;
|
|
38 my @sp = split('\t', $line);
|
|
39 # print "$line\n";
|
|
40 # exit;
|
2
|
41 $hSpecies{$sp[0]} = $sp[1]; # HASH -> key: N_Id ; val: name
|
0
|
42
|
|
43 }
|
|
44 my $nbr = keys (%hSpecies); #compter le nombre de souches max
|
|
45 # = taille de la table de hash
|
|
46 # print "J'ai $nbr clés\n";
|
|
47 # exit;
|
|
48
|
|
49 close (S);
|
|
50
|
|
51 #///////////////////////////////////////////////////////////////////////////////////////////////////
|
|
52
|
|
53 ############################################ LA MATRICE ############################################
|
|
54
|
|
55 open(M, $matrix_file);
|
|
56
|
|
57 my $first_line = <M>;
|
|
58 $first_line =~s/\n//g; $first_line =~s/\r//g; # ne garder que la première ligne du tableau
|
|
59 my @samples = split(/\t/,$first_line); # mettre dans une liste (@samples) chaque intitulé de colonne = N_Id
|
|
60 # print "$first_line\n";
|
|
61 # exit;
|
|
62
|
|
63 # Le but ici est de récupérer les combinaisons associées à chaque espèce : NF, NG et NL
|
|
64 my %hCombination =(); # HASH -> key: N_Id ; val: combinaison
|
|
65
|
|
66 for (my $i=1; $i <= $#samples; $i++){ # on parcourt chaque colonne ($i) mais on ne regarde que le N_Id
|
|
67 my $header = $samples[$i]; # on récupère le N_Id dans $header (soit le nom de la colonne i)
|
|
68 my $species = $hSpecies{$header}; # on regarde dans la table avec N_Id => Nom esp et on attribue à chaque header (qui est ici une clé) sa valeur donc son nom d'esp correspondant
|
|
69 $hCombination{$species} .= "_".$i; # à chaque tour de boucle, pour une $species spé va ajouter le n° de colonne $i pour avoir la combinaison spé à chaque esp
|
|
70 # print "$header\n";
|
|
71 # exit;
|
|
72 }
|
|
73
|
|
74
|
|
75 # foreach my $species (keys (%hCombination)){
|
|
76 # my $combination = $hCombination{$species};
|
|
77 # # print "$species $combination\n";
|
|
78 # }
|
|
79
|
|
80
|
|
81 # exit;
|
|
82
|
|
83 # orthogrp présents :
|
|
84 my %hCombination_prs = (); # HASH -> key: combinaison ; val: liste des orthogroupes
|
|
85 # orthogrp absents :
|
|
86 my %hCombination_abs = (); # idem
|
|
87
|
|
88
|
|
89
|
|
90 my %coregenes = (); # HASH -> key: gene ; val: orthogroupe (pour core-genome)
|
|
91 my %specificgenes = (); # HASH -> key: gene ; val: orthogroupe (pour gènes spécifiques)
|
|
92 my %accessorygenes = (); # HASH -> key: gene ; val: orthogroupe (pour génome accessoire)
|
|
93
|
2
|
94 my $coregene_line;
|
|
95 my %coregenes2 = (); # HASH -> key1: colonne i ; key2: gène ; val: orthogroupe
|
|
96
|
|
97 my %Genes_of_OG = (); # HASH -> key1: orthogroupe ; key2: colonne i ; val: gène
|
|
98
|
0
|
99
|
|
100 while(<M>) {
|
|
101
|
2
|
102 my $line = $_;
|
|
103 $line =~s/\n//g; $line =~s/\r//g;
|
0
|
104 my $nb_found = 0;
|
2
|
105 my @infos = split(/\t/,$line);
|
0
|
106 my $orthogroup = $infos[0]; # on récupère le nom de l'orthogroupe dans $orthogroup
|
|
107 my $first_column = $infos[1]; # ici on récupère les gènes de la première colonne qui vont nous servir pour le core-génome
|
|
108 my $combi_prs = "";
|
|
109 my $combi_abs = "";
|
|
110 my $val;
|
|
111 my $gene_random;
|
|
112
|
2
|
113
|
0
|
114 for (my $i=1; $i <= $#infos; $i++){ # on travaille par ligne puis dans chaque ligne (while(<M>)), cellule par cellule (cette boucle for)
|
|
115
|
|
116 $val = $infos[$i]; # on récupère l'information contenue dans la case $i
|
|
117
|
|
118 if ($val =~/\w/){ # s'il cette cellule contient qq chose...
|
|
119 $combi_prs .= "_".$i; # ...on va concaténer notre chaine $combi_prs pour que cela forme une combinaison
|
|
120 $nb_found++; # on incrémente le compteur qui permet de savoir cb de fois notre orthogroupe est présent (le but sera de l'utiliser quand nb_found == 9)
|
|
121 $gene_random=$val; # on récupère la valeur de la case (les gènes)
|
2
|
122
|
|
123 my @table_genes = split (',', $val);
|
|
124 my $premier_gene = $table_genes[0];
|
|
125 $Genes_of_OG{$i}{$orthogroup} = $premier_gene; # pour chaque orthorgoupe de chaque colonne, on récupère le premier gène
|
0
|
126 }
|
|
127
|
|
128 else { # si jamais il n'y a rien dans la cellule...
|
|
129 $combi_abs .= "_".$i; # ... on fait la même chose mais avec $combi_abs
|
|
130 }
|
|
131
|
|
132 }
|
|
133
|
|
134 # $hCount{$combi}++;
|
|
135 $hCombination_prs{$combi_prs}.=$orthogroup."\n"; # à la fin de chaque ligne, on va ajouter notre orthogroupe à la combinaison qui lui correspond
|
|
136 $hCombination_abs{$combi_abs}.=$orthogroup."\n";
|
|
137
|
|
138
|
|
139
|
|
140 if ($nb_found == $#infos){ # si nb_found = au nombre de souche, c'est qu'on a à faire à un core-génome
|
|
141 # print "$orthogroup\n";
|
|
142 # print "$nb_found\n=================\n";
|
2
|
143 for (my $i=1; $i <= $#infos; $i++){
|
|
144 my @list_of_genes = split (',', $infos[$i]); # ici va séparer tous les gènes (qui se présentent comme une liste, séparés par des ',')
|
|
145 my $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement !
|
|
146 $coregenes{$first_gene}= $orthogroup; # on va récupérer ce premier gène qu'on met dans un hash (pour y avoir accès facilement, d'où val = 1, ici ça n'a pas d'importance)
|
|
147 $coregenes2{$i}{$first_gene}= $orthogroup;
|
|
148
|
|
149
|
|
150
|
|
151
|
|
152 }
|
|
153 if (!$coregene_line){
|
|
154 $coregene_line = $line;
|
|
155 }
|
0
|
156 }
|
|
157 elsif ($nb_found == 1) { # si on a un gène spé
|
|
158 my @list_of_genes = split (',', $gene_random); # idem, on ne veut qu'un seul gène donc on crée la liste
|
|
159 my $first_gene = $list_of_genes[0]; # on ne prend que le premier
|
|
160 $specificgenes{$first_gene}= $orthogroup; # et pareil on crée la table de hash
|
|
161 }
|
|
162
|
|
163 else { # là c'est le génome accessoire, i.e tout le reste !
|
|
164 my @list_of_genes = split (',', $gene_random);
|
|
165 my $first_gene = $list_of_genes[0];
|
|
166 $accessorygenes{$first_gene}= $orthogroup;
|
|
167 }
|
|
168
|
|
169 }
|
2
|
170
|
|
171 my %hCol_Annotated = (); # HASH -> key: colonne ; val: 1 (colonnes pour lesquelles les GFF sont présents)
|
|
172
|
|
173 # Le but ici est de ne garder que les colonnes (donc les souches) qui ont un fichier GFF associé
|
|
174 my @list_column = split ('\t', $coregene_line);
|
|
175 for (my $i=1; $i <= $#list_column; $i++){
|
|
176 my @list_genes = split (', ', $list_column[$i]);
|
|
177 my $premier_gene = $list_genes[0];
|
|
178 my $strain = $samples[$i]; # récupérer le nom de la souche
|
|
179
|
|
180
|
|
181 foreach my $gff (@list_gff){
|
|
182 my $result_grep = `grep $premier_gene $gff`;
|
|
183
|
|
184 if ($result_grep){
|
|
185 $hCol_Annotated{$i}=$strain;
|
|
186
|
|
187 }
|
|
188 # print "$result_grep\n";
|
|
189 }
|
|
190 }
|
|
191 # exit;
|
|
192 # foreach my $i (sort keys (%coregenes2)){ # parcours de la table %hCount2 au niveau des catégories
|
|
193 # foreach my $gene (keys %{$coregenes2{$i} }){ # parcours de la table %hCount2 au niveau des espèces
|
|
194 # print "$i\t$gene\t".$coregenes2{$i}{$gene}."\n";
|
|
195 # }
|
|
196 # }
|
|
197
|
|
198 # while (my ($k,$v) = each(%strain_specie)) {
|
|
199 # print "i=$k strain=$v\n";
|
0
|
200 # }
|
|
201 # exit;
|
|
202 # foreach my $oups (keys (%coregenes)) {
|
|
203 # print "$oups\n";
|
|
204 # }
|
2
|
205 # exit;
|
0
|
206
|
|
207 close (M);
|
|
208
|
2
|
209 my %Hash_Specific = ();
|
|
210
|
0
|
211 open (OUT, '>', $output) or die $!;
|
|
212 print OUT "$annotation\n";
|
|
213 foreach my $species (keys (%hCombination)){ # parcours de la table de hash %hCombination (key: nom esp ; val: combi)
|
|
214 my $combination = $hCombination{$species}; # on récupère dans la variable $combination la valeur de chaque clé {species} (= nom esp) de la table de hash %hCombination
|
|
215 my $ortho_presents = $hCombination_prs{$combination}; # $ortho_presents prend la valeur de chaque clé {combination} (récupérée juste au-dessus) de la table de hash %hCombination
|
|
216 my $ortho_absents = $hCombination_abs{$combination}; # en somme on a 3 combi possibles (_1_2_3_4_5 | _6 | _7_8_9) donc pour ces 3 combi-là, qui sont les clés de %hCombination_prs ou_abs, on va retrouver la liste des orthogroupes qui correspondent
|
|
217
|
|
218 # open (OUT,">results.list.txt");
|
|
219
|
|
220 if ($ortho_presents){
|
|
221 print OUT "> $species - present\n";
|
|
222 print OUT "$ortho_presents\n";
|
2
|
223 my @orthogroups_name = split ('\n', $ortho_presents);
|
|
224 foreach my $ortho (@orthogroups_name){
|
|
225 $Hash_Specific{$ortho} = $species;
|
|
226 }
|
0
|
227 }
|
|
228
|
|
229 if ($ortho_absents){
|
|
230 # open (OUT2,">$species.$combination.absents.list.txt");
|
|
231 print OUT "> $species - absent\n";
|
|
232 print OUT "$ortho_absents\n";
|
|
233 }
|
|
234
|
|
235 # close(OUT2);
|
|
236 }
|
|
237
|
|
238 close(OUT);
|
|
239
|
|
240 #//////////////////////////////////////////////////////////////////////////////////////////////////
|
|
241
|
|
242 ############################################### COG ###############################################
|
|
243
|
|
244 # STEP 1 : CORRESPONDANCE ENTRE LES DIFFERENTS FICHIERS DE COG ET L'ORDRE --------------------------------------------
|
|
245 my @files = split(',', $annotation); # liste des différents fichiers COG (qui se retrouvent dans le dossier Naegleria)
|
|
246 my @list = split(',', $order); # liste de l'ordre des souches
|
5
|
247 #my ($f,$l);
|
0
|
248
|
|
249 my %hCorrespondance = (); #HASH -> key: un fichier COG ; val: un nom de souche (ces 2 données sont entrées en input = $annotation et $order)
|
|
250
|
|
251 # ++++++++++++ parcours de 2 listes en même temps ++++++++++++ #
|
8
|
252 my $l = 1;
|
5
|
253 foreach my $f (@files){
|
|
254 $hCorrespondance{$f} = $list[$l]; # on fait correspondre pour chaque fichier de COG, un nom de souche
|
|
255 $l++;
|
0
|
256 }
|
|
257
|
|
258 # #Affichage du hash
|
5
|
259 # foreach my $f (keys %hCorrespondance){
|
0
|
260 # print $f."=>".$hCorrespondance{$f}."\n"
|
|
261 # }
|
5
|
262 # exit;
|
0
|
263
|
|
264 # STEP 2 : POUR CHAQUE FICHIER DE COG, FAIRE CORRESPONDRE L'ESPECE (ET NON LA SOUCHE) -------------------------------------
|
|
265 my %hCorresp_file_species = (); # HASH -> key: un fichier de COG ; val: une espèce
|
|
266 my %species_names; # HASH -> key: nom d'espèce ; val: 1
|
|
267
|
|
268 foreach my $h (keys (%hCorrespondance)){ # parcours de la table de hash {fichier COG => nom souche}
|
|
269 my $smpl = $hCorrespondance{$h}; # $smpl prend la valeur de la clé (donc d'un nom de souche)
|
|
270 my $espece = $hSpecies{$smpl}; # on regarde la correspondance entre ce $smpl et les nom qu'on a dans notre table de hash %hSpecies (fichier "species.txt") pour avoir le nom de l'espèce dans $espece
|
|
271 $species_names{$espece} = 1; # on garde sous le coude nos nom d'espèce dans cette nouvelle table de hash
|
|
272 $hCorresp_file_species{$h} = $espece; # BUT ATTEINT : on donne pour chaque fichier de COG le nom de l'espèce qui lui correspond
|
|
273 }
|
|
274 # while (my ($k,$v) = each(%hCorresp_file_species)) {
|
|
275 # print "file=$k sp=$v\n";
|
|
276 # }
|
|
277 # exit;
|
|
278
|
|
279 # STEP 3 : COMPTAGE DES CATEGORIES DE COG ------------------------------------------------------------------------------
|
|
280 my %hCount2 = (); # HASH -> key1: catégorie de COG ; key2: espèce associée ; val: comptage
|
|
281
|
|
282 # comptage du core-genome / des gènes spé / du génome accessoire
|
|
283 my %hCore_Count = (); # HASH -> key: catégorie de COG ; val: comptage (ce hash ne sera utilisé que pour le core-genome)
|
|
284 my %hSpecific_Count = (); # HASH -> key: catégorie de COG ; val: comptage
|
|
285 my %hAccessory_Count = (); # HASH -> key: catégorie de COG ; val: comptage
|
|
286
|
|
287 # hash pour récupérer le gène
|
|
288 my %hCore_Cat = (); # HASH -> key: catégorie de COG ; val: gène
|
|
289 my %hAccessory_Cat = (); # HASH -> key: catégorie de COG ; val: gène
|
|
290 my %hSpecific_Cat = (); # HASH -> key: catégorie de COG ; val: gène
|
|
291
|
|
292 # hash pour récupérer le gène
|
|
293 my %hCore_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
|
|
294 my %hAccessory_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
|
|
295 my %hSpecific_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
|
|
296
|
2
|
297 my %Cog_of_gene = (); # HASH -> key: gène ; val: cat de COG
|
|
298 my %Specie_of_gene = (); # HASH -> key: gène ; val: souche
|
0
|
299
|
|
300 foreach my $file(@files){ # parcours de la liste des fichiers
|
|
301 my $esp = $hCorresp_file_species{$file}; # on récupère l'espèce pour chaque fichier de COG dans $esp
|
|
302 # print $esp."\n";
|
|
303 # exit;
|
|
304
|
|
305 my %hCount = (); # HASH -> key: catégorie de COG ; val: comptage
|
|
306
|
|
307
|
|
308 open (A, $file); # on va parcourir maintenant chaque fichier un à un
|
|
309
|
|
310 while (my $line2 = <A>){
|
|
311
|
|
312 $line2 =~s/\n//g; $line2 =~s/\r//g; # on procède ligne par ligne
|
|
313 my @Genes = split('\t', $line2);
|
|
314 my $gene = $Genes[0];
|
|
315 my $first_cat = $Genes[2];
|
|
316 $Cog_of_gene{$gene} = $first_cat;
|
|
317
|
|
318 for (my $j=2; $j <= $#Genes; $j++) {
|
|
319 my $cat = $Genes[$j]; # on récupère la ou les catégorie(s) de COG
|
|
320 $hCount{$cat}++; # pour la catégorie donnée, on incrémente son nb d'occurences
|
|
321
|
|
322 if ($coregenes{$gene}){ # si le $gene fait bien partie du core-genome (donc de notre table de hash %coregenes)
|
|
323 $hCore_Count{$cat}++; # on incrémente le hash
|
|
324 $hCore_Cat{$cat}=$gene; # on récupère le nom du gène
|
|
325 }
|
|
326 if ($accessorygenes{$gene}){ # s'il fait partie des gènes accessoires
|
|
327 $hAccessory_Count{$cat}++;
|
|
328 $hAccessory_Cat{$cat}=$gene;
|
|
329 }
|
|
330 if ($specificgenes{$gene}){ # s'il fait partie des gènes spécifiques
|
|
331 $hSpecific_Count{$cat}++;
|
|
332 $hSpecific_Cat{$cat}=$gene;
|
|
333 }
|
|
334 # $hCount2{$cat}{$esp}++; # TABLE DE HASH AVEC CLES=CAT DE COG + ESPECE VAL=COMPTAGE
|
|
335 }
|
|
336
|
|
337 }
|
|
338 close (A);
|
|
339
|
|
340 # print "$file $esp\n=============\n";
|
|
341 while (my ($k,$v) = each(%hCount)) { # parcours de la table de hash de comptage
|
|
342 # print "cat=$k nb=$v\n";
|
|
343 $hCount2{$k}{$esp}.= "$v,"; # pour un $k (= une catégorie de COG) on lui associe son espèce et on donne la valeur du comptage qui vient de %hCount
|
|
344 # le but ici est en fait pour une espèce et une catégorie données on veut le nombre d'occurences par souche (pour NF par ex on aura 5 valeurs car il y a 5 souches)
|
|
345 }
|
|
346
|
|
347 # Récupérer les gènes du core-génome
|
|
348 while (my ($cat_core,$gene_core) = each(%hCore_Cat)) {
|
|
349 $hCore_Cat_Esp{$cat_core}{$esp}=$gene_core;
|
|
350 }
|
|
351 # Récupérer les gènes du génome-accessoire
|
|
352 while (my ($cat_acc,$gene_acc) = each(%hAccessory_Cat)) {
|
|
353 $hAccessory_Cat_Esp{$cat_acc}{$esp}=$gene_acc;
|
|
354 }
|
|
355 # Récupérer les gènes spécifique
|
|
356 while (my ($cat_spe,$gene_spe) = each(%hSpecific_Cat)) {
|
|
357 $hSpecific_Cat_Esp{$cat_spe}{$esp}=$gene_spe;
|
|
358 }
|
|
359
|
|
360 }
|
|
361 # foreach my $category (sort keys (%hSpecific_Cat_Esp)) { # parcours au niveau de la 1ere clé
|
|
362
|
|
363 # foreach my $especeee (keys %{$hSpecific_Cat_Esp{$category} }) { # parcours au niveau de la 2e clé pour la $category donnée
|
|
364
|
|
365 # print "$category\t$especeee\t$hSpecific_Cat_Esp{$category}{$especeee}\n"; # on crée une sortie qui affiche en somme notre hash %hCount2
|
|
366 # }
|
|
367 # }
|
|
368 # exit;
|
|
369
|
|
370 # STEP 4 : AFFICHAGE DANS LE FICHIER DE SORTIE ------------------------------------------------------------------------------
|
|
371 open (OUT4, ">$output4") or die $!;
|
|
372
|
|
373 print OUT4 "Species"."\t"."COG categories"."\t"."Core-genome"."\t"."Accessory genome"."\t"."Specific genes"."\n";
|
|
374
|
|
375 foreach my $category (sort keys (%hCount2)){ # parcours de la table %hCount2 au niveau des catégories
|
|
376 foreach my $especeee (keys %{$hCount2{$category} }){ # parcours de la table %hCount2 au niveau des espèces
|
|
377 print OUT4 "$especeee\t$category\t"; # affichage des esp puis des cat
|
|
378
|
|
379 # if ($hCore_Cat_Esp{$category}{$especeee}) {
|
|
380 # print OUT4 "$hCore_Cat_Esp{$category}{$especeee}\t";
|
|
381 # }
|
|
382 my $c = 0;
|
|
383 if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome
|
|
384 $c = ($hCore_Count{$category}/scalar keys (%coregenes))*100; # calcul du % du comptage
|
|
385 }
|
|
386 print OUT4 "$c\t"; # affichage du %
|
|
387
|
|
388 # if ($hAccessory_Cat_Esp{$category}{$especeee}) {
|
|
389 # print OUT4 "$hAccessory_Cat_Esp{$category}{$especeee}\t";
|
|
390 # }
|
|
391 my $acc = 0;
|
|
392 if ($hAccessory_Count{$category}){ # si cette catégorie existe dans le génome accessoire
|
|
393 $acc = ($hAccessory_Count{$category}/scalar keys (%accessorygenes))*100; # calcul du % du comptage
|
|
394 }
|
|
395 print OUT4 "$acc\t"; # affichage du %
|
|
396
|
|
397 # if ($hSpecific_Cat_Esp{$category}{$especeee}) {
|
|
398 # print OUT4 "$hSpecific_Cat_Esp{$category}{$especeee}\t";
|
|
399 # }
|
|
400 my $s = 0;
|
|
401 if ($hSpecific_Count{$category}){ # si cette catégorie existe dans les gènes spécifiques
|
|
402 $s = ($hSpecific_Count{$category}/scalar keys (%specificgenes))*100; # calcul du % du comptage
|
|
403 }
|
|
404 print OUT4 "$s\n"; # affichage du %
|
|
405 }
|
|
406 }
|
|
407 close (OUT4);
|
|
408
|
|
409 open (OUT3, ">$output3") or die $!;
|
|
410 foreach my $category (sort keys (%hCount2)) { # parcours au niveau de la 1ere clé
|
|
411
|
|
412 foreach my $especeee (keys %{$hCount2{$category} }) { # parcours au niveau de la 2e clé pour la $category donnée
|
|
413
|
|
414 print OUT3 "$category\t$especeee\t$hCount2{$category}{$especeee}\n"; # on crée une sortie qui affiche en somme notre hash %hCount2
|
|
415 }
|
|
416 }
|
|
417
|
|
418 close (OUT3);
|
|
419
|
|
420
|
|
421 open (OUT2, ">$output2") or die $!;
|
|
422
|
|
423 print OUT2 "category";
|
|
424 foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces...
|
|
425 print OUT2 "\t".$e; #... où on récupère le nom de celles-ci
|
|
426 }
|
|
427 print OUT2 "\n";
|
|
428
|
|
429 foreach my $category (sort keys (%hCount2)) { # on parcourt de nouveau les catégories de notre hash à 2 clés
|
|
430 print OUT2 $category;
|
|
431
|
|
432 foreach my $especes (sort keys (%species_names)) { # on parcourt également le hash d'espèces
|
|
433
|
|
434 my $nbr = 0;
|
|
435 if ($hCount2{$category}{$especes}) { # si pour une catégorie et une espèce données, on a un nombre : $nbr prend la valeur de ce dernier
|
|
436 $nbr = $hCount2{$category}{$especes};
|
|
437 }
|
|
438 # $nbr =~s/\n//g; $nbr =~s/\r//g;
|
|
439
|
|
440
|
|
441 my @liste = split(',', $nbr); # vu qu'il peut y avoir plusieurs nombres on les dissocie
|
|
442
|
|
443 my $somme=0;
|
|
444 my $n=0;
|
|
445 my $moyenne=0;
|
|
446 #print "\nma liste de $nbr: ".join("%",@liste)."\n";
|
|
447 foreach my $x (@liste) { # on parcourt nos nombres
|
|
448 $somme=$somme+$x;
|
|
449 $n=$n+1;
|
|
450 }
|
|
451
|
|
452 if ($n>0){
|
|
453 $moyenne = $somme/$n; # on fait le calcul de la moyenne
|
|
454 }
|
|
455 # print "$category, $especes: $hCount2{$category}{$especes}\t";
|
|
456 # print "moyenne = $moyenne\n=============\n";
|
|
457
|
|
458 print OUT2 "\t".$moyenne; # fichier de sortie
|
|
459 }
|
|
460 print OUT2 "\n";
|
|
461 }
|
|
462
|
|
463 close (OUT2);
|
|
464
|
|
465 # foreach my $cat (keys (%hCore_Cat)){
|
|
466 # print OUT4 $c_gene."\t";
|
|
467 # }
|
|
468
|
|
469
|
|
470 #//////////////////////////////////////////////////////////////////////////////////////////////////
|
|
471
|
|
472 ############################################### GFF ###############################################
|
|
473
|
2
|
474
|
6
|
475 # my @order_gff = split(',', $order_GFF); # liste de l'ordre des souches
|
0
|
476 my ($g,$o);
|
|
477
|
|
478 my %hgff_order = (); #HASH -> key: un fichier GFF ; val: un nom de souche (ces 2 données sont entrées en input = $annotation_GFF et $order_GFF)
|
|
479 my %Gene_position = ();
|
|
480 my %Cat_genes = ();
|
|
481
|
2
|
482 my %hash_of_genes = ();
|
0
|
483
|
6
|
484
|
0
|
485 foreach $g (@list_gff){
|
2
|
486 # print "$g\n";
|
6
|
487 # $hgff_order{$g} = $order_gff[$o++]; # on fait correspondre pour chaque fichier GFF, un nom de souche
|
0
|
488 open (G, $g);
|
|
489 while (<G>) {
|
|
490 my @table_gff = split (/\t/, $_);
|
|
491 my $chr = $table_gff[0];
|
|
492 my $start = $table_gff[3];
|
|
493 my $end = $table_gff[4];
|
|
494 my $gene_name = $table_gff[8];
|
|
495 my $type = $table_gff[2];
|
|
496
|
|
497
|
|
498
|
10
|
499 if ($type && ($type eq "mRNA" or $type eq "CDS") && $gene_name =~ /ID=([^;]+);/){
|
0
|
500 my $gene = $1;
|
|
501 # print $gene."\n";
|
|
502 # exit;
|
2
|
503 $hash_of_genes{$gene}=1;
|
|
504
|
0
|
505 foreach my $cog (keys (%hCore_Cat)){
|
|
506 if ($hCore_Cat{$cog} eq $gene){
|
|
507 $Cat_genes{$gene}=$cog;
|
|
508 }
|
|
509 }
|
|
510
|
|
511 $Gene_position{$gene}="$chr\t$start\t$end";
|
|
512 }
|
|
513
|
2
|
514 # foreach my $gene (keys (%hash_of_genes)){
|
|
515 # my $orthogrp = $hGene_OG{$gene};
|
|
516 # print "$orthogrp\n";
|
|
517 # }
|
0
|
518 }
|
|
519
|
|
520 close (G);
|
|
521 }
|
|
522
|
4
|
523 my %Hash_Convert = ( "A"=>1, "B"=>2, "C"=>3, "D"=>4, "E"=>5, "F"=>6, "G"=>7, "H"=>8, "I"=>9, "J"=>10, "K"=>11, "L"=>12, "M"=>13, "N"=>14, "O"=>15, "P"=>16, "Q"=>17, "R"=>18,"S"=>19, "T"=>20, "U"=>21, "V"=>22, "W"=>23, "X"=>24, "Y"=>25, "Z"=>26, "unknown"=>27);
|
|
524
|
2
|
525 mkdir("Core");
|
|
526 foreach my $i (keys (%coregenes2)){
|
0
|
527
|
2
|
528 if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
|
|
529 next;
|
|
530 }
|
|
531
|
|
532
|
|
533 my $strain_name = $hCol_Annotated{$i};
|
|
534
|
|
535 my $specie_name = $hSpecies{$strain_name};
|
0
|
536
|
|
537
|
|
538
|
2
|
539 open (OUT5, "> Core/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
|
8
|
540 print OUT5 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n";
|
0
|
541
|
2
|
542 my $refcoregenes2 = $coregenes2{$i};
|
|
543 my %subhash = %$refcoregenes2;
|
|
544 foreach my $gene (keys (%subhash)){
|
|
545 # print "$gene\n";
|
|
546 my $cat = "unknown";
|
|
547 if ($Cog_of_gene{$gene}){
|
|
548 $cat = $Cog_of_gene{$gene};
|
|
549 }
|
|
550 # if (!$Gene_position{$gene}){
|
|
551 # print "$gene\n coucou"; exit;
|
|
552 # }
|
|
553
|
|
554 # if (!$subhash{$gene}){
|
|
555 # print "$gene\n";
|
|
556 # }
|
4
|
557 print OUT5 $subhash{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n";
|
2
|
558
|
|
559 }
|
|
560
|
|
561 close (OUT5);
|
0
|
562 }
|
|
563
|
2
|
564
|
|
565 mkdir("GroupSpecific");
|
|
566 foreach my $i (keys (%Genes_of_OG)){
|
|
567 if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
|
|
568 next;
|
|
569 }
|
|
570
|
|
571 my $strain_name = $hCol_Annotated{$i};
|
|
572
|
|
573 my $specie_name = $hSpecies{$strain_name};
|
|
574
|
|
575 open (OUT6, "> GroupSpecific/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
|
8
|
576 print OUT6 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n";
|
2
|
577
|
|
578 my $refGenes_of_OG = $Genes_of_OG{$i};
|
|
579 my %subhash = %$refGenes_of_OG;
|
|
580
|
|
581 foreach my $orthogroup (keys (%subhash)){
|
|
582 if ($Hash_Specific{$orthogroup} && $Hash_Specific{$orthogroup} eq $specie_name){
|
|
583 my $gene = $subhash{$orthogroup};
|
|
584
|
|
585 my $cat = "unknown";
|
|
586 if ($Cog_of_gene{$gene}){
|
|
587 $cat = $Cog_of_gene{$gene};
|
|
588 }
|
4
|
589 print OUT6 $orthogroup."\t".$subhash{$orthogroup}."\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n";
|
2
|
590 }
|
|
591
|
|
592 }
|
|
593 close (OUT6);
|
|
594 }
|