Mercurial > repos > mgarnier > pangenome_cog_analysis
changeset 13:45cc191a3290 draft
Uploaded
| author | mgarnier |
|---|---|
| date | Thu, 19 Aug 2021 13:39:14 +0000 |
| parents | b36506d26a43 |
| children | ef012cfc41e8 |
| files | pangenomeCogAnalysis.pl |
| diffstat | 1 files changed, 1242 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pangenomeCogAnalysis.pl Thu Aug 19 13:39:14 2021 +0000 @@ -0,0 +1,1242 @@ +#!/usr/bin/perl + +use strict; +use warnings; + +my $num_args = $#ARGV + 1; +if ($num_args != 11) { + print "Il n'y a pas le bon nombre d'arguments !\n"; + exit; +} + +# INPUT_ +my $matrix_file = $ARGV[0]; # fichier tabulé : une liste d'orthogroupes qui se retrouvent ou non dans les différentes souches +my $species_file = $ARGV[1]; # association de chaque souche à son espèce (fichier tabulé également) +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) +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 +my $annotation_GFF = $ARGV[4]; # fichiers avec les GFF +# my $order_GFF = $ARGV[5]; + +# OUTPUT_ +my $output = $ARGV[5]; # liste des espèces avec leurs orthogroupes (présence-absence) +my $output2 = $ARGV[6]; # fichier des moyennes +my $output3 = $ARGV[7]; # fichier de la liste des valeurs pour chaque catégorie de COG et pour chaque espèce +my $output4 = $ARGV[8]; # fichier avec les catégories de COG pour core-génome / génome accessoire / gènes spé +my $output7 = $ARGV[9]; +my $output8 = $ARGV[10]; +# my $output9 = $ARGV[11]; + + +#print "ok\n"; +#exit; + +my @list_gff = split(',', $annotation_GFF); # liste des différents fichiers GFF (qui se retrouvent dans le dossier Annotation Maker) +my %hSpecies = (); # HASH -> key: N_Id (ex NF_AR12) ; val: nom de l'esp (ex Naegleria Fowleri) + +######################## LE SPECIES_FILE ########################### +open (S, $species_file); +while (my $line = <S>){ + + $line =~s/\n//g; $line =~s/\r//g; + my @sp = split('\t', $line); + # print "$line\n"; + # exit; + $hSpecies{$sp[0]} = $sp[1]; # HASH -> key: N_Id ; val: name + +} + my $nbr = keys (%hSpecies); #compter le nombre de souches max + # = taille de la table de hash +# print "J'ai $nbr clés\n"; +# exit; + +close (S); + +#/////////////////////////////////////////////////////////////////////////////////////////////////// + +############################################ LA MATRICE ############################################ + +open(M, $matrix_file); + +my $first_line = <M>; +$first_line =~s/\n//g; $first_line =~s/\r//g; # ne garder que la première ligne du tableau +my @samples = split(/\t/,$first_line); # mettre dans une liste (@samples) chaque intitulé de colonne = N_Id +# print "$first_line\n"; +# exit; + +# Le but ici est de récupérer les combinaisons associées à chaque espèce : NF, NG et NL +my %hCombination =(); # HASH -> key: esp ; val: combinaison + +for (my $i=1; $i <= $#samples; $i++){ # on parcourt chaque colonne ($i) mais on ne regarde que le N_Id + my $header = $samples[$i]; # on récupère le N_Id dans $header (soit le nom de la colonne i) + 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 + $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 + # print "$header\n"; + # exit; +} + + +# foreach my $species (keys (%hCombination)){ +# my $combination = $hCombination{$species}; +# print "$species $combination\n"; +# } + + +# exit; + +# orthogrp présents : +my %hCombination_prs = (); # HASH -> key: combinaison ; val: liste des orthogroupes +# orthogrp absents : +my %hCombination_abs = (); # idem + + + +my %coregenes = (); # HASH -> key: gene ; val: orthogroupe (pour core-genome) +my %specificgenes = (); # HASH -> key: gene ; val: orthogroupe (pour gènes spécifiques) +my %accessorygenes = (); # HASH -> key: gene ; val: orthogroupe (pour génome accessoire) + +my $coregene_line; +my %coregenes2 = (); # HASH -> key1: colonne i ; key2: gène ; val: orthogroupe +my %specificgenes2 = (); # HASH -> key1: colonne i ; key2: gène ; val: orthogroupe + +my %Genes_of_OG = (); # HASH -> key1: orthogroupe ; key2: colonne i ; val: gène + +my %coregenes3 = (); #ligne complete +my %Type_count_byStrain = (); +my %OG_genes = (); +my $nb_genes_total = 0; +my %specificgenes3 = (); +my %Species_Total_Count = (); # HASH -> key: espèce ; val: comptage du nombre total de gènes pour cette espèce +my %Genes_Species_Total = (); +my %NonStrict_Spe = (); + +while(<M>) { + + my $line = $_; + $line =~s/\n//g; $line =~s/\r//g; + my $nb_found = 0; + my @infos = split(/\t/,$line); + my $orthogroup = $infos[0]; # on récupère le nom de l'orthogroupe dans $orthogroup + 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 + my $combi_prs = ""; + my $combi_abs = ""; + my $val; + my $gene_random; + my $unique_col_detected; + my %comptage_especes; + my $seule_espece; + + for (my $i=1; $i <= $#infos; $i++){ # on travaille par ligne puis dans chaque ligne (while(<M>)), cellule par cellule (cette boucle for) + + $val = $infos[$i]; # on récupère l'information contenue dans la case $i + + if ($val =~/\w/){ # s'il cette cellule contient qq chose... + $combi_prs .= "_".$i; # ...on va concaténer notre chaine $combi_prs pour que cela forme une combinaison + $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) + $gene_random=$val; # on récupère la valeur de la case (les gènes) + $unique_col_detected = $i; + my $espece = $hSpecies{$samples[$i]}; + + my @table_genes = split (',', $val); + my $premier_gene = $table_genes[0]; + # $Genes_of_OG{$i}{$orthogroup} = $premier_gene; # pour chaque orthorgoupe de chaque colonne, on récupère le premier gène + foreach my $genes (@table_genes){ + # $OG_genes{$orthogroup}{$genes} = 1; + $Genes_of_OG{$i}{$orthogroup} .= $genes; + $nb_genes_total++; + $Species_Total_Count{$espece}++; + $Genes_Species_Total{$genes} = $espece; + $seule_espece = $espece; + $comptage_especes{$espece} .= $genes; + + } + } + + else { # si jamais il n'y a rien dans la cellule... + $combi_abs .= "_".$i; # ... on fait la même chose mais avec $combi_abs + } + + } + + if (scalar keys(%comptage_especes) == 1){ + my $list = $comptage_especes{$seule_espece}; + my @table = split(" ",$list); + foreach my $gene (@table){ + $NonStrict_Spe{$gene} = $seule_espece; + } + # print $_; + # print $list."\n"; + } + + # $hCount{$combi}++; + $hCombination_prs{$combi_prs}.=$orthogroup."\n"; # à la fin de chaque ligne, on va ajouter notre orthogroupe à la combinaison qui lui correspond + $hCombination_abs{$combi_abs}.=$orthogroup."\n"; + + + + if ($nb_found == $#infos){ # si nb_found = au nombre de souche, c'est qu'on a à faire à un core-génome + # print "$orthogroup\n"; + # print "$nb_found\n=================\n"; + + for (my $i=1; $i <= $#infos; $i++){ + 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 ',') + foreach my $gene (@list_of_genes){ + $coregenes3{$samples[$i]}{$gene} = 1; + $Type_count_byStrain{"core"}{$samples[$i]}{"oui"}++; + $Type_count_byStrain{"accessory"}{$samples[$i]}{"non"}++; + $Type_count_byStrain{"unique"}{$samples[$i]}{"non"}++; + + } + + # $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) + # $coregenes2{$i}{$first_gene}= $orthogroup; + + } + my @liste_of_genes = split (',', $gene_random); + my $first_gene = $liste_of_genes[0]; + $coregenes{$first_gene}= $orthogroup; + + + if (!$coregene_line){ + $coregene_line = $line; + } + } + + elsif ($nb_found == 1) { # si on a un gène spé + + # # print "$gene_random\n"; + # # print "$line\n"; + # # print "$unique_col_detected\n"; + + + for (my $i=1; $i <= $#infos; $i++){ + 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 ',') + foreach my $gene (@list_of_genes){ + $specificgenes3{$samples[$i]}{$gene} = 1; + $Type_count_byStrain{"unique"}{$samples[$i]}{"oui"}++; + $Type_count_byStrain{"core"}{$samples[$i]}{"non"}++; + $Type_count_byStrain{"accessory"}{$samples[$i]}{"non"}++; + + } + + + } + # my @list_of_genes = split (',', $gene_random); # idem, on ne veut qu'un seul gène donc on crée la liste + # my $first_gene = $list_of_genes[0]; # on ne prend que le premier + # # print "$first_gene\n"; + # # exit; + # $specificgenes{$first_gene}= $orthogroup; # et pareil on crée la table de hash + # $specificgenes2{$unique_col_detected}{$first_gene}= $orthogroup; + } + + else { # là c'est le génome accessoire, i.e tout le reste ! + # for (my $i=1; $i <= $#infos; $i++){ + # 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 ',') + # my $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement ! + # $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) + + + # } + + for (my $i=1; $i <= $#infos; $i++){ + 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 ',') + foreach my $gene (@list_of_genes){ + # $coregenes3{$samples[$i]}{$gene} = 1; + $Type_count_byStrain{"accessory"}{$samples[$i]}{"oui"}++; + $Type_count_byStrain{"core"}{$samples[$i]}{"non"}++; + $Type_count_byStrain{"unique"}{$samples[$i]}{"non"}++; + } + + + + } + + + my @liste_of_genes = split (',', $gene_random); + my $first_gene = $liste_of_genes[0]; + $accessorygenes{$first_gene}= $orthogroup; + } + +} + +#print scalar keys(%Genes_of_OG);exit; +# print "$nb_genes_total\n"; + +# foreach my $og (keys %OG_genes) { +# foreach my $gene (keys %{$OG_genes{$og}}) { +# # print "$og\t$gene\n"; +# print $OG_genes{$og}."\n"; +# } +# } +# exit; +# foreach my $gene (keys (%Genes_Species_Total)) { +# print "$gene => ".$Genes_Species_Total{$gene}."\n"; +# } +# foreach my $strain (keys %specificgenes3) { +# foreach my $gene (keys %{$specificgenes3{$strain}}) { +# print "$strain\t$gene\n"; +# } +# } +# foreach my $gene (keys (%NonStrict_Spe)){ +# print $NonStrict_Spe{$gene}."\t$gene \n"; +# } +# exit; + +my %hCol_Annotated = (); # HASH -> key: colonne ; val: 1 (colonnes pour lesquelles les GFF sont présents) + +# Le but ici est de ne garder que les colonnes (donc les souches) qui ont un fichier GFF associé +my @list_column = split ('\t', $coregene_line); +for (my $i=1; $i <= $#list_column; $i++){ + my @list_genes = split (', ', $list_column[$i]); + my $premier_gene = $list_genes[0]; + my $strain = $samples[$i]; # récupérer le nom de la souche + + + foreach my $gff (@list_gff){ + my $result_grep = `grep $premier_gene $gff`; + + if ($result_grep){ + $hCol_Annotated{$i}=$strain; + + } + # print "$result_grep\n"; + } +} +# exit; +# foreach my $i (sort keys (%specificgenes2)){ # parcours de la table %hCount2 au niveau des catégories +# foreach my $gene (keys %{$specificgenes2{$i} }){ # parcours de la table %hCount2 au niveau des espèces +# print "$i\t$gene\t".$specificgenes2{$i}{$gene}."\n"; +# } +# } +# exit; +# while (my ($k,$v) = each(%accessorygenes)) { +# print "gene=$k OG=$v\n"; +# } +# exit; +# foreach my $oups (keys (%coregenes)) { +# print "$oups\n"; +# } +# exit; + +close (M); + +my %Hash_Specific = (); # HASH -> key: orthogroupe ; val: espèce + +open (OUT, '>', $output) or die $!; +# print OUT "$annotation\n"; +foreach my $species (keys (%hCombination)){ # parcours de la table de hash %hCombination (key: nom esp ; val: combi) + 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 + 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 + 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 + + # open (OUT,">results.list.txt"); + + if ($ortho_presents){ + print OUT "> $species - present\n"; + print OUT "$ortho_presents\n"; + my @orthogroups_name = split ('\n', $ortho_presents); + foreach my $ortho (@orthogroups_name){ + $Hash_Specific{$ortho} = $species; + } + } + + if ($ortho_absents){ + # open (OUT2,">$species.$combination.absents.list.txt"); + print OUT "> $species - absent\n"; + print OUT "$ortho_absents\n"; + } + +# close(OUT2); +} + +close(OUT); + +my %Gene_Specie_Spe = (); # HASH -> key: gène spé ; val: espèce +my %Species_Spe_Count = (); # HASH -> key: espèce ; val: comptage du nombre de gènes spécifiques à cette espèce + +foreach my $i (keys(%Genes_of_OG)){ + foreach my $ortho (keys %{$Genes_of_OG{$i}}){ + my $gene = $Genes_of_OG{$i}{$ortho}; + + if ($Hash_Specific{$ortho}){ + my $specie = $Hash_Specific{$ortho}; + + my @liste_genes = split(' ',$gene); + foreach my $g(@liste_genes){ + $Gene_Specie_Spe{$g} = $specie; + $Species_Spe_Count{$specie}++; + + } + + } + } +} + + +# exit; + +# my @table_keys = (); +my $nb_groupSpe_genes = 0; + +foreach my $gene (keys (%Gene_Specie_Spe)) { + my @table_keys = split (' ', $gene); + foreach my $unique_gene (@table_keys) { + $nb_groupSpe_genes++; + } +} + +# print scalar keys (%Gene_Specie_Spe)."\n"; +# while (my ($k,$v) = each(%Gene_Specie_Spe)) { +# if ($v =~/ruberi/) { +# print "gene=$k espece=$v\n"; +# } +# } +# foreach my $sp (keys (%Species_Spe_Count)){ +# print "$sp => ".$Species_Spe_Count{$sp}."\n"; +# } +# exit; +#////////////////////////////////////////////////////////////////////////////////////////////////// + +############################################### COG ############################################### + +# STEP 1 : CORRESPONDANCE ENTRE LES DIFFERENTS FICHIERS DE COG ET L'ORDRE -------------------------------------------- +my @files = split(',', $annotation); # liste des différents fichiers COG (qui se retrouvent dans le dossier Naegleria) +my @list = split(',', $order); # liste de l'ordre des souches +#my ($f,$l); + +my %hCorrespondance = (); #HASH -> key: un fichier COG ; val: un nom de souche (ces 2 données sont entrées en input = $annotation et $order) + +# ++++++++++++ parcours de 2 listes en même temps ++++++++++++ # +my $l = 1; +foreach my $f (@files){ + $hCorrespondance{$f} = $list[$l]; # on fait correspondre pour chaque fichier de COG, un nom de souche + $l++; +} + + + + +# #Affichage du hash +# foreach my $f (keys %hCorrespondance){ +# print $f."=>".$hCorrespondance{$f}."\n" +# } +# exit; + +# STEP 2 : POUR CHAQUE FICHIER DE COG, FAIRE CORRESPONDRE L'ESPECE (ET NON LA SOUCHE) ------------------------------------- +my %hCorresp_file_species = (); # HASH -> key: un fichier de COG ; val: une espèce +my %species_names; # HASH -> key: nom d'espèce ; val: 1 + +foreach my $h (keys (%hCorrespondance)){ # parcours de la table de hash {fichier COG => nom souche} + my $smpl = $hCorrespondance{$h}; # $smpl prend la valeur de la clé (donc d'un nom de souche) + 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 + $species_names{$espece} = 1; # on garde sous le coude nos nom d'espèce dans cette nouvelle table de hash + $hCorresp_file_species{$h} = $espece; # BUT ATTEINT : on donne pour chaque fichier de COG le nom de l'espèce qui lui correspond +} +# while (my ($k,$v) = each(%hCorresp_file_species)) { +# print "file=$k sp=$v\n"; +# } +# exit; + + + + +# STEP 3 : COMPTAGE DES CATEGORIES DE COG ------------------------------------------------------------------------------ +my %hCount2 = (); # HASH -> key1: catégorie de COG ; key2: espèce associée ; val: comptage + +# comptage du core-genome / des gènes spé / du génome accessoire +my %hCore_Count = (); # HASH -> key: catégorie de COG ; val: comptage (ce hash ne sera utilisé que pour le core-genome) +my %hSpecific_Count = (); # HASH -> key: catégorie de COG ; val: comptage +my %hAccessory_Count = (); # HASH -> key: catégorie de COG ; val: comptage + +# hash pour récupérer le gène +my %hCore_Cat = (); # HASH -> key: catégorie de COG ; val: gène +my %hAccessory_Cat = (); # HASH -> key: catégorie de COG ; val: gène +my %hSpecific_Cat = (); # HASH -> key: catégorie de COG ; val: gène + +# hash pour récupérer le gène +my %hCore_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène +my %hAccessory_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène +my %hSpecific_Cat_Esp = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: gène + +my %Acc_Cat_Esp_Count = (); # HASH -> key1: catégorie de COG ; key2: espèce ; val: comptage + +my %Cog_of_gene = (); # HASH -> key: gène ; val: cat de COG +my %Cogs_of_gene = (); # HASH -> key: gène ; val: cat de COG (plusieurs) +my %Specie_of_gene = (); # HASH -> key: gène ; val: souche + +my %Global_Count = (); +my %Species_Count = (); +my %Species_NonStrictSpe_Count = (); + +my %Genes_in_COG = (); +my %Count_Spe_Genes = (); +my %Count_Total_Species = (); +my %Count_NonStrictSpe_Genes = (); + +my %Nveau = (); + + +foreach my $file(@files){ # parcours de la liste des fichiers + my $esp = $hCorresp_file_species{$file}; # on récupère l'espèce pour chaque fichier de COG dans $esp + # print $esp."\n"; + # exit; + + my %hCount = (); # HASH -> key: catégorie de COG ; val: comptage + + + open (A, $file); # on va parcourir maintenant chaque fichier un à un + my $strain = $hCorrespondance{$file}; + + while (my $line2 = <A>){ + + $line2 =~s/\n//g; $line2 =~s/\r//g; # on procède ligne par ligne + my @Genes = split('\t', $line2); + my $cogs = $line2; + my $gene = $Genes[0]; + my $cog_id = $Genes[1]; + $cogs =~s/$gene//g; $cogs =~s/$cog_id//g; + my $first_cat = $Genes[2]; + $Cog_of_gene{$gene} = $first_cat; + $Cogs_of_gene{$gene} = $cogs; + + $Genes_in_COG{$gene} = $esp; + + for (my $j=2; $j <= $#Genes; $j++) { + my $cat = $Genes[$j]; # on récupère la ou les catégorie(s) de COG + $hCount{$cat}++; # pour la catégorie donnée, on incrémente son nb d'occurences + + + if ($coregenes{$gene}){ # si le $gene fait bien partie du core-genome (donc de notre table de hash %coregenes) + $hCore_Count{$cat}++; # on incrémente le hash + $hCore_Cat{$cat}=$gene; # on récupère le nom du gène + } + + if ($accessorygenes{$gene}){ # s'il fait partie des gènes accessoires + + $hAccessory_Count{$cat}++; + $hAccessory_Cat{$cat}=$gene; + + # if ($accessorygenes{$gene} && $Gene_Specie_Spe{$gene}){ + # print "$gene\n"; + # # my $espece = $Gene_Specie_Spe{$gene}; + # # print "$espece\n"; + # # $Nveau{$cat}{$espece}++; + # } + + } + if ($coregenes3{$strain}{$gene}){ + $Global_Count{"core"}{$cat}{$strain}{"oui"}++; + $Global_Count{"accessory"}{$cat}{$strain}{"non"}++; + $Global_Count{"unique"}{$cat}{$strain}{"non"}++; + } + elsif ($specificgenes3{$strain}{$gene}){ + $Global_Count{"unique"}{$cat}{$strain}{"oui"}++; + $Global_Count{"core"}{$cat}{$strain}{"non"}++; + $Global_Count{"accessory"}{$cat}{$strain}{"non"}++; + } + else { + $Global_Count{"accessory"}{$cat}{$strain}{"oui"}++; + $Global_Count{"core"}{$cat}{$strain}{"non"}++; + $Global_Count{"unique"}{$cat}{$strain}{"non"}++; + } + + + + + + + if ($Gene_Specie_Spe{$gene}) { + $Species_Count{$esp}{$cat}{$strain}{"oui"}++; + } + else { + $Species_Count{$esp}{$cat}{$strain}{"non"}++; + } + + + + if ($NonStrict_Spe{$gene}) { + $Species_NonStrictSpe_Count{$esp}{$cat}{$strain}{"oui"}++; + } + else { + $Species_NonStrictSpe_Count{$esp}{$cat}{$strain}{"non"}++; + } + + # $Global_Count{$cat}{"accessory"}{$strain}++; + + + # if ($specificgenes{$gene}){ # s'il fait partie des gènes spécifiques + # $hSpecific_Count{$cat}++; + # $hSpecific_Cat{$cat}=$gene; + # } + # $hCount2{$cat}{$esp}++; # TABLE DE HASH AVEC CLES=CAT DE COG + ESPECE VAL=COMPTAGE + + + } + } + + + + close (A); + + # foreach my $espece (sort keys (%Species_NonStrictSpe_Count)) { + # foreach my $cat (sort keys %{$Species_NonStrictSpe_Count{$espece}}) { + # foreach my $strain (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}}) { + # foreach my $choix (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}{$strain}}) { + + # print "$espece - $cat - $strain - $choix ". $Species_NonStrictSpe_Count{$espece}{$cat}{$strain}{$choix}."\n"; + # } + # } + + # } + # } + # exit; + + # while (my ($k,$v) = each(%hCore_Cat)) { + # print "cat=$k gene=$v\n"; + # } + # exit; + + # print "$file $esp\n=============\n"; + while (my ($k,$v) = each(%hCount)) { # parcours de la table de hash de comptage + # print "cat=$k nb=$v\n"; + $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 + # 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) + } + + # Récupérer les gènes du core-génome + while (my ($cat_core,$gene_core) = each(%hCore_Cat)) { + $hCore_Cat_Esp{$cat_core}{$esp}=$gene_core; + } + # Récupérer les gènes du génome-accessoire + while (my ($cat_acc,$gene_acc) = each(%hAccessory_Cat)) { + $hAccessory_Cat_Esp{$cat_acc}{$esp}=$gene_acc; + } + # Récupérer les gènes spécifique + while (my ($cat_spe,$gene_spe) = each(%hSpecific_Cat)) { + $hSpecific_Cat_Esp{$cat_spe}{$esp}=$gene_spe; + } + +### + while (my ($cat,$count) = each(%hAccessory_Count)) { + $Acc_Cat_Esp_Count{$cat}{$esp}=$count; + } +} +# foreach my $type (sort keys (%Global_Count)) { +# foreach my $cat (sort keys %{$Global_Count{$type}}) { +# foreach my $strain (sort keys %{$Global_Count{$type}{$cat}}) { +# foreach my $choix (sort keys %{$Global_Count{$type}{$cat}{$strain}}) { + +# print "$type - $cat - $strain - $choix ". $Global_Count{$type}{$cat}{$strain}{$choix}."\n"; +# } +# } + +# } +# } +# foreach my $espece (sort keys (%Global_Count)) { +# foreach my $cat (sort keys %{$Species_NonStrictSpe_Count{$espece}}) { +# foreach my $strain (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}}) { +# foreach my $choix (sort keys %{$Species_NonStrictSpe_Count{$espece}{$cat}{$strain}}) { + +# print "$espece - $cat - $strain - $choix ". $Species_NonStrictSpe_Count{$espece}{$cat}{$strain}{$choix}."\n"; +# } +# } + +# } +# } +# exit; + +foreach my $gene (keys (%Genes_in_COG)){ + my $espece = $Genes_in_COG{$gene}; + if ($Gene_Specie_Spe{$gene}) { + $Count_Spe_Genes{$espece}++; + } + if ($Genes_Species_Total{$gene}) { + $Count_Total_Species{$espece}++; + } + if ($NonStrict_Spe{$gene}) { + $Count_NonStrictSpe_Genes{$espece}++; + } +} + + +# ############################################# +# # p / (1-p) p * (1-q) # +# # odds ratio = ----------- = ----------- # +# # q / (1-q) q * (1-p) # +# ############################################# +# # où p : proba qu'un E arrive au groupe A +# # où q : proba que ce même E arrive au groupe B + + +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"); +################### +open (OUT7, ">$output7") or die $!; + +# my $nb_files = scalar keys @files; + + +print OUT7 "\t"; + +foreach my $category(@orders){ +# foreach my $category (sort keys (%Acc_Cat_Esp_Count)) { + # my $cat = $category."\t"; + print OUT7 $category."\t"; +} + +print OUT7 "\n"; + + +# foreach my $category (sort keys (%Global_Count)){ +foreach my $type (sort keys (%Global_Count)){ + + print OUT7 "$type\t"; + #foreach my $category (sort keys (%{$Global_Count{$type}})){ + foreach my $category(@orders){ + + foreach my $strain (sort keys (%{$Global_Count{$type}{$category}})){ + + # foreach my $type (sort keys (%{$Global_Count{$category}{$strain}})){ + my $nb_type1; my $nb_type2; + + if ($Global_Count{$type}{$category}{$strain}{"non"} && $Global_Count{$type}{$category}{$strain}{"oui"}) { + $nb_type1 = $Type_count_byStrain{$type}{$strain}{"non"} - $Global_Count{$type}{$category}{$strain}{"non"}; + $nb_type2 = $Type_count_byStrain{$type}{$strain}{"oui"} - $Global_Count{$type}{$category}{$strain}{"oui"}; + } + # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$type}{$strain}."\t"."$nb_type\n"; + my $ratio1; my $ratio2; + if ($nb_type1 && $nb_type2) { + $ratio1 = $Global_Count{$type}{$category}{$strain}{"non"}/ $nb_type1; + $ratio2 = $Global_Count{$type}{$category}{$strain}{"oui"} / $nb_type2; + } + my $odds_ratio; + + if ($ratio1 && $ratio2) { + $odds_ratio = $ratio2 / $ratio1; + } + # print $strain." ".$Type_count_byStrain{$type}{$strain}{"oui"}."\n"; + # 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"; + if ($odds_ratio) { + print OUT7 "$odds_ratio;"; + } + + } + print OUT7 "\t"; + } + + print OUT7 "\n"; +} + + +print OUT7 "\n"; +close (OUT7); + + +#////////////////////////////////////////////// +open (OUT8, ">$output8") or die $!; + + + +print OUT8 "\t"; + +# +foreach my $category(@orders){ + print OUT8 $category."\t"; +} + +print OUT8 "\n"; + + +# foreach my $category (sort keys (%Global_Count)){ +foreach my $specie (sort keys (%Species_Count)){ + # my $nb_genes_nonSpe = $Species_Total_Count{$specie} - $Species_Spe_Count{$specie}; + my $nb_genes_nonSpe = $Count_Total_Species{$specie} - $Count_Spe_Genes{$specie}; + + print OUT8 "$specie\t"; + foreach my $category (sort keys (%{$Species_Count{$specie}})){ + + foreach my $strain (sort keys (%{$Species_Count{$specie}{$category}})){ + + + my $nb_type1; my $nb_type2; + + if ($Species_Count{$specie}{$category}{$strain}{"non"} && $Species_Count{$specie}{$category}{$strain}{"oui"}) { + $nb_type1 = $nb_genes_nonSpe - $Species_Count{$specie}{$category}{$strain}{"non"}; # 1-q + $nb_type2 = $Count_Spe_Genes{$specie} - $Species_Count{$specie}{$category}{$strain}{"oui"}; # 1-p + # $nb_type2 = $Species_Spe_Count{$specie} - $Species_Count{$specie}{$category}{$strain}{"oui"}; # 1-p + # 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; + } + + # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$type}{$strain}."\t"."$nb_type\n"; + my $ratio1; my $ratio2; + if ($nb_type1 && $nb_type2) { + $ratio1 = $Species_Count{$specie}{$category}{$strain}{"non"}/ $nb_type1; + $ratio2 = $Species_Count{$specie}{$category}{$strain}{"oui"} / $nb_type2; + } + my $odds_ratio; + + if ($ratio1 && $ratio2) { + $odds_ratio = $ratio2 / $ratio1; + } + # print $strain." ".$Type_count_byStrain{$type}{$strain}{"oui"}."\n"; + # 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"; + if ($odds_ratio) { + print OUT8 "$odds_ratio;"; + } + + } + print OUT8 "\t"; + } +print OUT8 "\n"; +} + +# print OUT9 "\n"; + +close (OUT8); + +#/////////////////////////////////////////////////// +# open (OUT9, '>', $output9) or die $!; + + + +# print OUT9 "\t"; + +# # +# foreach my $category(@orders){ +# print OUT9 $category."\t"; +# } + +# print OUT9 "\n"; + + +# # foreach my $category (sort keys (%Global_Count)){ +# foreach my $specie (sort keys (%Species_NonStrictSpe_Count)){ +# # my $nb_genes_nonSpe = $Species_Total_Count{$specie} - $Species_Spe_Count{$specie}; +# my $nb_genes_nonSpeNS = $Count_Total_Species{$specie} - $Count_NonStrictSpe_Genes{$specie}; + +# print OUT9 "$specie\t"; +# foreach my $category (sort keys (%{$Species_NonStrictSpe_Count{$specie}})){ + +# foreach my $strain (sort keys (%{$Species_NonStrictSpe_Count{$specie}{$category}})){ + + +# my $nb_type1; my $nb_type2; + +# if ($Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"non"} && $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"oui"}) { +# $nb_type1 = $nb_genes_nonSpeNS - $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"non"}; # 1-q +# $nb_type2 = $Count_NonStrictSpe_Genes{$specie} - $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"oui"}; # 1-p +# # $nb_type2 = $Species_Spe_Count{$specie} - $Species_Count{$specie}{$category}{$strain}{"oui"}; # 1-p +# # 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; +# } + +# # print OUT8 "$category\t$type\t$strain\t".$Global_Count{$category}{$type}{$strain}."\t"."$nb_type\n"; +# my $ratio1; my $ratio2; +# if ($nb_type1 && $nb_type2) { +# $ratio1 = $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"non"}/ $nb_type1; +# $ratio2 = $Species_NonStrictSpe_Count{$specie}{$category}{$strain}{"oui"} / $nb_type2; +# } +# my $odds_ratio; + +# if ($ratio1 && $ratio2) { +# $odds_ratio = $ratio2 / $ratio1; +# } +# # print $strain." ".$Type_count_byStrain{$type}{$strain}{"oui"}."\n"; +# # 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"; +# if ($odds_ratio) { +# print OUT9 "$odds_ratio;"; +# } + +# } +# print OUT9 "\t"; +# } +# print OUT9 "\n"; +# } + +# # print OUT9 "\n"; + +# close (OUT9); + +# exit; +########################## sortie de pourcentages ########################## +# my $somme_core = 0; +# my $somme_acc = 0; + + +# foreach my $cat (keys(%hCore_Count)){ +# $somme_core = $somme_core + $hCore_Count{$cat}; +# } +# foreach my $category (sort keys (%Acc_Cat_Esp_Count)) { + +# foreach my $especeee (keys %{$Acc_Cat_Esp_Count{$category}}) { +# $somme_acc = $somme_acc + $Acc_Cat_Esp_Count{$category}{$especeee}; + +# } + +# } + + +# print "COG categories\tCore-genome\tAccessory genome\n"."\n"; +# # foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces... +# # print $e."\t"; #... où on récupère le nom de celles-ci +# # } +# # print "\n"; + +# foreach my $category (sort keys (%Acc_Cat_Esp_Count)) { # parcours au niveau de la 1ere clé +# my $nb_core = 0; +# my $somme_totale = 0; +# my $number = 0; +# print $category."\t"; +# my $c = 0; +# if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome +# $c = $hCore_Count{$category}; +# # $hash_core_pc{$c} = 1; +# $somme_totale = $somme_totale + $c; +# $nb_core = ($c/$somme_core)*100; +# # print "$nb_core\t"; +# } + +# foreach my $especes (sort keys (%species_names)) { +# my $nb_acc = 0; +# my $acc = 0; + +# 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 +# $acc = $Acc_Cat_Esp_Count{$category}{$especes}; +# $number = $number + $acc; +# $somme_totale = $somme_totale + $acc; +# $nb_acc = ($acc/$somme_acc)*100; +# # print "$nb_acc\t"; +# } + +# } + + +# print "|\t"; +# my $pourcentage_core = ($c/$somme_totale)*100; +# print "$pourcentage_core\t"; +# my $pourcentage_acc = ($number/$somme_totale)*100; +# print "$pourcentage_acc\n"; + +# } + +### +# exit; + + +########################## sortie de comptage ########################## +# print "COG categories\tCore-genome\tAccessory genome\n"."\t\t"; +# foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces... +# print $e."\t"; #... où on récupère le nom de celles-ci +# } +# print "\n"; + +# foreach my $category (sort keys (%Acc_Cat_Esp_Count)) { # parcours au niveau de la 1ere clé +# print $category."\t\t"; +# my $c = 0; +# if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome +# $c = $hCore_Count{$category}; +# print "$c\t"; +# } +# foreach my $especes (sort keys (%species_names)) { +# 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 +# print $Acc_Cat_Esp_Count{$category}{$especes}."\t"; + +# } +# } + + +# print "\n"; +# } + +### + +# foreach my $category (sort keys (%hCount2)) { # on parcourt de nouveau les catégories de notre hash à 2 clés +# print OUT2 $category; + +# foreach my $especes (sort keys (%species_names)) { # on parcourt également le hash d'espèces + +# my $nbr = 0; +# 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 +# $nbr = $hCount2{$category}{$especes}; +# } +# STEP 4 : AFFICHAGE DANS LE FICHIER DE SORTIE ------------------------------------------------------------------------------ +open (OUT4, ">$output4") or die $!; + +print OUT4 "COG categories"."\t"."Core-genome"."\t"."Accessory genome"."\t"."Strain Specific genes"."\n"; + +foreach my $category (sort keys (%hCount2)){ # parcours de la table %hCount2 au niveau des catégories + my $c = 0; + if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome + $c = $hCore_Count{$category}; + } + my $acc = 0; + if ($hAccessory_Count{$category}){ # si cette catégorie existe dans le génome accessoire + $acc = $hAccessory_Count{$category}; + } + # my $s = 0; + # if ($hSpecific_Count{$category}){ # si cette catégorie existe dans les gènes spécifiques + # $s = $hSpecific_Count{$category}; + # } + print OUT4 "$category\t".$c."\t".$acc."\n";#.$s."\n"; + + foreach my $especeee (keys %{$hCount2{$category} }){ # parcours de la table %hCount2 au niveau des espèces + # print OUT4 "$especeee\t$category\t"; # affichage des esp puis des cat + + # if ($hCore_Cat_Esp{$category}{$especeee}) { + # print OUT4 "$hCore_Cat_Esp{$category}{$especeee}\t"; + # } + my $c = 0; + if ($hCore_Count{$category}){ # si cette catégorie existe dans le core-génome + $c = ($hCore_Count{$category}/scalar keys (%coregenes))*100; # calcul du % du comptage + } + # print OUT4 "$c\t"; # affichage du % + + # if ($hAccessory_Cat_Esp{$category}{$especeee}) { + # print OUT4 "$hAccessory_Cat_Esp{$category}{$especeee}\t"; + # } + my $acc = 0; + if ($hAccessory_Count{$category}){ # si cette catégorie existe dans le génome accessoire + $acc = ($hAccessory_Count{$category}/scalar keys (%accessorygenes))*100; # calcul du % du comptage + } + # print OUT4 "$acc\t"; # affichage du % + + # # if ($hSpecific_Cat_Esp{$category}{$especeee}) { + # # print OUT4 "$hSpecific_Cat_Esp{$category}{$especeee}\t"; + # # } + # my $s = 0; + # if ($hSpecific_Count{$category}){ # si cette catégorie existe dans les gènes spécifiques + # $s = ($hSpecific_Count{$category}/scalar keys (%specificgenes))*100; # calcul du % du comptage + # } + # # print OUT4 "$s\n"; # affichage du % + } +} +close (OUT4); + +open (OUT3, ">$output3") or die $!; +foreach my $category (sort keys (%hCount2)) { # parcours au niveau de la 1ere clé + + foreach my $especeee (keys %{$hCount2{$category} }) { # parcours au niveau de la 2e clé pour la $category donnée + + print OUT3 "$category\t$especeee\t$hCount2{$category}{$especeee}\n"; # on crée une sortie qui affiche en somme notre hash %hCount2 + } +} + +close (OUT3); + + +open (OUT2, ">$output2") or die $!; + +print OUT2 "category"; +foreach my $e (sort keys (%species_names)){ # on parcours le hash d'espèces... + print OUT2 "\t".$e; #... où on récupère le nom de celles-ci +} +print OUT2 "\n"; + +foreach my $category (sort keys (%hCount2)) { # on parcourt de nouveau les catégories de notre hash à 2 clés + print OUT2 $category; + + foreach my $especes (sort keys (%species_names)) { # on parcourt également le hash d'espèces + + my $nbr = 0; + 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 + $nbr = $hCount2{$category}{$especes}; + } + # $nbr =~s/\n//g; $nbr =~s/\r//g; + + + my @liste = split(',', $nbr); # vu qu'il peut y avoir plusieurs nombres on les dissocie + + my $somme=0; + my $n=0; + my $moyenne=0; + #print "\nma liste de $nbr: ".join("%",@liste)."\n"; + foreach my $x (@liste) { # on parcourt nos nombres + $somme=$somme+$x; + $n=$n+1; + } + + if ($n>0){ + $moyenne = $somme/$n; # on fait le calcul de la moyenne + } + # print "$category, $especes: $hCount2{$category}{$especes}\t"; + # print "moyenne = $moyenne\n=============\n"; + + print OUT2 "\t".$moyenne; # fichier de sortie + } +print OUT2 "\n"; +} + +close (OUT2); + +# foreach my $cat (keys (%hCore_Cat)){ +# print OUT4 $c_gene."\t"; +# } + + +#////////////////////////////////////////////////////////////////////////////////////////////////// + +############################################### GFF ############################################### + + +# my @order_gff = split(',', $order_GFF); # liste de l'ordre des souches +my ($g,$o); + +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) +my %Gene_position = (); +my %Cat_genes = (); +my %Cat_genes2 = (); + +my %hash_of_genes = (); + + +foreach $g (@list_gff){ + # print "$g\n"; + # $hgff_order{$g} = $order_gff[$o++]; # on fait correspondre pour chaque fichier GFF, un nom de souche + open (G, $g); + while (<G>) { + my @table_gff = split (/\t/, $_); + my $chr = $table_gff[0]; + my $start = $table_gff[3]; + my $end = $table_gff[4]; + my $gene_name = $table_gff[8]; + my $type = $table_gff[2]; + + + + if ($type && $type eq "mRNA" && $gene_name =~ /ID=([^;]+);/){ #or $type eq "CDS" + my $gene = $1; + # print $gene."\n"; + # exit; + $hash_of_genes{$gene}=1; + + foreach my $cog (keys (%hCore_Cat)){ + if ($hCore_Cat{$cog} eq $gene){ + $Cat_genes{$gene}=$cog; + } + } + foreach my $cog_bis (keys (%hSpecific_Cat)){ + if ($hSpecific_Cat{$cog_bis} eq $gene){ + $Cat_genes2{$gene}=$cog_bis; + } + } + + $Gene_position{$gene}="$chr\t$start\t$end"; + } + + # foreach my $gene (keys (%hash_of_genes)){ + # my $orthogrp = $hGene_OG{$gene}; + # print "$orthogrp\n"; + # } + } + + close (G); +} + +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); + +mkdir("Core"); +foreach my $i (keys (%coregenes2)){ + + if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas + next; + } + + + my $strain_name = $hCol_Annotated{$i}; + + my $specie_name = $hSpecies{$strain_name}; + + + + open (OUT5, "> Core/$strain_name.$specie_name.txt") or die "Cannot create file $!\n"; + print OUT5 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n"; + + my $refcoregenes2 = $coregenes2{$i}; + my %subhash = %$refcoregenes2; + foreach my $gene (keys (%subhash)){ + # print "$gene\n"; + my $cat = "unknown"; + if ($Cog_of_gene{$gene}){ + $cat = $Cog_of_gene{$gene}; + } + # if (!$Gene_position{$gene}){ + # print "$gene\n coucou"; exit; + # } + + # if (!$subhash{$gene}){ + # print "$gene\n"; + # } + print OUT5 $subhash{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n"; + + } + + close (OUT5); +} + +mkdir("StrainSpecific"); +foreach my $i (keys (%specificgenes2)){ + + if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas + next; + } + + + my $strain_name = $hCol_Annotated{$i}; + + my $specie_name = $hSpecies{$strain_name}; + + + + open (OUT7, "> StrainSpecific/$strain_name.$specie_name.txt") or die "Cannot create file $!\n"; + print OUT7 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n"; + + my $refspecificgenes2 = $specificgenes2{$i}; + my %subhash = %$refspecificgenes2; + foreach my $gene (keys (%subhash)){ + # print "$gene\n"; exit; + my $cat = "unknown"; + if ($Cog_of_gene{$gene}){ + $cat = $Cog_of_gene{$gene}; + } + # if (!$Gene_position{$gene}){ + # print "$gene\n coucou"; exit; + # } + + # if (!$subhash{$gene}){ + # print "$gene\n"; + # } + print OUT7 $subhash{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n"; + + } + + close (OUT7); +} + + +mkdir("GroupSpecific"); +foreach my $i (keys (%Genes_of_OG)){ + if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas + next; + } + + my $strain_name = $hCol_Annotated{$i}; + + my $specie_name = $hSpecies{$strain_name}; + + open (OUT6, "> GroupSpecific/$strain_name.$specie_name.txt") or die "Cannot create file $!\n"; + print OUT6 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\tNumber assigned\n"; + + my $refGenes_of_OG = $Genes_of_OG{$i}; + my %subhash = %$refGenes_of_OG; + + foreach my $orthogroup (keys (%subhash)){ + if ($Hash_Specific{$orthogroup} && $Hash_Specific{$orthogroup} eq $specie_name){ + my $gene = $subhash{$orthogroup}; + + my $cat = "unknown"; + if ($Cog_of_gene{$gene}){ + $cat = $Cog_of_gene{$gene}; + } + print OUT6 $orthogroup."\t".$subhash{$orthogroup}."\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n"; + } + + } + close (OUT6); +}