# HG changeset patch # User mgarnier # Date 1629400340 0 # Node ID 574fece473bf3693485ab5260a3b71858eced6dd # Parent 3b12b1b1d67f8e4a39d611f90c231d17ccfb1c20 Uploaded diff -r 3b12b1b1d67f -r 574fece473bf pangenomeCogAnalysis.pl --- 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