Mercurial > repos > mgarnier > pangenome_cog_analysis
comparison pangenomeCogAnalysis_V1.pl @ 0:731fb6cb324b draft
Uploaded
| author | mgarnier |
|---|---|
| date | Wed, 30 Jun 2021 13:30:19 +0000 |
| parents | |
| children | 1f75641c2ee8 |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:731fb6cb324b |
|---|---|
| 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[6]; # liste des espèces avec leurs orthogroupes (présence-absence) | |
| 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 | |
| 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 | |
| 27 | |
| 28 # print "ok\n"; | |
| 29 # exit; | |
| 30 | |
| 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; | |
| 41 $hSpecies{$sp[0]} = $sp[1]; # key = N_Id ; val = name | |
| 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 | |
| 94 | |
| 95 while(<M>) { | |
| 96 | |
| 97 my $nb_found = 0; | |
| 98 my @infos = split(/\t/,$_); | |
| 99 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 | |
| 101 my $combi_prs = ""; | |
| 102 my $combi_abs = ""; | |
| 103 my $val; | |
| 104 my $gene_random; | |
| 105 | |
| 106 for (my $i=1; $i <= $#infos; $i++){ # on travaille par ligne puis dans chaque ligne (while(<M>)), cellule par cellule (cette boucle for) | |
| 107 | |
| 108 $val = $infos[$i]; # on récupère l'information contenue dans la case $i | |
| 109 | |
| 110 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 | |
| 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) | |
| 113 $gene_random=$val; # on récupère la valeur de la case (les gènes) | |
| 114 } | |
| 115 | |
| 116 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 | |
| 118 } | |
| 119 | |
| 120 } | |
| 121 | |
| 122 # $hCount{$combi}++; | |
| 123 $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"; | |
| 125 | |
| 126 | |
| 127 | |
| 128 | |
| 129 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"; | |
| 131 # 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 ',') | |
| 133 my $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement ! | |
| 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) | |
| 135 # foreach my $oups (@list_of_genes) { | |
| 136 # print "$orthogroup\n"; | |
| 137 # print "$first_gene\n--------------------------\n"; | |
| 138 # } | |
| 139 | |
| 140 } | |
| 141 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 | |
| 143 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 | |
| 145 } | |
| 146 | |
| 147 else { # là c'est le génome accessoire, i.e tout le reste ! | |
| 148 my @list_of_genes = split (',', $gene_random); | |
| 149 my $first_gene = $list_of_genes[0]; | |
| 150 $accessorygenes{$first_gene}= $orthogroup; | |
| 151 } | |
| 152 | |
| 153 } | |
| 154 # while (my ($k,$v) = each(%coregenes)) { | |
| 155 # print "gene=$k OG=$v\n"; | |
| 156 # } | |
| 157 # exit; | |
| 158 # foreach my $oups (keys (%coregenes)) { | |
| 159 # print "$oups\n"; | |
| 160 # } | |
| 161 # exit; | |
| 162 | |
| 163 close (M); | |
| 164 | |
| 165 open (OUT, '>', $output) or die $!; | |
| 166 print OUT "$annotation\n"; | |
| 167 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 | |
| 169 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 | |
| 170 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 | |
| 171 | |
| 172 # open (OUT,">results.list.txt"); | |
| 173 | |
| 174 if ($ortho_presents){ | |
| 175 print OUT "> $species - present\n"; | |
| 176 print OUT "$ortho_presents\n"; | |
| 177 } | |
| 178 | |
| 179 if ($ortho_absents){ | |
| 180 # open (OUT2,">$species.$combination.absents.list.txt"); | |
| 181 print OUT "> $species - absent\n"; | |
| 182 print OUT "$ortho_absents\n"; | |
| 183 } | |
| 184 | |
| 185 # close(OUT2); | |
| 186 } | |
| 187 | |
| 188 close(OUT); | |
| 189 | |
| 190 #////////////////////////////////////////////////////////////////////////////////////////////////// | |
| 191 | |
| 192 ############################################### COG ############################################### | |
| 193 | |
| 194 # STEP 1 : CORRESPONDANCE ENTRE LES DIFFERENTS FICHIERS DE COG ET L'ORDRE -------------------------------------------- | |
| 195 my @files = split(',', $annotation); # liste des différents fichiers COG (qui se retrouvent dans le dossier Naegleria) | |
| 196 my @list = split(',', $order); # liste de l'ordre des souches | |
| 197 my ($f,$l); | |
| 198 | |
| 199 my %hCorrespondance = (); #HASH -> key: un fichier COG ; val: un nom de souche (ces 2 données sont entrées en input = $annotation et $order) | |
| 200 | |
| 201 # ++++++++++++ parcours de 2 listes en même temps ++++++++++++ # | |
| 202 foreach $f (@files){ | |
| 203 $hCorrespondance{$f} = $list[$l++]; # on fait correspondre pour chaque fichier de COG, un nom de souche | |
| 204 } | |
| 205 | |
| 206 # #Affichage du hash | |
| 207 # foreach $f (keys %hCorrespondance){ | |
| 208 # print $f."=>".$hCorrespondance{$f}."\n" | |
| 209 # } | |
| 210 | |
| 211 # STEP 2 : POUR CHAQUE FICHIER DE COG, FAIRE CORRESPONDRE L'ESPECE (ET NON LA SOUCHE) ------------------------------------- | |
| 212 my %hCorresp_file_species = (); # HASH -> key: un fichier de COG ; val: une espèce | |
| 213 my %species_names; # HASH -> key: nom d'espèce ; val: 1 | |
| 214 | |
| 215 foreach my $h (keys (%hCorrespondance)){ # parcours de la table de hash {fichier COG => nom souche} | |
| 216 my $smpl = $hCorrespondance{$h}; # $smpl prend la valeur de la clé (donc d'un nom de souche) | |
| 217 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 | |
| 218 $species_names{$espece} = 1; # on garde sous le coude nos nom d'espèce dans cette nouvelle table de hash | |
| 219 $hCorresp_file_species{$h} = $espece; # BUT ATTEINT : on donne pour chaque fichier de COG le nom de l'espèce qui lui correspond | |
| 220 } | |
| 221 # while (my ($k,$v) = each(%hCorresp_file_species)) { | |
| 222 # print "file=$k sp=$v\n"; | |
| 223 # } | |
| 224 # exit; | |
| 225 | |
| 226 # STEP 3 : COMPTAGE DES CATEGORIES DE COG ------------------------------------------------------------------------------ | |
| 227 my %hCount2 = (); # HASH -> key1: catégorie de COG ; key2: espèce associée ; val: comptage | |
| 228 | |
| 229 # comptage du core-genome / des gènes spé / du génome accessoire | |
| 230 my %hCore_Count = (); # HASH -> key: catégorie de COG ; val: comptage (ce hash ne sera utilisé que pour le core-genome) | |
| 231 my %hSpecific_Count = (); # HASH -> key: catégorie de COG ; val: comptage | |
| 232 my %hAccessory_Count = (); # HASH -> key: catégorie de COG ; val: comptage | |
| 233 | |
| 234 # hash pour récupérer le gène | |
| 235 my %hCore_Cat = (); # HASH -> key: catégorie de COG ; val: gène | |
| 236 my %hAccessory_Cat = (); # HASH -> key: catégorie de COG ; val: gène | |
| 237 my %hSpecific_Cat = (); # HASH -> key: catégorie de COG ; val: gène | |
| 238 | |
| 239 # 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 | |
| 241 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 | |
| 243 | |
| 244 my %Cog_of_gene = (); | |
| 245 | |
| 246 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 | |
| 248 # print $esp."\n"; | |
| 249 # exit; | |
| 250 | |
| 251 my %hCount = (); # HASH -> key: catégorie de COG ; val: comptage | |
| 252 | |
| 253 | |
| 254 open (A, $file); # on va parcourir maintenant chaque fichier un à un | |
| 255 | |
| 256 while (my $line2 = <A>){ | |
| 257 | |
| 258 $line2 =~s/\n//g; $line2 =~s/\r//g; # on procède ligne par ligne | |
| 259 my @Genes = split('\t', $line2); | |
| 260 my $gene = $Genes[0]; | |
| 261 my $first_cat = $Genes[2]; | |
| 262 $Cog_of_gene{$gene} = $first_cat; | |
| 263 | |
| 264 for (my $j=2; $j <= $#Genes; $j++) { | |
| 265 my $cat = $Genes[$j]; # on récupère la ou les catégorie(s) de COG | |
| 266 $hCount{$cat}++; # pour la catégorie donnée, on incrémente son nb d'occurences | |
| 267 | |
| 268 if ($coregenes{$gene}){ # si le $gene fait bien partie du core-genome (donc de notre table de hash %coregenes) | |
| 269 $hCore_Count{$cat}++; # on incrémente le hash | |
| 270 $hCore_Cat{$cat}=$gene; # on récupère le nom du gène | |
| 271 } | |
| 272 if ($accessorygenes{$gene}){ # s'il fait partie des gènes accessoires | |
| 273 $hAccessory_Count{$cat}++; | |
| 274 $hAccessory_Cat{$cat}=$gene; | |
| 275 } | |
| 276 if ($specificgenes{$gene}){ # s'il fait partie des gènes spécifiques | |
| 277 $hSpecific_Count{$cat}++; | |
| 278 $hSpecific_Cat{$cat}=$gene; | |
| 279 } | |
| 280 # $hCount2{$cat}{$esp}++; # TABLE DE HASH AVEC CLES=CAT DE COG + ESPECE VAL=COMPTAGE | |
| 281 } | |
| 282 | |
| 283 } | |
| 284 close (A); | |
| 285 | |
| 286 # print "$file $esp\n=============\n"; | |
| 287 while (my ($k,$v) = each(%hCount)) { # parcours de la table de hash de comptage | |
| 288 # print "cat=$k nb=$v\n"; | |
| 289 $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 | |
| 290 # 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) | |
| 291 } | |
| 292 | |
| 293 # Récupérer les gènes du core-génome | |
| 294 while (my ($cat_core,$gene_core) = each(%hCore_Cat)) { | |
| 295 $hCore_Cat_Esp{$cat_core}{$esp}=$gene_core; | |
| 296 } | |
| 297 # Récupérer les gènes du génome-accessoire | |
| 298 while (my ($cat_acc,$gene_acc) = each(%hAccessory_Cat)) { | |
| 299 $hAccessory_Cat_Esp{$cat_acc}{$esp}=$gene_acc; | |
| 300 } | |
| 301 # Récupérer les gènes spécifique | |
| 302 while (my ($cat_spe,$gene_spe) = each(%hSpecific_Cat)) { | |
| 303 $hSpecific_Cat_Esp{$cat_spe}{$esp}=$gene_spe; | |
| 304 } | |
| 305 | |
| 306 } | |
| 307 # foreach my $category (sort keys (%hSpecific_Cat_Esp)) { # parcours au niveau de la 1ere clé | |
| 308 | |
| 309 # foreach my $especeee (keys %{$hSpecific_Cat_Esp{$category} }) { # parcours au niveau de la 2e clé pour la $category donnée | |
| 310 | |
| 311 # print "$category\t$especeee\t$hSpecific_Cat_Esp{$category}{$especeee}\n"; # on crée une sortie qui affiche en somme notre hash %hCount2 | |
| 312 # } | |
| 313 # } | |
| 314 # exit; | |
| 315 | |
| 316 # STEP 4 : AFFICHAGE DANS LE FICHIER DE SORTIE ------------------------------------------------------------------------------ | |
| 317 open (OUT4, ">$output4") or die $!; | |
| 318 | |
| 319 print OUT4 "Species"."\t"."COG categories"."\t"."Core-genome"."\t"."Accessory genome"."\t"."Specific genes"."\n"; | |
| 320 | |
| 321 foreach my $category (sort keys (%hCount2)){ # parcours de la table %hCount2 au niveau des catégories | |
| 322 foreach my $especeee (keys %{$hCount2{$category} }){ # parcours de la table %hCount2 au niveau des espèces | |
| 323 print OUT4 "$especeee\t$category\t"; # affichage des esp puis des cat | |
| 324 | |
| 325 # if ($hCore_Cat_Esp{$category}{$especeee}) { | |
| 326 # print OUT4 "$hCore_Cat_Esp{$category}{$especeee}\t"; | |
| 327 # } | |
| 328 my $c = 0; | |
| 329 if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome | |
| 330 $c = ($hCore_Count{$category}/scalar keys (%coregenes))*100; # calcul du % du comptage | |
| 331 } | |
| 332 print OUT4 "$c\t"; # affichage du % | |
| 333 | |
| 334 # if ($hAccessory_Cat_Esp{$category}{$especeee}) { | |
| 335 # print OUT4 "$hAccessory_Cat_Esp{$category}{$especeee}\t"; | |
| 336 # } | |
| 337 my $acc = 0; | |
| 338 if ($hAccessory_Count{$category}){ # si cette catégorie existe dans le génome accessoire | |
| 339 $acc = ($hAccessory_Count{$category}/scalar keys (%accessorygenes))*100; # calcul du % du comptage | |
| 340 } | |
| 341 print OUT4 "$acc\t"; # affichage du % | |
| 342 | |
| 343 # if ($hSpecific_Cat_Esp{$category}{$especeee}) { | |
| 344 # print OUT4 "$hSpecific_Cat_Esp{$category}{$especeee}\t"; | |
| 345 # } | |
| 346 my $s = 0; | |
| 347 if ($hSpecific_Count{$category}){ # si cette catégorie existe dans les gènes spécifiques | |
| 348 $s = ($hSpecific_Count{$category}/scalar keys (%specificgenes))*100; # calcul du % du comptage | |
| 349 } | |
| 350 print OUT4 "$s\n"; # affichage du % | |
| 351 } | |
| 352 } | |
| 353 close (OUT4); | |
| 354 | |
| 355 open (OUT3, ">$output3") or die $!; | |
| 356 foreach my $category (sort keys (%hCount2)) { # parcours au niveau de la 1ere clé | |
| 357 | |
| 358 foreach my $especeee (keys %{$hCount2{$category} }) { # parcours au niveau de la 2e clé pour la $category donnée | |
| 359 | |
| 360 print OUT3 "$category\t$especeee\t$hCount2{$category}{$especeee}\n"; # on crée une sortie qui affiche en somme notre hash %hCount2 | |
| 361 } | |
| 362 } | |
| 363 | |
| 364 close (OUT3); | |
| 365 | |
| 366 | |
| 367 open (OUT2, ">$output2") or die $!; | |
| 368 | |
| 369 print OUT2 "category"; | |
| 370 foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces... | |
| 371 print OUT2 "\t".$e; #... où on récupère le nom de celles-ci | |
| 372 } | |
| 373 print OUT2 "\n"; | |
| 374 | |
| 375 foreach my $category (sort keys (%hCount2)) { # on parcourt de nouveau les catégories de notre hash à 2 clés | |
| 376 print OUT2 $category; | |
| 377 | |
| 378 foreach my $especes (sort keys (%species_names)) { # on parcourt également le hash d'espèces | |
| 379 | |
| 380 my $nbr = 0; | |
| 381 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 | |
| 382 $nbr = $hCount2{$category}{$especes}; | |
| 383 } | |
| 384 # $nbr =~s/\n//g; $nbr =~s/\r//g; | |
| 385 | |
| 386 | |
| 387 my @liste = split(',', $nbr); # vu qu'il peut y avoir plusieurs nombres on les dissocie | |
| 388 | |
| 389 my $somme=0; | |
| 390 my $n=0; | |
| 391 my $moyenne=0; | |
| 392 #print "\nma liste de $nbr: ".join("%",@liste)."\n"; | |
| 393 foreach my $x (@liste) { # on parcourt nos nombres | |
| 394 $somme=$somme+$x; | |
| 395 $n=$n+1; | |
| 396 } | |
| 397 | |
| 398 if ($n>0){ | |
| 399 $moyenne = $somme/$n; # on fait le calcul de la moyenne | |
| 400 } | |
| 401 # print "$category, $especes: $hCount2{$category}{$especes}\t"; | |
| 402 # print "moyenne = $moyenne\n=============\n"; | |
| 403 | |
| 404 print OUT2 "\t".$moyenne; # fichier de sortie | |
| 405 } | |
| 406 print OUT2 "\n"; | |
| 407 } | |
| 408 | |
| 409 close (OUT2); | |
| 410 | |
| 411 # foreach my $cat (keys (%hCore_Cat)){ | |
| 412 # print OUT4 $c_gene."\t"; | |
| 413 # } | |
| 414 | |
| 415 | |
| 416 #////////////////////////////////////////////////////////////////////////////////////////////////// | |
| 417 | |
| 418 ############################################### GFF ############################################### | |
| 419 | |
| 420 my @list_gff = split(',', $annotation_GFF); # liste des différents fichiers GFF (qui se retrouvent dans le dossier Annotation Maker) | |
| 421 my @order_gff = split(',', $order_GFF); # liste de l'ordre des souches | |
| 422 my ($g,$o); | |
| 423 | |
| 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) | |
| 425 my %Gene_position = (); | |
| 426 my %Cat_genes = (); | |
| 427 | |
| 428 # my %hCore_Cat = reverse (%hCore_Cat); | |
| 429 | |
| 430 # ++++++++++++ parcours de 2 listes en même temps ++++++++++++ # | |
| 431 foreach $g (@list_gff){ | |
| 432 $hgff_order{$g} = $order_gff[$o++]; # on fait correspondre pour chaque fichier GFF, un nom de souche | |
| 433 open (G, $g); | |
| 434 while (<G>) { | |
| 435 my @table_gff = split (/\t/, $_); | |
| 436 my $chr = $table_gff[0]; | |
| 437 my $start = $table_gff[3]; | |
| 438 my $end = $table_gff[4]; | |
| 439 my $gene_name = $table_gff[8]; | |
| 440 my $type = $table_gff[2]; | |
| 441 | |
| 442 | |
| 443 | |
| 444 if ($type && $type eq "mRNA" && $gene_name =~ /ID=([^;]+);/){ | |
| 445 my $gene = $1; | |
| 446 # print $gene."\n"; | |
| 447 # exit; | |
| 448 foreach my $cog (keys (%hCore_Cat)){ | |
| 449 if ($hCore_Cat{$cog} eq $gene){ | |
| 450 $Cat_genes{$gene}=$cog; | |
| 451 } | |
| 452 } | |
| 453 | |
| 454 $Gene_position{$gene}="$chr\t$start\t$end"; | |
| 455 } | |
| 456 | |
| 457 | |
| 458 } | |
| 459 | |
| 460 close (G); | |
| 461 } | |
| 462 | |
| 463 open (OUT5, "> $output5"); | |
| 464 print OUT5 "Orthogroups\tGenes\tChromosomes\tStart\tEnd\tCOG categories\n"; | |
| 465 | |
| 466 | |
| 467 | |
| 468 | |
| 469 # foreach my $gene1 (keys (%Gene_position)){ | |
| 470 # if ($hCore_Cat{$gene1}){ | |
| 471 # # print "$hCore_Cat{$gene1}\n"; | |
| 472 # my $cog = $hCore_Cat{$gene1}; | |
| 473 # $Hash_genes{$gene1}=$cog; | |
| 474 # } | |
| 475 # } | |
| 476 # while (my ($k,$v) = each(%Hash_genes)) { | |
| 477 # print "gene=$k cog=$v\n"; | |
| 478 # } | |
| 479 # exit; | |
| 480 | |
| 481 foreach my $gene (keys (%coregenes)){ | |
| 482 # print "$gene\n"; | |
| 483 my $cat = "unknown"; | |
| 484 if ($Cog_of_gene{$gene}){ | |
| 485 $cat = $Cog_of_gene{$gene}; | |
| 486 } | |
| 487 print OUT5 $coregenes{$gene}."\t"."$gene\t"."\t".$Gene_position{$gene}.$cat."\n"; | |
| 488 | |
| 489 } | |
| 490 | |
| 491 close (OUT5); |