Mercurial > repos > mgarnier > pangenome_cog_analysis

comparison pangenomeCogAnalysis_V1.pl @ 2:0428ce25da81 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mgarnier
date Fri, 02 Jul 2021 14:53:33 +0000
parents 1f75641c2ee8
children db4e1e6850b0
comparison
equal deleted inserted replaced
1:1f75641c2ee8 2:0428ce25da81
2 2
3 use strict; 3 use strict;
4 use warnings; 4 use warnings;
5 5
6 my $num_args = $#ARGV + 1; 6 my $num_args = $#ARGV + 1;
7 if ($num_args != 11) { 7 if ($num_args != 10) {
8 print "Il n'y a pas le bon nombre d'arguments !\n"; 8 print "Il n'y a pas le bon nombre d'arguments !\n";
9 exit; 9 exit;
10 } 10 }
11 11
12 # INPUT_ 12 # INPUT_
20 # OUTPUT_ 20 # OUTPUT_
21 my $output = $ARGV[6]; # liste des espèces avec leurs orthogroupes (présence-absence) 21 my $output = $ARGV[6]; # liste des espèces avec leurs orthogroupes (présence-absence)
22 my $output2 = $ARGV[7]; # fichier des moyennes 22 my $output2 = $ARGV[7]; # fichier des moyennes
23 my $output3 = $ARGV[8]; # fichier de la liste des valeurs pour chaque catégorie de COG et pour chaque espèce 23 my $output3 = $ARGV[8]; # fichier de la liste des valeurs pour chaque catégorie de COG et pour chaque espèce
24 my $output4 = $ARGV[9]; # fichier avec les catégories de COG pour core-génome / génome accessoire / gènes spé 24 my $output4 = $ARGV[9]; # fichier avec les catégories de COG pour core-génome / génome accessoire / gènes spé
25 my $output5 = $ARGV[10]; # sortie qui affiche les GFF
26 25
27 26
28 # print "ok\n"; 27 # print "ok\n";
29 # exit; 28 # exit;
30 29
30 my @list_gff = split(',', $annotation_GFF); # liste des différents fichiers GFF (qui se retrouvent dans le dossier Annotation Maker)
31 my %hSpecies = (); # HASH -> key: N_Id (ex NF_AR12) ; val: nom de l'esp (ex Naegleria Fowleri) 31 my %hSpecies = (); # HASH -> key: N_Id (ex NF_AR12) ; val: nom de l'esp (ex Naegleria Fowleri)
32 32
33 ######################## LE SPECIES_FILE ########################### 33 ######################## LE SPECIES_FILE ###########################
34 open (S, $species_file); 34 open (S, $species_file);
35 while (my $line = <S>){ 35 while (my $line = <S>){
36 36
37 $line =~s/\n//g; $line =~s/\r//g; 37 $line =~s/\n//g; $line =~s/\r//g;
38 my @sp = split('\t', $line); 38 my @sp = split('\t', $line);
39 # print "$line\n"; 39 # print "$line\n";
40 # exit; 40 # exit;
41 $hSpecies{$sp[0]} = $sp[1]; # key = N_Id ; val = name 41 $hSpecies{$sp[0]} = $sp[1]; # HASH -> key: N_Id ; val: name
42 42
43 } 43 }
44 my $nbr = keys (%hSpecies); #compter le nombre de souches max 44 my $nbr = keys (%hSpecies); #compter le nombre de souches max
45 # = taille de la table de hash 45 # = taille de la table de hash
46 # print "J'ai $nbr clés\n"; 46 # print "J'ai $nbr clés\n";
89 89
90 my %coregenes = (); # HASH -> key: gene ; val: orthogroupe (pour core-genome) 90 my %coregenes = (); # HASH -> key: gene ; val: orthogroupe (pour core-genome)
91 my %specificgenes = (); # HASH -> key: gene ; val: orthogroupe (pour gènes spécifiques) 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) 92 my %accessorygenes = (); # HASH -> key: gene ; val: orthogroupe (pour génome accessoire)
93 93
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
94 99
95 while(<M>) { 100 while(<M>) {
96 101
102 my $line = $_;
103 $line =~s/\n//g; $line =~s/\r//g;
97 my $nb_found = 0; 104 my $nb_found = 0;
98 my @infos = split(/\t/,$_); 105 my @infos = split(/\t/,$line);
99 my $orthogroup = $infos[0]; # on récupère le nom de l'orthogroupe dans $orthogroup 106 my $orthogroup = $infos[0]; # on récupère le nom de l'orthogroupe dans $orthogroup
100 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 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
101 my $combi_prs = ""; 108 my $combi_prs = "";
102 my $combi_abs = ""; 109 my $combi_abs = "";
103 my $val; 110 my $val;
104 my $gene_random; 111 my $gene_random;
105 112
113
106 for (my $i=1; $i <= $#infos; $i++){ # on travaille par ligne puis dans chaque ligne (while(<M>)), cellule par cellule (cette boucle for) 114 for (my $i=1; $i <= $#infos; $i++){ # on travaille par ligne puis dans chaque ligne (while(<M>)), cellule par cellule (cette boucle for)
107 115
108 $val = $infos[$i]; # on récupère l'information contenue dans la case $i 116 $val = $infos[$i]; # on récupère l'information contenue dans la case $i
109 117
110 if ($val =~/\w/){ # s'il cette cellule contient qq chose... 118 if ($val =~/\w/){ # s'il cette cellule contient qq chose...
111 $combi_prs .= "_".$i; # ...on va concaténer notre chaine $combi_prs pour que cela forme une combinaison 119 $combi_prs .= "_".$i; # ...on va concaténer notre chaine $combi_prs pour que cela forme une combinaison
112 $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) 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)
113 $gene_random=$val; # on récupère la valeur de la case (les gènes) 121 $gene_random=$val; # on récupère la valeur de la case (les gènes)
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
114 } 126 }
115 127
116 else { # si jamais il n'y a rien dans la cellule... 128 else { # si jamais il n'y a rien dans la cellule...
117 $combi_abs .= "_".$i; # ... on fait la même chose mais avec $combi_abs 129 $combi_abs .= "_".$i; # ... on fait la même chose mais avec $combi_abs
118 } 130 }
122 # $hCount{$combi}++; 134 # $hCount{$combi}++;
123 $hCombination_prs{$combi_prs}.=$orthogroup."\n"; # à la fin de chaque ligne, on va ajouter notre orthogroupe à la combinaison qui lui correspond 135 $hCombination_prs{$combi_prs}.=$orthogroup."\n"; # à la fin de chaque ligne, on va ajouter notre orthogroupe à la combinaison qui lui correspond
124 $hCombination_abs{$combi_abs}.=$orthogroup."\n"; 136 $hCombination_abs{$combi_abs}.=$orthogroup."\n";
125 137
126 138
127
128 139
129 if ($nb_found == $#infos){ # si nb_found = au nombre de souche, c'est qu'on a à faire à un core-génome 140 if ($nb_found == $#infos){ # si nb_found = au nombre de souche, c'est qu'on a à faire à un core-génome
130 # print "$orthogroup\n"; 141 # print "$orthogroup\n";
131 # print "$nb_found\n=================\n"; 142 # print "$nb_found\n=================\n";
132 my @list_of_genes = split (',', $first_column); # ici va séparer tous les gènes (qui se présentent comme une liste, séparés par des ',') 143 for (my $i=1; $i <= $#infos; $i++){
133 my $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement ! 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 ',')
134 $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) 145 my $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement !
135 # foreach my $oups (@list_of_genes) { 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)
136 # print "$orthogroup\n"; 147 $coregenes2{$i}{$first_gene}= $orthogroup;
137 # print "$first_gene\n--------------------------\n"; 148
138 # } 149
139 150
151
152 }
153 if (!$coregene_line){
154 $coregene_line = $line;
155 }
140 } 156 }
141 elsif ($nb_found == 1) { # si on a un gène spé 157 elsif ($nb_found == 1) { # si on a un gène spé
142 my @list_of_genes = split (',', $gene_random); # idem, on ne veut qu'un seul gène donc on crée la liste 158 my @list_of_genes = split (',', $gene_random); # idem, on ne veut qu'un seul gène donc on crée la liste
143 my $first_gene = $list_of_genes[0]; # on ne prend que le premier 159 my $first_gene = $list_of_genes[0]; # on ne prend que le premier
144 $specificgenes{$first_gene}= $orthogroup; # et pareil on crée la table de hash 160 $specificgenes{$first_gene}= $orthogroup; # et pareil on crée la table de hash
149 my $first_gene = $list_of_genes[0]; 165 my $first_gene = $list_of_genes[0];
150 $accessorygenes{$first_gene}= $orthogroup; 166 $accessorygenes{$first_gene}= $orthogroup;
151 } 167 }
152 168
153 } 169 }
154 # while (my ($k,$v) = each(%coregenes)) { 170
155 # print "gene=$k OG=$v\n"; 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";
156 # } 200 # }
157 # exit; 201 # exit;
158 # foreach my $oups (keys (%coregenes)) { 202 # foreach my $oups (keys (%coregenes)) {
159 # print "$oups\n"; 203 # print "$oups\n";
160 # } 204 # }
161 # exit; 205 # exit;
162 206
163 close (M); 207 close (M);
208
209 my %Hash_Specific = ();
164 210
165 open (OUT, '>', $output) or die $!; 211 open (OUT, '>', $output) or die $!;
166 print OUT "$annotation\n"; 212 print OUT "$annotation\n";
167 foreach my $species (keys (%hCombination)){ # parcours de la table de hash %hCombination (key: nom esp ; val: combi) 213 foreach my $species (keys (%hCombination)){ # parcours de la table de hash %hCombination (key: nom esp ; val: combi)
168 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 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
172 # open (OUT,">results.list.txt"); 218 # open (OUT,">results.list.txt");
173 219
174 if ($ortho_presents){ 220 if ($ortho_presents){
175 print OUT "> $species - present\n"; 221 print OUT "> $species - present\n";
176 print OUT "$ortho_presents\n"; 222 print OUT "$ortho_presents\n";
223 my @orthogroups_name = split ('\n', $ortho_presents);
224 foreach my $ortho (@orthogroups_name){
225 $Hash_Specific{$ortho} = $species;
226 }
177 } 227 }
178 228
179 if ($ortho_absents){ 229 if ($ortho_absents){
180 # open (OUT2,">$species.$combination.absents.list.txt"); 230 # open (OUT2,">$species.$combination.absents.list.txt");
181 print OUT "> $species - absent\n"; 231 print OUT "> $species - absent\n";
239 # hash pour récupérer le gène 289 # hash pour récupérer le gène
240 my %hCore_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène 290 my %hCore_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
241 my %hAccessory_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène 291 my %hAccessory_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
242 my %hSpecific_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène 292 my %hSpecific_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène
243 293
244 my %Cog_of_gene = (); 294 my %Cog_of_gene = (); # HASH -> key: gène ; val: cat de COG
295 my %Specie_of_gene = (); # HASH -> key: gène ; val: souche
245 296
246 foreach my $file(@files){ # parcours de la liste des fichiers 297 foreach my $file(@files){ # parcours de la liste des fichiers
247 my $esp = $hCorresp_file_species{$file}; # on récupère l'espèce pour chaque fichier de COG dans $esp 298 my $esp = $hCorresp_file_species{$file}; # on récupère l'espèce pour chaque fichier de COG dans $esp
248 # print $esp."\n"; 299 # print $esp."\n";
249 # exit; 300 # exit;
415 466
416 #////////////////////////////////////////////////////////////////////////////////////////////////// 467 #//////////////////////////////////////////////////////////////////////////////////////////////////
417 468
418 ############################################### GFF ############################################### 469 ############################################### GFF ###############################################
419 470
420 my @list_gff = split(',', $annotation_GFF); # liste des différents fichiers GFF (qui se retrouvent dans le dossier Annotation Maker) 471
421 my @order_gff = split(',', $order_GFF); # liste de l'ordre des souches 472 my @order_gff = split(',', $order_GFF); # liste de l'ordre des souches
422 my ($g,$o); 473 my ($g,$o);
423 474
424 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) 475 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)
425 my %Gene_position = (); 476 my %Gene_position = ();
426 my %Cat_genes = (); 477 my %Cat_genes = ();
427 478
428 # my %hCore_Cat = reverse (%hCore_Cat); 479 my %hash_of_genes = ();
429 480
430 # ++++++++++++ parcours de 2 listes en même temps ++++++++++++ # 481 # ++++++++++++ parcours de 2 listes en même temps ++++++++++++ #
431 foreach $g (@list_gff){ 482 foreach $g (@list_gff){
483 # print "$g\n";
432 $hgff_order{$g} = $order_gff[$o++]; # on fait correspondre pour chaque fichier GFF, un nom de souche 484 $hgff_order{$g} = $order_gff[$o++]; # on fait correspondre pour chaque fichier GFF, un nom de souche
433 open (G, $g); 485 open (G, $g);
434 while (<G>) { 486 while (<G>) {
435 my @table_gff = split (/\t/, $_); 487 my @table_gff = split (/\t/, $_);
436 my $chr = $table_gff[0]; 488 my $chr = $table_gff[0];
443 495
444 if ($type && $type eq "mRNA" && $gene_name =~ /ID=([^;]+);/){ 496 if ($type && $type eq "mRNA" && $gene_name =~ /ID=([^;]+);/){
445 my $gene = $1; 497 my $gene = $1;
446 # print $gene."\n"; 498 # print $gene."\n";
447 # exit; 499 # exit;
500 $hash_of_genes{$gene}=1;
501
448 foreach my $cog (keys (%hCore_Cat)){ 502 foreach my $cog (keys (%hCore_Cat)){
449 if ($hCore_Cat{$cog} eq $gene){ 503 if ($hCore_Cat{$cog} eq $gene){
450 $Cat_genes{$gene}=$cog; 504 $Cat_genes{$gene}=$cog;
451 } 505 }
452 } 506 }
453 507
454 $Gene_position{$gene}="$chr\t$start\t$end"; 508 $Gene_position{$gene}="$chr\t$start\t$end";
455 } 509 }
456 510
457 511 # foreach my $gene (keys (%hash_of_genes)){
512 # my $orthogrp = $hGene_OG{$gene};
513 # print "$orthogrp\n";
514 # }
458 } 515 }
459 516
460 close (G); 517 close (G);
461 } 518 }
462 519
463 open (OUT5, "> $output5"); 520 mkdir("Core");
464 print OUT5 "Orthogroups\tGenes\tChromosomes\tStart\tEnd\tCOG categories\n"; 521 foreach my $i (keys (%coregenes2)){
465 522
466 523 if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
467 524 next;
468 525 }
469 # foreach my $gene1 (keys (%Gene_position)){ 526
470 # if ($hCore_Cat{$gene1}){ 527
471 # # print "$hCore_Cat{$gene1}\n"; 528 my $strain_name = $hCol_Annotated{$i};
472 # my $cog = $hCore_Cat{$gene1}; 529
473 # $Hash_genes{$gene1}=$cog; 530 my $specie_name = $hSpecies{$strain_name};
474 # } 531
475 # } 532
476 # while (my ($k,$v) = each(%Hash_genes)) { 533
477 # print "gene=$k cog=$v\n"; 534 open (OUT5, "> Core/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
478 # } 535 print OUT5 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\n";
479 # exit; 536
480 537 my $refcoregenes2 = $coregenes2{$i};
481 foreach my $gene (keys (%coregenes)){ 538 my %subhash = %$refcoregenes2;
482 # print "$gene\n"; 539 foreach my $gene (keys (%subhash)){
483 my $cat = "unknown"; 540 # print "$gene\n";
484 if ($Cog_of_gene{$gene}){ 541 my $cat = "unknown";
485 $cat = $Cog_of_gene{$gene}; 542 if ($Cog_of_gene{$gene}){
486 } 543 $cat = $Cog_of_gene{$gene};
487 print OUT5 $coregenes{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\n"; 544 }
488 545 # if (!$Gene_position{$gene}){
489 } 546 # print "$gene\n coucou"; exit;
490 547 # }
491 close (OUT5); 548
549 # if (!$subhash{$gene}){
550 # print "$gene\n";
551 # }
552 print OUT5 $subhash{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\n";
553
554 }
555
556 close (OUT5);
557 }
558
559
560 mkdir("GroupSpecific");
561 foreach my $i (keys (%Genes_of_OG)){
562 if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
563 next;
564 }
565
566 my $strain_name = $hCol_Annotated{$i};
567
568 my $specie_name = $hSpecies{$strain_name};
569
570 open (OUT6, "> GroupSpecific/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
571 print OUT6 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\n";
572
573 my $refGenes_of_OG = $Genes_of_OG{$i};
574 my %subhash = %$refGenes_of_OG;
575
576 foreach my $orthogroup (keys (%subhash)){
577 if ($Hash_Specific{$orthogroup} && $Hash_Specific{$orthogroup} eq $specie_name){
578 my $gene = $subhash{$orthogroup};
579
580 my $cat = "unknown";
581 if ($Cog_of_gene{$gene}){
582 $cat = $Cog_of_gene{$gene};
583 }
584 print OUT6 $orthogroup."\t".$subhash{$orthogroup}."\t".$Gene_position{$gene}."\t".$cat."\n";
585 }
586
587 }
588 close (OUT6);
589 }