13
|
1 #!/usr/bin/perl
|
|
2
|
|
3 use strict;
|
|
4 use warnings;
|
|
5
|
|
6 my $num_args = $#ARGV + 1;
|
|
7 if ($num_args != 11) {
|
|
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
|
|
18 # my $order_GFF = $ARGV[5];
|
|
19
|
|
20 # OUTPUT_
|
|
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é
|
|
25 my $output7 = $ARGV[9];
|
|
26 my $output8 = $ARGV[10];
|
|
27 # my $output9 = $ARGV[11];
|
|
28
|
|
29
|
|
30 #print "ok\n";
|
|
31 #exit;
|
|
32
|
|
33 my @list_gff = split(',', $annotation_GFF); # liste des différents fichiers GFF (qui se retrouvent dans le dossier Annotation Maker)
|
|
34 my %hSpecies = (); # HASH -> key: N_Id (ex NF_AR12) ; val: nom de l'esp (ex Naegleria Fowleri)
|
|
35
|
|
36 ######################## LE SPECIES_FILE ###########################
|
|
37 open (S, $species_file);
|
|
38 while (my $line = <S>){
|
|
39
|
|
40 $line =~s/\n//g; $line =~s/\r//g;
|
|
41 my @sp = split('\t', $line);
|
|
42 # print "$line\n";
|
|
43 # exit;
|
|
44 $hSpecies{$sp[0]} = $sp[1]; # HASH -> key: N_Id ; val: name
|
|
45
|
|
46 }
|
|
47 my $nbr = keys (%hSpecies); #compter le nombre de souches max
|
|
48 # = taille de la table de hash
|
|
49 # print "J'ai $nbr clés\n";
|
|
50 # exit;
|
|
51
|
|
52 close (S);
|
|
53
|
|
54 #///////////////////////////////////////////////////////////////////////////////////////////////////
|
|
55
|
|
56 ############################################ LA MATRICE ############################################
|
|
57
|
|
58 open(M, $matrix_file);
|
|
59
|
|
60 my $first_line = <M>;
|
|
61 $first_line =~s/\n//g; $first_line =~s/\r//g; # ne garder que la première ligne du tableau
|
|
62 my @samples = split(/\t/,$first_line); # mettre dans une liste (@samples) chaque intitulé de colonne = N_Id
|
|
63 # print "$first_line\n";
|
|
64 # exit;
|
|
65
|
|
66 # Le but ici est de récupérer les combinaisons associées à chaque espèce : NF, NG et NL
|
|
67 my %hCombination =(); # HASH -> key: esp ; val: combinaison
|
|
68
|
|
69 for (my $i=1; $i <= $#samples; $i++){ # on parcourt chaque colonne ($i) mais on ne regarde que le N_Id
|
|
70 my $header = $samples[$i]; # on récupère le N_Id dans $header (soit le nom de la colonne i)
|
|
71 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
|
|
72 $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
|
|
73 # print "$header\n";
|
|
74 # exit;
|
|
75 }
|
|
76
|
|
77
|
|
78 # foreach my $species (keys (%hCombination)){
|
|
79 # my $combination = $hCombination{$species};
|
|
80 # print "$species $combination\n";
|
|
81 # }
|
|
82
|
|
83
|
|
84 # exit;
|
|
85
|
|
86 # orthogrp présents :
|
|
87 my %hCombination_prs = (); # HASH -> key: combinaison ; val: liste des orthogroupes
|
|
88 # orthogrp absents :
|
|
89 my %hCombination_abs = (); # idem
|
|
90
|
|
91
|
|
92
|
|
93 my %coregenes = (); # HASH -> key: gene ; val: orthogroupe (pour core-genome)
|
|
94 my %specificgenes = (); # HASH -> key: gene ; val: orthogroupe (pour gènes spécifiques)
|
|
95 my %accessorygenes = (); # HASH -> key: gene ; val: orthogroupe (pour génome accessoire)
|
|
96
|
|
97 my $coregene_line;
|
|
98 my %coregenes2 = (); # HASH -> key1: colonne i ; key2: gène ; val: orthogroupe
|
|
99 my %specificgenes2 = (); # HASH -> key1: colonne i ; key2: gène ; val: orthogroupe
|
|
100
|
|
101 my %Genes_of_OG = (); # HASH -> key1: orthogroupe ; key2: colonne i ; val: gène
|
|
102
|
|
103 my %coregenes3 = (); #ligne complete
|
|
104 my %Type_count_byStrain = ();
|
|
105 my %OG_genes = ();
|
|
106 my $nb_genes_total = 0;
|
|
107 my %specificgenes3 = ();
|
|
108 my %Species_Total_Count = (); # HASH -> key: espèce ; val: comptage du nombre total de gènes pour cette espèce
|
|
109 my %Genes_Species_Total = ();
|
|
110 my %NonStrict_Spe = ();
|
|
111
|
|
112 while(<M>) {
|
|
113
|
|
114 my $line = $_;
|
|
115 $line =~s/\n//g; $line =~s/\r//g;
|
|
116 my $nb_found = 0;
|
|
117 my @infos = split(/\t/,$line);
|
|
118 my $orthogroup = $infos[0]; # on récupère le nom de l'orthogroupe dans $orthogroup
|
|
119 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
|
|
120 my $combi_prs = "";
|
|
121 my $combi_abs = "";
|
|
122 my $val;
|
|
123 my $gene_random;
|
|
124 my $unique_col_detected;
|
|
125 my %comptage_especes;
|
|
126 my $seule_espece;
|
|
127
|
|
128 for (my $i=1; $i <= $#infos; $i++){ # on travaille par ligne puis dans chaque ligne (while(<M>)), cellule par cellule (cette boucle for)
|
|
129
|
|
130 $val = $infos[$i]; # on récupère l'information contenue dans la case $i
|
|
131
|
|
132 if ($val =~/\w/){ # s'il cette cellule contient qq chose...
|
|
133 $combi_prs .= "_".$i; # ...on va concaténer notre chaine $combi_prs pour que cela forme une combinaison
|
|
134 $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)
|
|
135 $gene_random=$val; # on récupère la valeur de la case (les gènes)
|
|
136 $unique_col_detected = $i;
|
|
137 my $espece = $hSpecies{$samples[$i]};
|
|
138
|
|
139 my @table_genes = split (',', $val);
|
|
140 my $premier_gene = $table_genes[0];
|
|
141 # $Genes_of_OG{$i}{$orthogroup} = $premier_gene; # pour chaque orthorgoupe de chaque colonne, on récupère le premier gène
|
|
142 foreach my $genes (@table_genes){
|
|
143 # $OG_genes{$orthogroup}{$genes} = 1;
|
|
144 $Genes_of_OG{$i}{$orthogroup} .= $genes;
|
|
145 $nb_genes_total++;
|
|
146 $Species_Total_Count{$espece}++;
|
|
147 $Genes_Species_Total{$genes} = $espece;
|
|
148 $seule_espece = $espece;
|
|
149 $comptage_especes{$espece} .= $genes;
|
|
150
|
|
151 }
|
|
152 }
|
|
153
|
|
154 else { # si jamais il n'y a rien dans la cellule...
|
|
155 $combi_abs .= "_".$i; # ... on fait la même chose mais avec $combi_abs
|
|
156 }
|
|
157
|
|
158 }
|
|
159
|
|
160 if (scalar keys(%comptage_especes) == 1){
|
|
161 my $list = $comptage_especes{$seule_espece};
|
|
162 my @table = split(" ",$list);
|
|
163 foreach my $gene (@table){
|
|
164 $NonStrict_Spe{$gene} = $seule_espece;
|
|
165 }
|
|
166 # print $_;
|
|
167 # print $list."\n";
|
|
168 }
|
|
169
|
|
170 # $hCount{$combi}++;
|
|
171 $hCombination_prs{$combi_prs}.=$orthogroup."\n"; # à la fin de chaque ligne, on va ajouter notre orthogroupe à la combinaison qui lui correspond
|
|
172 $hCombination_abs{$combi_abs}.=$orthogroup."\n";
|
|
173
|
|
174
|
|
175
|
|
176 if ($nb_found == $#infos){ # si nb_found = au nombre de souche, c'est qu'on a à faire à un core-génome
|
|
177 # print "$orthogroup\n";
|
|
178 # print "$nb_found\n=================\n";
|
|
179
|
|
180 for (my $i=1; $i <= $#infos; $i++){
|
|
181 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 ',')
|
|
182 foreach my $gene (@list_of_genes){
|
|
183 $coregenes3{$samples[$i]}{$gene} = 1;
|
|
184 $Type_count_byStrain{"core"}{$samples[$i]}{"oui"}++;
|
|
185 $Type_count_byStrain{"accessory"}{$samples[$i]}{"non"}++;
|
|
186 $Type_count_byStrain{"unique"}{$samples[$i]}{"non"}++;
|
|
187
|
|
188 }
|
|
189
|
|
190 # $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)
|
|
191 # $coregenes2{$i}{$first_gene}= $orthogroup;
|
|
192
|
|
193 }
|
|
194 my @liste_of_genes = split (',', $gene_random);
|
|
195 my $first_gene = $liste_of_genes[0];
|
|
196 $coregenes{$first_gene}= $orthogroup;
|
|
197
|
|
198
|
|
199 if (!$coregene_line){
|
|
200 $coregene_line = $line;
|
|
201 }
|
|
202 }
|
|
203
|
|
204 elsif ($nb_found == 1) { # si on a un gène spé
|
|
205
|
|
206 # # print "$gene_random\n";
|
|
207 # # print "$line\n";
|
|
208 # # print "$unique_col_detected\n";
|
|
209
|
|
210
|
|
211 for (my $i=1; $i <= $#infos; $i++){
|
|
212 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 ',')
|
|
213 foreach my $gene (@list_of_genes){
|
|
214 $specificgenes3{$samples[$i]}{$gene} = 1;
|
|
215 $Type_count_byStrain{"unique"}{$samples[$i]}{"oui"}++;
|
|
216 $Type_count_byStrain{"core"}{$samples[$i]}{"non"}++;
|
|
217 $Type_count_byStrain{"accessory"}{$samples[$i]}{"non"}++;
|
|
218
|
|
219 }
|
|
220
|
|
221
|
|
222 }
|
|
223 # my @list_of_genes = split (',', $gene_random); # idem, on ne veut qu'un seul gène donc on crée la liste
|
|
224 # my $first_gene = $list_of_genes[0]; # on ne prend que le premier
|
|
225 # # print "$first_gene\n";
|
|
226 # # exit;
|
|
227 # $specificgenes{$first_gene}= $orthogroup; # et pareil on crée la table de hash
|
|
228 # $specificgenes2{$unique_col_detected}{$first_gene}= $orthogroup;
|
|
229 }
|
|
230
|
|
231 else { # là c'est le génome accessoire, i.e tout le reste !
|
|
232 # for (my $i=1; $i <= $#infos; $i++){
|
|
233 # 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 ',')
|
|
234 # my $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement !
|
|
235 # $accessorygenes{$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)
|
|
236
|
|
237
|
|
238 # }
|
|
239
|
|
240 for (my $i=1; $i <= $#infos; $i++){
|
|
241 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 ',')
|
|
242 foreach my $gene (@list_of_genes){
|
|
243 # $coregenes3{$samples[$i]}{$gene} = 1;
|
|
244 $Type_count_byStrain{"accessory"}{$samples[$i]}{"oui"}++;
|
|
245 $Type_count_byStrain{"core"}{$samples[$i]}{"non"}++;
|
|
246 $Type_count_byStrain{"unique"}{$samples[$i]}{"non"}++;
|
|
247 }
|
|
248
|
|
249
|
|
250
|
|
251 }
|
|
252
|
|
253
|
|
254 my @liste_of_genes = split (',', $gene_random);
|
|
255 my $first_gene = $liste_of_genes[0];
|
|
256 $accessorygenes{$first_gene}= $orthogroup;
|
|
257 }
|
|
258
|
|
259 }
|
|
260
|
|
261 #print scalar keys(%Genes_of_OG);exit;
|
|
262 # print "$nb_genes_total\n";
|
|
263
|
|
264 # foreach my $og (keys %OG_genes) {
|
|
265 # foreach my $gene (keys %{$OG_genes{$og}}) {
|
|
266 # # print "$og\t$gene\n";
|
|
267 # print $OG_genes{$og}."\n";
|
|
268 # }
|
|
269 # }
|
|
270 # exit;
|
|
271 # foreach my $gene (keys (%Genes_Species_Total)) {
|
|
272 # print "$gene => ".$Genes_Species_Total{$gene}."\n";
|
|
273 # }
|
|
274 # foreach my $strain (keys %specificgenes3) {
|
|
275 # foreach my $gene (keys %{$specificgenes3{$strain}}) {
|
|
276 # print "$strain\t$gene\n";
|
|
277 # }
|
|
278 # }
|
|
279 # foreach my $gene (keys (%NonStrict_Spe)){
|
|
280 # print $NonStrict_Spe{$gene}."\t$gene \n";
|
|
281 # }
|
|
282 # exit;
|
|
283
|
|
284 my %hCol_Annotated = (); # HASH -> key: colonne ; val: 1 (colonnes pour lesquelles les GFF sont présents)
|
|
285
|
|
286 # Le but ici est de ne garder que les colonnes (donc les souches) qui ont un fichier GFF associé
|
|
287 my @list_column = split ('\t', $coregene_line);
|
|
288 for (my $i=1; $i <= $#list_column; $i++){
|
|
289 my @list_genes = split (', ', $list_column[$i]);
|
|
290 my $premier_gene = $list_genes[0];
|
|
291 my $strain = $samples[$i]; # récupérer le nom de la souche
|
|
292
|
|
293
|
|
294 foreach my $gff (@list_gff){
|
|
295 my $result_grep = `grep $premier_gene $gff`;
|
|
296
|
|
297 if ($result_grep){
|
|
298 $hCol_Annotated{$i}=$strain;
|
|
299
|
|
300 }
|
|
301 # print "$result_grep\n";
|
|
302 }
|
|
303 }
|
|
304 # exit;
|
|
305 # foreach my $i (sort keys (%specificgenes2)){ # parcours de la table %hCount2 au niveau des catégories
|
|
306 # foreach my $gene (keys %{$specificgenes2{$i} }){ # parcours de la table %hCount2 au niveau des espèces
|
|
307 # print "$i\t$gene\t".$specificgenes2{$i}{$gene}."\n";
|
|
308 # }
|
|
309 # }
|
|
310 # exit;
|
|
311 # while (my ($k,$v) = each(%accessorygenes)) {
|
|
312 # print "gene=$k OG=$v\n";
|
|
313 # }
|
|
314 # exit;
|
|
315 # foreach my $oups (keys (%coregenes)) {
|
|
316 # print "$oups\n";
|
|
317 # }
|
|
318 # exit;
|
|
319
|
|
320 close (M);
|
|
321
|
|
322 my %Hash_Specific = (); # HASH -> key: orthogroupe ; val: espèce
|
|
323
|
|
324 open (OUT, '>', $output) or die $!;
|
|
325 # print OUT "$annotation\n";
|
|
326 foreach my $species (keys (%hCombination)){ # parcours de la table de hash %hCombination (key: nom esp ; val: combi)
|
|
327 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
|
|
328 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
|
|
329 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
|
|
330
|
|
331 # open (OUT,">results.list.txt");
|
|
332
|
|
333 if ($ortho_presents){
|
|
334 print OUT "> $species - present\n";
|
|
335 print OUT "$ortho_presents\n";
|
|
336 my @orthogroups_name = split ('\n', $ortho_presents);
|
|
337 foreach my $ortho (@orthogroups_name){
|
|
338 $Hash_Specific{$ortho} = $species;
|
|
339 }
|
|
340 }
|
|
341
|
|
342 if ($ortho_absents){
|
|
343 # open (OUT2,">$species.$combination.absents.list.txt");
|
|
344 print OUT "> $species - absent\n";
|
|
345 print OUT "$ortho_absents\n";
|
|
346 }
|
|
347
|
|
348 # close(OUT2);
|
|
349 }
|
|
350
|
|
351 close(OUT);
|
|
352
|
|
353 my %Gene_Specie_Spe = (); # HASH -> key: gène spé ; val: espèce
|
|
354 my %Species_Spe_Count = (); # HASH -> key: espèce ; val: comptage du nombre de gènes spécifiques à cette espèce
|
|
355
|
|
356 foreach my $i (keys(%Genes_of_OG)){
|
|
357 foreach my $ortho (keys %{$Genes_of_OG{$i}}){
|
|
358 my $gene = $Genes_of_OG{$i}{$ortho};
|
|
359
|
|
360 if ($Hash_Specific{$ortho}){
|
|
361 my $specie = $Hash_Specific{$ortho};
|
|
362
|
|
363 my @liste_genes = split(' ',$gene);
|
|
364 foreach my $g(@liste_genes){
|
|
365 $Gene_Specie_Spe{$g} = $specie;
|
|
366 $Species_Spe_Count{$specie}++;
|
|
367
|
|
368 }
|
|
369
|
|
370 }
|
|
371 }
|
|
372 }
|
|
373
|
|
374
|
|
375 # exit;
|
|
376
|
|
377 # my @table_keys = ();
|
|
378 my $nb_groupSpe_genes = 0;
|
|
379
|
|
380 foreach my $gene (keys (%Gene_Specie_Spe)) {
|
|
381 my @table_keys = split (' ', $gene);
|
|
382 foreach my $unique_gene (@table_keys) {
|
|
383 $nb_groupSpe_genes++;
|
|
384 }
|
|
385 }
|
|
386
|
|
387 # print scalar keys (%Gene_Specie_Spe)."\n";
|
|
388 # while (my ($k,$v) = each(%Gene_Specie_Spe)) {
|
|
389 # if ($v =~/ruberi/) {
|
|
390 # print "gene=$k espece=$v\n";
|
|
391 # }
|
|
392 # }
|
|
393 # foreach my $sp (keys (%Species_Spe_Count)){
|
|
394 # print "$sp => ".$Species_Spe_Count{$sp}."\n";
|
|
395 # }
|
|
396 # exit;
|
|
397 #//////////////////////////////////////////////////////////////////////////////////////////////////
|
|
398
|
|
399 ############################################### COG ###############################################
|
|
400
|
|
401 # STEP 1 : CORRESPONDANCE ENTRE LES DIFFERENTS FICHIERS DE COG ET L'ORDRE --------------------------------------------
|
|
402 my @files = split(',', $annotation); # liste des différents fichiers COG (qui se retrouvent dans le dossier Naegleria)
|
|
403 my @list = split(',', $order); # liste de l'ordre des souches
|
|
404 #my ($f,$l);
|
|
405
|
|
406 my %hCorrespondance = (); #HASH -> key: un fichier COG ; val: un nom de souche (ces 2 données sont entrées en input = $annotation et $order)
|
|
407
|
|
408 # ++++++++++++ parcours de 2 listes en même temps ++++++++++++ #
|
|
409 my $l = 1;
|
|
410 foreach my $f (@files){
|
|
411 $hCorrespondance{$f} = $list[$l]; # on fait correspondre pour chaque fichier de COG, un nom de souche
|
|
412 $l++;
|
|
413 }
|
|
414
|
|
415
|
|
416
|
|
417
|
|
418 # #Affichage du hash
|
|
419 # foreach my $f (keys %hCorrespondance){
|
|
420 # print $f."=>".$hCorrespondance{$f}."\n"
|
|
421 # }
|
|
422 # exit;
|
|
423
|
|
424 # STEP 2 : POUR CHAQUE FICHIER DE COG, FAIRE CORRESPONDRE L'ESPECE (ET NON LA SOUCHE) -------------------------------------
|
|
425 my %hCorresp_file_species = (); # HASH -> key: un fichier de COG ; val: une espèce
|
|
426 my %species_names; # HASH -> key: nom d'espèce ; val: 1
|
|
427
|
|
428 foreach my $h (keys (%hCorrespondance)){ # parcours de la table de hash {fichier COG => nom souche}
|
|
429 my $smpl = $hCorrespondance{$h}; # $smpl prend la valeur de la clé (donc d'un nom de souche)
|
|
430 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
|
|
431 $species_names{$espece} = 1; # on garde sous le coude nos nom d'espèce dans cette nouvelle table de hash
|
|
432 $hCorresp_file_species{$h} = $espece; # BUT ATTEINT : on donne pour chaque fichier de COG le nom de l'espèce qui lui correspond
|
|
433 }
|
|
434 # while (my ($k,$v) = each(%hCorresp_file_species)) {
|
|
435 # print "file=$k sp=$v\n";
|
|
436 # }
|
|
437 # exit;
|
|
438
|
|
439
|
|
440
|
|
441
|
|
442 # STEP 3 : COMPTAGE DES CATEGORIES DE COG ------------------------------------------------------------------------------
|
|
443 my %hCount2 = (); # HASH -> key1: catégorie de COG ; key2: espèce associée ; val: comptage
|
|
444
|
|
445 # comptage du core-genome / des gènes spé / du génome accessoire
|
|
446 my %hCore_Count = (); # HASH -> key: catégorie de COG ; val: comptage (ce hash ne sera utilisé que pour le core-genome)
|
|
447 my %hSpecific_Count = (); # HASH -> key: catégorie de COG ; val: comptage
|
|
448 my %hAccessory_Count = (); # HASH -> key: catégorie de COG ; val: comptage
|
|
449
|
|
450 # hash pour récupérer le gène
|
|
451 my %hCore_Cat = (); # HASH -> key: catégorie de COG ; val: gène
|
|
452 my %hAccessory_Cat = (); # HASH -> key: catégorie de COG ; val: gène
|
|
453 my %hSpecific_Cat = (); # HASH -> key: catégorie de COG ; val: gène
|
|
454
|
|
455 # hash pour récupérer le gène
|
|
456 my %hCore_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
|
|
457 my %hAccessory_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
|
|
458 my %hSpecific_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
|
|
459
|
|
460 my %Acc_Cat_Esp_Count = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: comptage
|
|
461
|
|
462 my %Cog_of_gene = (); # HASH -> key: gène ; val: cat de COG
|
|
463 my %Cogs_of_gene = (); # HASH -> key: gène ; val: cat de COG (plusieurs)
|
|
464 my %Specie_of_gene = (); # HASH -> key: gène ; val: souche
|
|
465
|
|
466 my %Global_Count = ();
|
|
467 my %Species_Count = ();
|
|
468 my %Species_NonStrictSpe_Count = ();
|
|
469
|
|
470 my %Genes_in_COG = ();
|
|
471 my %Count_Spe_Genes = ();
|
|
472 my %Count_Total_Species = ();
|
|
473 my %Count_NonStrictSpe_Genes = ();
|
|
474
|
|
475 my %Nveau = ();
|
|
476
|
|
477
|
|
478 foreach my $file(@files){ # parcours de la liste des fichiers
|
|
479 my $esp = $hCorresp_file_species{$file}; # on récupère l'espèce pour chaque fichier de COG dans $esp
|
|
480 # print $esp."\n";
|
|
481 # exit;
|
|
482
|
|
483 my %hCount = (); # HASH -> key: catégorie de COG ; val: comptage
|
|
484
|
|
485
|
|
486 open (A, $file); # on va parcourir maintenant chaque fichier un à un
|
|
487 my $strain = $hCorrespondance{$file};
|
|
488
|
|
489 while (my $line2 = <A>){
|
|
490
|
|
491 $line2 =~s/\n//g; $line2 =~s/\r//g; # on procède ligne par ligne
|
|
492 my @Genes = split('\t', $line2);
|
|
493 my $cogs = $line2;
|
|
494 my $gene = $Genes[0];
|
|
495 my $cog_id = $Genes[1];
|
|
496 $cogs =~s/$gene//g; $cogs =~s/$cog_id//g;
|
|
497 my $first_cat = $Genes[2];
|
|
498 $Cog_of_gene{$gene} = $first_cat;
|
|
499 $Cogs_of_gene{$gene} = $cogs;
|
|
500
|
|
501 $Genes_in_COG{$gene} = $esp;
|
|
502
|
|
503 for (my $j=2; $j <= $#Genes; $j++) {
|
|
504 my $cat = $Genes[$j]; # on récupère la ou les catégorie(s) de COG
|
|
505 $hCount{$cat}++; # pour la catégorie donnée, on incrémente son nb d'occurences
|
|
506
|
|
507
|
|
508 if ($coregenes{$gene}){ # si le $gene fait bien partie du core-genome (donc de notre table de hash %coregenes)
|
|
509 $hCore_Count{$cat}++; # on incrémente le hash
|
|
510 $hCore_Cat{$cat}=$gene; # on récupère le nom du gène
|
|
511 }
|
|
512
|
|
513 if ($accessorygenes{$gene}){ # s'il fait partie des gènes accessoires
|
|
514
|
|
515 $hAccessory_Count{$cat}++;
|
|
516 $hAccessory_Cat{$cat}=$gene;
|
|
517
|
|
518 # if ($accessorygenes{$gene} && $Gene_Specie_Spe{$gene}){
|
|
519 # print "$gene\n";
|
|
520 # # my $espece = $Gene_Specie_Spe{$gene};
|
|
521 # # print "$espece\n";
|
|
522 # # $Nveau{$cat}{$espece}++;
|
|
523 # }
|
|
524
|
|
525 }
|
|
526 if ($coregenes3{$strain}{$gene}){
|
|
527 $Global_Count{"core"}{$cat}{$strain}{"oui"}++;
|
|
528 $Global_Count{"accessory"}{$cat}{$strain}{"non"}++;
|
|
529 $Global_Count{"unique"}{$cat}{$strain}{"non"}++;
|
|
530 }
|
|
531 elsif ($specificgenes3{$strain}{$gene}){
|
|
532 $Global_Count{"unique"}{$cat}{$strain}{"oui"}++;
|
|
533 $Global_Count{"core"}{$cat}{$strain}{"non"}++;
|
|
534 $Global_Count{"accessory"}{$cat}{$strain}{"non"}++;
|
|
535 }
|
|
536 else {
|
|
537 $Global_Count{"accessory"}{$cat}{$strain}{"oui"}++;
|
|
538 $Global_Count{"core"}{$cat}{$strain}{"non"}++;
|
|
539 $Global_Count{"unique"}{$cat}{$strain}{"non"}++;
|
|
540 }
|
|
541
|
|
542
|
|
543
|
|
544
|
|
545
|
|
546
|
|
547 if ($Gene_Specie_Spe{$gene}) {
|
|
548 $Species_Count{$esp}{$cat}{$strain}{"oui"}++;
|
|
549 }
|
|
550 else {
|
|
551 $Species_Count{$esp}{$cat}{$strain}{"non"}++;
|
|
552 }
|
|
553
|
|
554
|
|
555
|
|
556 if ($NonStrict_Spe{$gene}) {
|
|
557 $Species_NonStrictSpe_Count{$esp}{$cat}{$strain}{"oui"}++;
|
|
558 }
|
|
559 else {
|
|
560 $Species_NonStrictSpe_Count{$esp}{$cat}{$strain}{"non"}++;
|
|
561 }
|
|
562
|
|
563 # $Global_Count{$cat}{"accessory"}{$strain}++;
|
|
564
|
|
565
|
|
566 # if ($specificgenes{$gene}){ # s'il fait partie des gènes spécifiques
|
|
567 # $hSpecific_Count{$cat}++;
|
|
568 # $hSpecific_Cat{$cat}=$gene;
|
|
569 # }
|
|
570 # $hCount2{$cat}{$esp}++; # TABLE DE HASH AVEC CLES=CAT DE COG + ESPECE VAL=COMPTAGE
|
|
571
|
|
572
|
|
573 }
|
|
574 }
|
|
575
|
|
576
|
|
577
|
|
578 close (A);
|
|
579
|
|
580 # foreach my $espece (sort keys (%Species_NonStrictSpe_Count)) {
|
|
581 # foreach my $cat (sort keys %{$Species_NonStrictSpe_Count{$espece}}) {
|
|
582 # foreach my $strain (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}}) {
|
|
583 # foreach my $choix (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}{$strain}}) {
|
|
584
|
|
585 # print "$espece - $cat - $strain - $choix ". $Species_NonStrictSpe_Count{$espece}{$cat}{$strain}{$choix}."\n";
|
|
586 # }
|
|
587 # }
|
|
588
|
|
589 # }
|
|
590 # }
|
|
591 # exit;
|
|
592
|
|
593 # while (my ($k,$v) = each(%hCore_Cat)) {
|
|
594 # print "cat=$k gene=$v\n";
|
|
595 # }
|
|
596 # exit;
|
|
597
|
|
598 # print "$file $esp\n=============\n";
|
|
599 while (my ($k,$v) = each(%hCount)) { # parcours de la table de hash de comptage
|
|
600 # print "cat=$k nb=$v\n";
|
|
601 $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
|
|
602 # 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)
|
|
603 }
|
|
604
|
|
605 # Récupérer les gènes du core-génome
|
|
606 while (my ($cat_core,$gene_core) = each(%hCore_Cat)) {
|
|
607 $hCore_Cat_Esp{$cat_core}{$esp}=$gene_core;
|
|
608 }
|
|
609 # Récupérer les gènes du génome-accessoire
|
|
610 while (my ($cat_acc,$gene_acc) = each(%hAccessory_Cat)) {
|
|
611 $hAccessory_Cat_Esp{$cat_acc}{$esp}=$gene_acc;
|
|
612 }
|
|
613 # Récupérer les gènes spécifique
|
|
614 while (my ($cat_spe,$gene_spe) = each(%hSpecific_Cat)) {
|
|
615 $hSpecific_Cat_Esp{$cat_spe}{$esp}=$gene_spe;
|
|
616 }
|
|
617
|
|
618 ###
|
|
619 while (my ($cat,$count) = each(%hAccessory_Count)) {
|
|
620 $Acc_Cat_Esp_Count{$cat}{$esp}=$count;
|
|
621 }
|
|
622 }
|
|
623 # foreach my $type (sort keys (%Global_Count)) {
|
|
624 # foreach my $cat (sort keys %{$Global_Count{$type}}) {
|
|
625 # foreach my $strain (sort keys %{$Global_Count{$type}{$cat}}) {
|
|
626 # foreach my $choix (sort keys %{$Global_Count{$type}{$cat}{$strain}}) {
|
|
627
|
|
628 # print "$type - $cat - $strain - $choix ". $Global_Count{$type}{$cat}{$strain}{$choix}."\n";
|
|
629 # }
|
|
630 # }
|
|
631
|
|
632 # }
|
|
633 # }
|
|
634 # foreach my $espece (sort keys (%Global_Count)) {
|
|
635 # foreach my $cat (sort keys %{$Species_NonStrictSpe_Count{$espece}}) {
|
|
636 # foreach my $strain (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}}) {
|
|
637 # foreach my $choix (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}{$strain}}) {
|
|
638
|
|
639 # print "$espece - $cat - $strain - $choix ". $Species_NonStrictSpe_Count{$espece}{$cat}{$strain}{$choix}."\n";
|
|
640 # }
|
|
641 # }
|
|
642
|
|
643 # }
|
|
644 # }
|
|
645 # exit;
|
|
646
|
|
647 foreach my $gene (keys (%Genes_in_COG)){
|
|
648 my $espece = $Genes_in_COG{$gene};
|
|
649 if ($Gene_Specie_Spe{$gene}) {
|
|
650 $Count_Spe_Genes{$espece}++;
|
|
651 }
|
|
652 if ($Genes_Species_Total{$gene}) {
|
|
653 $Count_Total_Species{$espece}++;
|
|
654 }
|
|
655 if ($NonStrict_Spe{$gene}) {
|
|
656 $Count_NonStrictSpe_Genes{$espece}++;
|
|
657 }
|
|
658 }
|
|
659
|
|
660
|
|
661 # #############################################
|
|
662 # # p / (1-p) p * (1-q) #
|
|
663 # # odds ratio = ----------- = ----------- #
|
|
664 # # q / (1-q) q * (1-p) #
|
|
665 # #############################################
|
|
666 # # où p : proba qu'un E arrive au groupe A
|
|
667 # # où q : proba que ce même E arrive au groupe B
|
|
668
|
|
669
|
|
670 my @orders = ("D","M","N","O","T","U","V","Y","Z","A","B","J","K","L","C","E","F","G","H","I","P","Q","R","S");
|
|
671 ###################
|
|
672 open (OUT7, ">$output7") or die $!;
|
|
673
|
|
674 # my $nb_files = scalar keys @files;
|
|
675
|
|
676
|
|
677 print OUT7 "\t";
|
|
678
|
|
679 foreach my $category(@orders){
|
|
680 # foreach my $category (sort keys (%Acc_Cat_Esp_Count)) {
|
|
681 # my $cat = $category."\t";
|
|
682 print OUT7 $category."\t";
|
|
683 }
|
|
684
|
|
685 print OUT7 "\n";
|
|
686
|
|
687
|
|
688 # foreach my $category (sort keys (%Global_Count)){
|
|
689 foreach my $type (sort keys (%Global_Count)){
|
|
690
|
|
691 print OUT7 "$type\t";
|
|
692 #foreach my $category (sort keys (%{$Global_Count{$type}})){
|
|
693 foreach my $category(@orders){
|
|
694
|
|
695 foreach my $strain (sort keys (%{$Global_Count{$type}{$category}})){
|
|
696
|
|
697 # foreach my $type (sort keys (%{$Global_Count{$category}{$strain}})){
|
|
698 my $nb_type1; my $nb_type2;
|
|
699
|
|
700 if ($Global_Count{$type}{$category}{$strain}{"non"} && $Global_Count{$type}{$category}{$strain}{"oui"}) {
|
|
701 $nb_type1 = $Type_count_byStrain{$type}{$strain}{"non"} - $Global_Count{$type}{$category}{$strain}{"non"};
|
|
702 $nb_type2 = $Type_count_byStrain{$type}{$strain}{"oui"} - $Global_Count{$type}{$category}{$strain}{"oui"};
|
|
703 }
|
|
704 # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$type}{$strain}."\t"."$nb_type\n";
|
|
705 my $ratio1; my $ratio2;
|
|
706 if ($nb_type1 && $nb_type2) {
|
|
707 $ratio1 = $Global_Count{$type}{$category}{$strain}{"non"}/ $nb_type1;
|
|
708 $ratio2 = $Global_Count{$type}{$category}{$strain}{"oui"} / $nb_type2;
|
|
709 }
|
|
710 my $odds_ratio;
|
|
711
|
|
712 if ($ratio1 && $ratio2) {
|
|
713 $odds_ratio = $ratio2 / $ratio1;
|
|
714 }
|
|
715 # print $strain." ".$Type_count_byStrain{$type}{$strain}{"oui"}."\n";
|
|
716 # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$strain}{$type}{"oui"}."\t"."$nb_type2\t".$Global_Count{$category}{$strain}{$type}{"non"}."\t"."$nb_type1\t$odds_ratio\n";
|
|
717 if ($odds_ratio) {
|
|
718 print OUT7 "$odds_ratio;";
|
|
719 }
|
|
720
|
|
721 }
|
|
722 print OUT7 "\t";
|
|
723 }
|
|
724
|
|
725 print OUT7 "\n";
|
|
726 }
|
|
727
|
|
728
|
|
729 print OUT7 "\n";
|
|
730 close (OUT7);
|
|
731
|
|
732
|
|
733 #//////////////////////////////////////////////
|
|
734 open (OUT8, ">$output8") or die $!;
|
|
735
|
|
736
|
|
737
|
|
738 print OUT8 "\t";
|
|
739
|
|
740 #
|
|
741 foreach my $category(@orders){
|
|
742 print OUT8 $category."\t";
|
|
743 }
|
|
744
|
|
745 print OUT8 "\n";
|
|
746
|
|
747
|
|
748 # foreach my $category (sort keys (%Global_Count)){
|
|
749 foreach my $specie (sort keys (%Species_Count)){
|
|
750 # my $nb_genes_nonSpe = $Species_Total_Count{$specie} - $Species_Spe_Count{$specie};
|
|
751 my $nb_genes_nonSpe = $Count_Total_Species{$specie} - $Count_Spe_Genes{$specie};
|
|
752
|
|
753 print OUT8 "$specie\t";
|
|
754 foreach my $category (sort keys (%{$Species_Count{$specie}})){
|
|
755
|
|
756 foreach my $strain (sort keys (%{$Species_Count{$specie}{$category}})){
|
|
757
|
|
758
|
|
759 my $nb_type1; my $nb_type2;
|
|
760
|
|
761 if ($Species_Count{$specie}{$category}{$strain}{"non"} && $Species_Count{$specie}{$category}{$strain}{"oui"}) {
|
|
762 $nb_type1 = $nb_genes_nonSpe - $Species_Count{$specie}{$category}{$strain}{"non"}; # 1-q
|
|
763 $nb_type2 = $Count_Spe_Genes{$specie} - $Species_Count{$specie}{$category}{$strain}{"oui"}; # 1-p
|
|
764 # $nb_type2 = $Species_Spe_Count{$specie} - $Species_Count{$specie}{$category}{$strain}{"oui"}; # 1-p
|
|
765 # print "$nb_genes_nonSpe - ".$Species_Count{$specie}{$category}{$strain}{"non"}. " $nb_type1\n$nb_groupSpe_genes - ".$Species_Count{$specie}{$category}{$strain}{"oui"}. " $nb_type2\n"; exit;
|
|
766 }
|
|
767
|
|
768 # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$type}{$strain}."\t"."$nb_type\n";
|
|
769 my $ratio1; my $ratio2;
|
|
770 if ($nb_type1 && $nb_type2) {
|
|
771 $ratio1 = $Species_Count{$specie}{$category}{$strain}{"non"}/ $nb_type1;
|
|
772 $ratio2 = $Species_Count{$specie}{$category}{$strain}{"oui"} / $nb_type2;
|
|
773 }
|
|
774 my $odds_ratio;
|
|
775
|
|
776 if ($ratio1 && $ratio2) {
|
|
777 $odds_ratio = $ratio2 / $ratio1;
|
|
778 }
|
|
779 # print $strain." ".$Type_count_byStrain{$type}{$strain}{"oui"}."\n";
|
|
780 # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$strain}{$type}{"oui"}."\t"."$nb_type2\t".$Global_Count{$category}{$strain}{$type}{"non"}."\t"."$nb_type1\t$odds_ratio\n";
|
|
781 if ($odds_ratio) {
|
|
782 print OUT8 "$odds_ratio;";
|
|
783 }
|
|
784
|
|
785 }
|
|
786 print OUT8 "\t";
|
|
787 }
|
|
788 print OUT8 "\n";
|
|
789 }
|
|
790
|
|
791 # print OUT9 "\n";
|
|
792
|
|
793 close (OUT8);
|
|
794
|
|
795 #///////////////////////////////////////////////////
|
|
796 # open (OUT9, '>', $output9) or die $!;
|
|
797
|
|
798
|
|
799
|
|
800 # print OUT9 "\t";
|
|
801
|
|
802 # #
|
|
803 # foreach my $category(@orders){
|
|
804 # print OUT9 $category."\t";
|
|
805 # }
|
|
806
|
|
807 # print OUT9 "\n";
|
|
808
|
|
809
|
|
810 # # foreach my $category (sort keys (%Global_Count)){
|
|
811 # foreach my $specie (sort keys (%Species_NonStrictSpe_Count)){
|
|
812 # # my $nb_genes_nonSpe = $Species_Total_Count{$specie} - $Species_Spe_Count{$specie};
|
|
813 # my $nb_genes_nonSpeNS = $Count_Total_Species{$specie} - $Count_NonStrictSpe_Genes{$specie};
|
|
814
|
|
815 # print OUT9 "$specie\t";
|
|
816 # foreach my $category (sort keys (%{$Species_NonStrictSpe_Count{$specie}})){
|
|
817
|
|
818 # foreach my $strain (sort keys (%{$Species_NonStrictSpe_Count{$specie}{$category}})){
|
|
819
|
|
820
|
|
821 # my $nb_type1; my $nb_type2;
|
|
822
|
|
823 # if ($Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"non"} && $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"oui"}) {
|
|
824 # $nb_type1 = $nb_genes_nonSpeNS - $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"non"}; # 1-q
|
|
825 # $nb_type2 = $Count_NonStrictSpe_Genes{$specie} - $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"oui"}; # 1-p
|
|
826 # # $nb_type2 = $Species_Spe_Count{$specie} - $Species_Count{$specie}{$category}{$strain}{"oui"}; # 1-p
|
|
827 # # print "$nb_genes_nonSpe - ".$Species_Count{$specie}{$category}{$strain}{"non"}. " $nb_type1\n$nb_groupSpe_genes - ".$Species_Count{$specie}{$category}{$strain}{"oui"}. " $nb_type2\n"; exit;
|
|
828 # }
|
|
829
|
|
830 # # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$type}{$strain}."\t"."$nb_type\n";
|
|
831 # my $ratio1; my $ratio2;
|
|
832 # if ($nb_type1 && $nb_type2) {
|
|
833 # $ratio1 = $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"non"}/ $nb_type1;
|
|
834 # $ratio2 = $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"oui"} / $nb_type2;
|
|
835 # }
|
|
836 # my $odds_ratio;
|
|
837
|
|
838 # if ($ratio1 && $ratio2) {
|
|
839 # $odds_ratio = $ratio2 / $ratio1;
|
|
840 # }
|
|
841 # # print $strain." ".$Type_count_byStrain{$type}{$strain}{"oui"}."\n";
|
|
842 # # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$strain}{$type}{"oui"}."\t"."$nb_type2\t".$Global_Count{$category}{$strain}{$type}{"non"}."\t"."$nb_type1\t$odds_ratio\n";
|
|
843 # if ($odds_ratio) {
|
|
844 # print OUT9 "$odds_ratio;";
|
|
845 # }
|
|
846
|
|
847 # }
|
|
848 # print OUT9 "\t";
|
|
849 # }
|
|
850 # print OUT9 "\n";
|
|
851 # }
|
|
852
|
|
853 # # print OUT9 "\n";
|
|
854
|
|
855 # close (OUT9);
|
|
856
|
|
857 # exit;
|
|
858 ########################## sortie de pourcentages ##########################
|
|
859 # my $somme_core = 0;
|
|
860 # my $somme_acc = 0;
|
|
861
|
|
862
|
|
863 # foreach my $cat (keys(%hCore_Count)){
|
|
864 # $somme_core = $somme_core + $hCore_Count{$cat};
|
|
865 # }
|
|
866 # foreach my $category (sort keys (%Acc_Cat_Esp_Count)) {
|
|
867
|
|
868 # foreach my $especeee (keys %{$Acc_Cat_Esp_Count{$category}}) {
|
|
869 # $somme_acc = $somme_acc + $Acc_Cat_Esp_Count{$category}{$especeee};
|
|
870
|
|
871 # }
|
|
872
|
|
873 # }
|
|
874
|
|
875
|
|
876 # print "COG categories\tCore-genome\tAccessory genome\n"."\n";
|
|
877 # # foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces...
|
|
878 # # print $e."\t"; #... où on récupère le nom de celles-ci
|
|
879 # # }
|
|
880 # # print "\n";
|
|
881
|
|
882 # foreach my $category (sort keys (%Acc_Cat_Esp_Count)) { # parcours au niveau de la 1ere clé
|
|
883 # my $nb_core = 0;
|
|
884 # my $somme_totale = 0;
|
|
885 # my $number = 0;
|
|
886 # print $category."\t";
|
|
887 # my $c = 0;
|
|
888 # if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome
|
|
889 # $c = $hCore_Count{$category};
|
|
890 # # $hash_core_pc{$c} = 1;
|
|
891 # $somme_totale = $somme_totale + $c;
|
|
892 # $nb_core = ($c/$somme_core)*100;
|
|
893 # # print "$nb_core\t";
|
|
894 # }
|
|
895
|
|
896 # foreach my $especes (sort keys (%species_names)) {
|
|
897 # my $nb_acc = 0;
|
|
898 # my $acc = 0;
|
|
899
|
|
900 # if ($Acc_Cat_Esp_Count{$category}{$especes}) { # si pour une catégorie et une espèce données, on a un nombre : $nbr prend la valeur de ce dernier
|
|
901 # $acc = $Acc_Cat_Esp_Count{$category}{$especes};
|
|
902 # $number = $number + $acc;
|
|
903 # $somme_totale = $somme_totale + $acc;
|
|
904 # $nb_acc = ($acc/$somme_acc)*100;
|
|
905 # # print "$nb_acc\t";
|
|
906 # }
|
|
907
|
|
908 # }
|
|
909
|
|
910
|
|
911 # print "|\t";
|
|
912 # my $pourcentage_core = ($c/$somme_totale)*100;
|
|
913 # print "$pourcentage_core\t";
|
|
914 # my $pourcentage_acc = ($number/$somme_totale)*100;
|
|
915 # print "$pourcentage_acc\n";
|
|
916
|
|
917 # }
|
|
918
|
|
919 ###
|
|
920 # exit;
|
|
921
|
|
922
|
|
923 ########################## sortie de comptage ##########################
|
|
924 # print "COG categories\tCore-genome\tAccessory genome\n"."\t\t";
|
|
925 # foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces...
|
|
926 # print $e."\t"; #... où on récupère le nom de celles-ci
|
|
927 # }
|
|
928 # print "\n";
|
|
929
|
|
930 # foreach my $category (sort keys (%Acc_Cat_Esp_Count)) { # parcours au niveau de la 1ere clé
|
|
931 # print $category."\t\t";
|
|
932 # my $c = 0;
|
|
933 # if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome
|
|
934 # $c = $hCore_Count{$category};
|
|
935 # print "$c\t";
|
|
936 # }
|
|
937 # foreach my $especes (sort keys (%species_names)) {
|
|
938 # if ($Acc_Cat_Esp_Count{$category}{$especes}) { # si pour une catégorie et une espèce données, on a un nombre : $nbr prend la valeur de ce dernier
|
|
939 # print $Acc_Cat_Esp_Count{$category}{$especes}."\t";
|
|
940
|
|
941 # }
|
|
942 # }
|
|
943
|
|
944
|
|
945 # print "\n";
|
|
946 # }
|
|
947
|
|
948 ###
|
|
949
|
|
950 # foreach my $category (sort keys (%hCount2)) { # on parcourt de nouveau les catégories de notre hash à 2 clés
|
|
951 # print OUT2 $category;
|
|
952
|
|
953 # foreach my $especes (sort keys (%species_names)) { # on parcourt également le hash d'espèces
|
|
954
|
|
955 # my $nbr = 0;
|
|
956 # 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
|
|
957 # $nbr = $hCount2{$category}{$especes};
|
|
958 # }
|
|
959 # STEP 4 : AFFICHAGE DANS LE FICHIER DE SORTIE ------------------------------------------------------------------------------
|
|
960 open (OUT4, ">$output4") or die $!;
|
|
961
|
|
962 print OUT4 "COG categories"."\t"."Core-genome"."\t"."Accessory genome"."\t"."Strain Specific genes"."\n";
|
|
963
|
|
964 foreach my $category (sort keys (%hCount2)){ # parcours de la table %hCount2 au niveau des catégories
|
|
965 my $c = 0;
|
|
966 if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome
|
|
967 $c = $hCore_Count{$category};
|
|
968 }
|
|
969 my $acc = 0;
|
|
970 if ($hAccessory_Count{$category}){ # si cette catégorie existe dans le génome accessoire
|
|
971 $acc = $hAccessory_Count{$category};
|
|
972 }
|
|
973 # my $s = 0;
|
|
974 # if ($hSpecific_Count{$category}){ # si cette catégorie existe dans les gènes spécifiques
|
|
975 # $s = $hSpecific_Count{$category};
|
|
976 # }
|
|
977 print OUT4 "$category\t".$c."\t".$acc."\n";#.$s."\n";
|
|
978
|
|
979 foreach my $especeee (keys %{$hCount2{$category} }){ # parcours de la table %hCount2 au niveau des espèces
|
|
980 # print OUT4 "$especeee\t$category\t"; # affichage des esp puis des cat
|
|
981
|
|
982 # if ($hCore_Cat_Esp{$category}{$especeee}) {
|
|
983 # print OUT4 "$hCore_Cat_Esp{$category}{$especeee}\t";
|
|
984 # }
|
|
985 my $c = 0;
|
|
986 if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome
|
|
987 $c = ($hCore_Count{$category}/scalar keys (%coregenes))*100; # calcul du % du comptage
|
|
988 }
|
|
989 # print OUT4 "$c\t"; # affichage du %
|
|
990
|
|
991 # if ($hAccessory_Cat_Esp{$category}{$especeee}) {
|
|
992 # print OUT4 "$hAccessory_Cat_Esp{$category}{$especeee}\t";
|
|
993 # }
|
|
994 my $acc = 0;
|
|
995 if ($hAccessory_Count{$category}){ # si cette catégorie existe dans le génome accessoire
|
|
996 $acc = ($hAccessory_Count{$category}/scalar keys (%accessorygenes))*100; # calcul du % du comptage
|
|
997 }
|
|
998 # print OUT4 "$acc\t"; # affichage du %
|
|
999
|
|
1000 # # if ($hSpecific_Cat_Esp{$category}{$especeee}) {
|
|
1001 # # print OUT4 "$hSpecific_Cat_Esp{$category}{$especeee}\t";
|
|
1002 # # }
|
|
1003 # my $s = 0;
|
|
1004 # if ($hSpecific_Count{$category}){ # si cette catégorie existe dans les gènes spécifiques
|
|
1005 # $s = ($hSpecific_Count{$category}/scalar keys (%specificgenes))*100; # calcul du % du comptage
|
|
1006 # }
|
|
1007 # # print OUT4 "$s\n"; # affichage du %
|
|
1008 }
|
|
1009 }
|
|
1010 close (OUT4);
|
|
1011
|
|
1012 open (OUT3, ">$output3") or die $!;
|
|
1013 foreach my $category (sort keys (%hCount2)) { # parcours au niveau de la 1ere clé
|
|
1014
|
|
1015 foreach my $especeee (keys %{$hCount2{$category} }) { # parcours au niveau de la 2e clé pour la $category donnée
|
|
1016
|
|
1017 print OUT3 "$category\t$especeee\t$hCount2{$category}{$especeee}\n"; # on crée une sortie qui affiche en somme notre hash %hCount2
|
|
1018 }
|
|
1019 }
|
|
1020
|
|
1021 close (OUT3);
|
|
1022
|
|
1023
|
|
1024 open (OUT2, ">$output2") or die $!;
|
|
1025
|
|
1026 print OUT2 "category";
|
|
1027 foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces...
|
|
1028 print OUT2 "\t".$e; #... où on récupère le nom de celles-ci
|
|
1029 }
|
|
1030 print OUT2 "\n";
|
|
1031
|
|
1032 foreach my $category (sort keys (%hCount2)) { # on parcourt de nouveau les catégories de notre hash à 2 clés
|
|
1033 print OUT2 $category;
|
|
1034
|
|
1035 foreach my $especes (sort keys (%species_names)) { # on parcourt également le hash d'espèces
|
|
1036
|
|
1037 my $nbr = 0;
|
|
1038 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
|
|
1039 $nbr = $hCount2{$category}{$especes};
|
|
1040 }
|
|
1041 # $nbr =~s/\n//g; $nbr =~s/\r//g;
|
|
1042
|
|
1043
|
|
1044 my @liste = split(',', $nbr); # vu qu'il peut y avoir plusieurs nombres on les dissocie
|
|
1045
|
|
1046 my $somme=0;
|
|
1047 my $n=0;
|
|
1048 my $moyenne=0;
|
|
1049 #print "\nma liste de $nbr: ".join("%",@liste)."\n";
|
|
1050 foreach my $x (@liste) { # on parcourt nos nombres
|
|
1051 $somme=$somme+$x;
|
|
1052 $n=$n+1;
|
|
1053 }
|
|
1054
|
|
1055 if ($n>0){
|
|
1056 $moyenne = $somme/$n; # on fait le calcul de la moyenne
|
|
1057 }
|
|
1058 # print "$category, $especes: $hCount2{$category}{$especes}\t";
|
|
1059 # print "moyenne = $moyenne\n=============\n";
|
|
1060
|
|
1061 print OUT2 "\t".$moyenne; # fichier de sortie
|
|
1062 }
|
|
1063 print OUT2 "\n";
|
|
1064 }
|
|
1065
|
|
1066 close (OUT2);
|
|
1067
|
|
1068 # foreach my $cat (keys (%hCore_Cat)){
|
|
1069 # print OUT4 $c_gene."\t";
|
|
1070 # }
|
|
1071
|
|
1072
|
|
1073 #//////////////////////////////////////////////////////////////////////////////////////////////////
|
|
1074
|
|
1075 ############################################### GFF ###############################################
|
|
1076
|
|
1077
|
|
1078 # my @order_gff = split(',', $order_GFF); # liste de l'ordre des souches
|
|
1079 my ($g,$o);
|
|
1080
|
|
1081 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)
|
|
1082 my %Gene_position = ();
|
|
1083 my %Cat_genes = ();
|
|
1084 my %Cat_genes2 = ();
|
|
1085
|
|
1086 my %hash_of_genes = ();
|
|
1087
|
|
1088
|
|
1089 foreach $g (@list_gff){
|
|
1090 # print "$g\n";
|
|
1091 # $hgff_order{$g} = $order_gff[$o++]; # on fait correspondre pour chaque fichier GFF, un nom de souche
|
|
1092 open (G, $g);
|
|
1093 while (<G>) {
|
|
1094 my @table_gff = split (/\t/, $_);
|
|
1095 my $chr = $table_gff[0];
|
|
1096 my $start = $table_gff[3];
|
|
1097 my $end = $table_gff[4];
|
|
1098 my $gene_name = $table_gff[8];
|
|
1099 my $type = $table_gff[2];
|
|
1100
|
|
1101
|
|
1102
|
|
1103 if ($type && $type eq "mRNA" && $gene_name =~ /ID=([^;]+);/){ #or $type eq "CDS"
|
|
1104 my $gene = $1;
|
|
1105 # print $gene."\n";
|
|
1106 # exit;
|
|
1107 $hash_of_genes{$gene}=1;
|
|
1108
|
|
1109 foreach my $cog (keys (%hCore_Cat)){
|
|
1110 if ($hCore_Cat{$cog} eq $gene){
|
|
1111 $Cat_genes{$gene}=$cog;
|
|
1112 }
|
|
1113 }
|
|
1114 foreach my $cog_bis (keys (%hSpecific_Cat)){
|
|
1115 if ($hSpecific_Cat{$cog_bis} eq $gene){
|
|
1116 $Cat_genes2{$gene}=$cog_bis;
|
|
1117 }
|
|
1118 }
|
|
1119
|
|
1120 $Gene_position{$gene}="$chr\t$start\t$end";
|
|
1121 }
|
|
1122
|
|
1123 # foreach my $gene (keys (%hash_of_genes)){
|
|
1124 # my $orthogrp = $hGene_OG{$gene};
|
|
1125 # print "$orthogrp\n";
|
|
1126 # }
|
|
1127 }
|
|
1128
|
|
1129 close (G);
|
|
1130 }
|
|
1131
|
|
1132 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);
|
|
1133
|
|
1134 mkdir("Core");
|
|
1135 foreach my $i (keys (%coregenes2)){
|
|
1136
|
|
1137 if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
|
|
1138 next;
|
|
1139 }
|
|
1140
|
|
1141
|
|
1142 my $strain_name = $hCol_Annotated{$i};
|
|
1143
|
|
1144 my $specie_name = $hSpecies{$strain_name};
|
|
1145
|
|
1146
|
|
1147
|
|
1148 open (OUT5, "> Core/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
|
|
1149 print OUT5 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n";
|
|
1150
|
|
1151 my $refcoregenes2 = $coregenes2{$i};
|
|
1152 my %subhash = %$refcoregenes2;
|
|
1153 foreach my $gene (keys (%subhash)){
|
|
1154 # print "$gene\n";
|
|
1155 my $cat = "unknown";
|
|
1156 if ($Cog_of_gene{$gene}){
|
|
1157 $cat = $Cog_of_gene{$gene};
|
|
1158 }
|
|
1159 # if (!$Gene_position{$gene}){
|
|
1160 # print "$gene\n coucou"; exit;
|
|
1161 # }
|
|
1162
|
|
1163 # if (!$subhash{$gene}){
|
|
1164 # print "$gene\n";
|
|
1165 # }
|
|
1166 print OUT5 $subhash{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n";
|
|
1167
|
|
1168 }
|
|
1169
|
|
1170 close (OUT5);
|
|
1171 }
|
|
1172
|
|
1173 mkdir("StrainSpecific");
|
|
1174 foreach my $i (keys (%specificgenes2)){
|
|
1175
|
|
1176 if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
|
|
1177 next;
|
|
1178 }
|
|
1179
|
|
1180
|
|
1181 my $strain_name = $hCol_Annotated{$i};
|
|
1182
|
|
1183 my $specie_name = $hSpecies{$strain_name};
|
|
1184
|
|
1185
|
|
1186
|
|
1187 open (OUT7, "> StrainSpecific/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
|
|
1188 print OUT7 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n";
|
|
1189
|
|
1190 my $refspecificgenes2 = $specificgenes2{$i};
|
|
1191 my %subhash = %$refspecificgenes2;
|
|
1192 foreach my $gene (keys (%subhash)){
|
|
1193 # print "$gene\n"; exit;
|
|
1194 my $cat = "unknown";
|
|
1195 if ($Cog_of_gene{$gene}){
|
|
1196 $cat = $Cog_of_gene{$gene};
|
|
1197 }
|
|
1198 # if (!$Gene_position{$gene}){
|
|
1199 # print "$gene\n coucou"; exit;
|
|
1200 # }
|
|
1201
|
|
1202 # if (!$subhash{$gene}){
|
|
1203 # print "$gene\n";
|
|
1204 # }
|
|
1205 print OUT7 $subhash{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n";
|
|
1206
|
|
1207 }
|
|
1208
|
|
1209 close (OUT7);
|
|
1210 }
|
|
1211
|
|
1212
|
|
1213 mkdir("GroupSpecific");
|
|
1214 foreach my $i (keys (%Genes_of_OG)){
|
|
1215 if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
|
|
1216 next;
|
|
1217 }
|
|
1218
|
|
1219 my $strain_name = $hCol_Annotated{$i};
|
|
1220
|
|
1221 my $specie_name = $hSpecies{$strain_name};
|
|
1222
|
|
1223 open (OUT6, "> GroupSpecific/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
|
|
1224 print OUT6 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n";
|
|
1225
|
|
1226 my $refGenes_of_OG = $Genes_of_OG{$i};
|
|
1227 my %subhash = %$refGenes_of_OG;
|
|
1228
|
|
1229 foreach my $orthogroup (keys (%subhash)){
|
|
1230 if ($Hash_Specific{$orthogroup} && $Hash_Specific{$orthogroup} eq $specie_name){
|
|
1231 my $gene = $subhash{$orthogroup};
|
|
1232
|
|
1233 my $cat = "unknown";
|
|
1234 if ($Cog_of_gene{$gene}){
|
|
1235 $cat = $Cog_of_gene{$gene};
|
|
1236 }
|
|
1237 print OUT6 $orthogroup."\t".$subhash{$orthogroup}."\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n";
|
|
1238 }
|
|
1239
|
|
1240 }
|
|
1241 close (OUT6);
|
|
1242 }
|