changeset 2:0428ce25da81 draft

Uploaded

--- a/pangenomeCogAnalysis_V1.pl	Wed Jun 30 15:43:43 2021 +0000
+++ b/pangenomeCogAnalysis_V1.pl	Fri Jul 02 14:53:33 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;
 }
@@ -22,12 +22,12 @@
 my $output2 = $ARGV[7]; # fichier des moyennes
 my $output3 = $ARGV[8]; # fichier de la liste des valeurs pour chaque catégorie de COG et pour chaque espèce 
 my $output4 = $ARGV[9]; # fichier avec les catégories de COG pour core-génome / génome accessoire / gènes spé
-my $output5 = $ARGV[10]; # sortie qui affiche les GFF
 
 
 # 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 ###########################
@@ -38,7 +38,7 @@
     my @sp = split('\t', $line);
     # print "$line\n"; 
     # exit; 
-    $hSpecies{$sp[0]} = $sp[1]; # key = N_Id ; val = name
+    $hSpecies{$sp[0]} = $sp[1]; # HASH -> key: N_Id ; val: name
 
 }
  my $nbr = keys (%hSpecies); #compter le nombre de souches max 
@@ -91,11 +91,18 @@
 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 %Genes_of_OG = (); # HASH -> key1: orthogroupe ; key2: colonne i ; val: gène
+
 
 while(<M>) {
 
+    my $line = $_;
+    $line =~s/\n//g; $line =~s/\r//g;
     my $nb_found = 0;
-	my @infos = split(/\t/,$_);
+	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 = "";
@@ -103,6 +110,7 @@
     my $val;
 	my $gene_random;
 
+
     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
@@ -111,6 +119,10 @@
             $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)
+
+            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
         }
 
         else { # si jamais il n'y a rien dans la cellule...
@@ -124,19 +136,23 @@
     $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";
-        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 ',')
-        my $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement !
-        $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)
-        # foreach my $oups (@list_of_genes) {
-        #     print "$orthogroup\n";
-        #     print "$first_gene\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 $first_gene = $list_of_genes[0]; # prend la valeur du premier gène uniquement !
+            $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;
+            
+            
+            
+            
+        }
+        if (!$coregene_line){
+            $coregene_line = $line;
+        } 
     }
     elsif ($nb_found == 1) { # si on a un gène spé
         my @list_of_genes = split (',', $gene_random); # idem, on ne veut qu'un seul gène donc on crée la liste
@@ -151,17 +167,47 @@
     } 
        
 }
-# while (my ($k,$v) = each(%coregenes)) {
-#     print "gene=$k OG=$v\n";
+
+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 (%coregenes2)){ # parcours de la table %hCount2 au niveau des catégories
+#     foreach my $gene (keys %{$coregenes2{$i} }){ # parcours de la table %hCount2 au niveau des espèces 
+#         print "$i\t$gene\t".$coregenes2{$i}{$gene}."\n";
+#     }
+# }        
+
+# while (my ($k,$v) = each(%strain_specie)) {
+#     print "i=$k strain=$v\n";
 # }
 # exit;
 # foreach my $oups (keys (%coregenes)) {
 #     print "$oups\n";
 # }
-#     exit;
+    # exit;
 
 close (M);
 
+my %Hash_Specific = ();
+
 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)
@@ -174,6 +220,10 @@
     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){
@@ -241,7 +291,8 @@
 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 %Cog_of_gene = ();
+my %Cog_of_gene = (); # HASH -> key: gène ; val: cat de COG
+my %Specie_of_gene = (); # HASH -> key: gène ; val: souche
 
 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
@@ -417,7 +468,7 @@
 
 ############################################### GFF ###############################################
 
-my @list_gff = split(',', $annotation_GFF); # liste des différents fichiers GFF (qui se retrouvent dans le dossier Annotation Maker)
+
 my @order_gff = split(',', $order_GFF); # liste de l'ordre des souches 
 my ($g,$o);
 
@@ -425,10 +476,11 @@
 my %Gene_position = ();
 my %Cat_genes = ();
 
-# my %hCore_Cat = reverse (%hCore_Cat);
+my %hash_of_genes = ();
 
 # ++++++++++++ parcours de 2 listes en même temps ++++++++++++ # 
 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>) {
@@ -445,6 +497,8 @@
             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;  
@@ -454,38 +508,82 @@
             $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);
 }
 
-open (OUT5, "> $output5");
-print OUT5 "Orthogroups\tGenes\tChromosomes\tStart\tEnd\tCOG categories\n";
+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};
 
 
 
-# foreach my $gene1 (keys (%Gene_position)){
-#     if ($hCore_Cat{$gene1}){
-#         # print "$hCore_Cat{$gene1}\n";
-#         my $cog = $hCore_Cat{$gene1};
-#         $Hash_genes{$gene1}=$cog;
-#     } 
-# }
-# while (my ($k,$v) = each(%Hash_genes)) {
-#             print "gene=$k cog=$v\n";
-# }
-# exit;
+    open (OUT5, "> Core/$strain_name.$specie_name.txt") or die "Cannot create file $!\n";
+    print OUT5 "Orthogroup\tGene\tChromosome\tStart\tEnd\tCOG categories\n";
 
-foreach my $gene (keys (%coregenes)){
-# print "$gene\n";
-    my $cat = "unknown";
-    if ($Cog_of_gene{$gene}){
-        $cat = $Cog_of_gene{$gene};
-    }     
-    print OUT5 $coregenes{$gene}."\t"."$gene\t".$Gene_position{$gene}."\t".$cat."\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."\n";
+       
+    }
+
+    close (OUT5);
 }
 
-close (OUT5);
+
+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\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."\n";
+        }
+        
+    }
+    close (OUT6);
+}