Mercurial > repos > mgarnier > pangenome_cog_analysis

--- a/pangenomeCogAnalysis.pl	Thu Aug 19 19:12:05 2021 +0000
+++ b/pangenomeCogAnalysis.pl	Thu Aug 19 19:12:20 2021 +0000
@@ -4,7 +4,7 @@
 use warnings;

 my $num_args = $#ARGV + 1;
-if ($num_args != 11) {
+if ($num_args != 10) {
     print "Il n'y a pas le bon nombre d'arguments !\n";
     exit;
 }
@@ -21,9 +21,9 @@
 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 $output4 = $ARGV[8]; # fichier avec les catégories de COG pour core-génome / génome accessoire / gènes spé
+my $output8 = $ARGV[8];
+my $output9 = $ARGV[9];
 # my $output9 = $ARGV[11];


@@ -136,17 +136,17 @@
             $unique_col_detected = $i;
             my $espece = $hSpecies{$samples[$i]};

-            my @table_genes = split (',', $val);
+            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;
+                $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;
+                $comptage_especes{$espece} .= ','.$genes;

             }
         }
@@ -159,7 +159,7 @@

     if (scalar keys(%comptage_especes) == 1){
         my $list = $comptage_especes{$seule_espece};
-        my @table = split(" ",$list);
+        my @table = split(",",$list);
         foreach my $gene (@table){
             $NonStrict_Spe{$gene} = $seule_espece;
         }
@@ -178,22 +178,24 @@
         # 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 ',')
+            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"}++;
-
+                $coregenes2{$i}{$gene}= $orthogroup;
+                $coregenes{$gene}= $orthogroup;
             }

-            # $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;
+        # my @liste_of_genes = split (',', $gene_random);
+        # my $first_gene = $liste_of_genes[0];
+        # $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;
+        # $coregenes{$first_gene}= $orthogroup;


         if (!$coregene_line){
@@ -209,7 +211,7 @@


         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 @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"}++;
@@ -220,12 +222,12 @@


         }
-    #     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;
+        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 !
@@ -238,7 +240,7 @@
         # }

         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 @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"}++;
@@ -251,7 +253,7 @@
         }


-        my @liste_of_genes = split (',', $gene_random);
+        my @liste_of_genes = split (', ', $gene_random);
         my $first_gene = $liste_of_genes[0];
         $accessorygenes{$first_gene}= $orthogroup;
     }
@@ -276,8 +278,8 @@
 #         print "$strain\t$gene\n";
 #     }
 # }
-# foreach my $gene (keys (%NonStrict_Spe)){
-#     print $NonStrict_Spe{$gene}."\t$gene \n";
+# foreach my $gene (keys (%coregenes)){
+#     print "$gene\n";
 # }
 # exit;

@@ -352,6 +354,7 @@

 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
+my %GrpSpe = ();

 foreach my $i (keys(%Genes_of_OG)){
     foreach my $ortho (keys %{$Genes_of_OG{$i}}){
@@ -360,10 +363,11 @@
         if ($Hash_Specific{$ortho}){
             my $specie = $Hash_Specific{$ortho};

-            my @liste_genes = split(' ',$gene);
+            my @liste_genes = split(',',$gene);
             foreach my $g(@liste_genes){
                 $Gene_Specie_Spe{$g} = $specie;
                 $Species_Spe_Count{$specie}++;
+                $GrpSpe{$i}{$ortho} = $g;

             }

@@ -371,14 +375,18 @@
     }
 }

-
+# foreach my $i (keys(%GrpSpe)) {
+#     foreach my $ortho (keys %{$GrpSpe{$i}}) {
+#         print $GrpSpe{$i}{$ortho}."\n";
+#     }
+# }
 # exit;

 # my @table_keys = ();
 my $nb_groupSpe_genes = 0;

 foreach my $gene (keys (%Gene_Specie_Spe)) {
-    my @table_keys = split (' ', $gene);
+    my @table_keys = split (',', $gene);
     foreach my $unique_gene (@table_keys) {
         $nb_groupSpe_genes++;
     }
@@ -669,26 +677,26 @@

 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 $!;
+open (OUT8, ">$output8") or die $!;

 # my $nb_files = scalar keys @files;


-print OUT7 "\t";
+print OUT8 "\t";

 foreach my $category(@orders){
 # foreach my $category (sort keys (%Acc_Cat_Esp_Count)) {
     # my $cat = $category."\t";
-    print OUT7  $category."\t";
+    print OUT8  $category."\t";
 }

-print OUT7 "\n";
+print OUT8 "\n";


 # foreach my $category (sort keys (%Global_Count)){
 foreach my $type (sort keys (%Global_Count)){

-    print OUT7 "$type\t";
+    print OUT8 "$type\t";
     #foreach my $category (sort keys (%{$Global_Count{$type}})){
     foreach my $category(@orders){

@@ -715,34 +723,34 @@
             # 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 OUT8 "$odds_ratio;";
             }

         }
-        print OUT7 "\t";
+        print OUT8 "\t";
     }

-    print OUT7 "\n";
+    print OUT8 "\n";
 }


-print OUT7 "\n";
-close (OUT7);
+print OUT8 "\n";
+close (OUT8);


 #//////////////////////////////////////////////
-open (OUT8, ">$output8") or die $!;
+open (OUT9, ">$output9") or die $!;


-print OUT8 "\t";
+print OUT9 "\t";

 #
 foreach my $category(@orders){
-    print OUT8  $category."\t";
+    print OUT9  $category."\t";
 }

-print OUT8 "\n";
+print OUT9 "\n";


 # foreach my $category (sort keys (%Global_Count)){
@@ -750,7 +758,7 @@
     # 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";
+    print OUT9 "$specie\t";
     foreach my $category (sort keys (%{$Species_Count{$specie}})){

         foreach my $strain (sort keys (%{$Species_Count{$specie}{$category}})){
@@ -779,18 +787,18 @@
             # 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 OUT9 "$odds_ratio;";
             }

         }
-        print OUT8 "\t";
+        print OUT9 "\t";
     }
-print OUT8 "\n";
+print OUT9 "\n";
 }

 # print OUT9 "\n";

-close (OUT8);
+close (OUT9);

 #///////////////////////////////////////////////////
 # open (OUT9, '>', $output9) or die $!;
@@ -957,57 +965,57 @@
 #             $nbr = $hCount2{$category}{$especes};
 #         }
 # STEP 4 : AFFICHAGE DANS LE FICHIER DE SORTIE ------------------------------------------------------------------------------
-open (OUT4, ">$output4") or die $!;
+# open (OUT4, ">$output4") or die $!;

-print OUT4 "COG categories"."\t"."Core-genome"."\t"."Accessory genome"."\t"."Strain Specific genes"."\n";
+# 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 $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
+#     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 ($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 ($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);
+#         # # 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é
@@ -1129,6 +1137,8 @@
     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");
@@ -1163,13 +1173,15 @@
         # 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);
 }
-
+# exit;
 mkdir("StrainSpecific");
 foreach my $i (keys (%specificgenes2)){

@@ -1211,7 +1223,7 @@


 mkdir("GroupSpecific");
-foreach my $i (keys (%Genes_of_OG)){
+foreach my $i (keys (%GrpSpe)){
     if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas
         next;
     }
@@ -1223,8 +1235,8 @@
     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;
+    my $refGrpSpe = $GrpSpe{$i};
+    my %subhash = %$refGrpSpe;

     foreach my $orthogroup (keys (%subhash)){
         if ($Hash_Specific{$orthogroup} && $Hash_Specific{$orthogroup} eq $specie_name){
@@ -1234,9 +1246,13 @@
             if ($Cog_of_gene{$gene}){
                 $cat = $Cog_of_gene{$gene};
             }
+            # if (!$Gene_position{$gene}){
+            #     print "$gene\n coucou"; exit;
+            # }
             print OUT6 $orthogroup."\t".$subhash{$orthogroup}."\t".$Gene_position{$gene}."\t".$cat."\t".$Hash_Convert{$cat}."\n";
         }

     }
     close (OUT6);
 }
+# exit;
\ No newline at end of file