# HG changeset patch
# User mgarnier
# Date 1625059819 0
# Node ID 731fb6cb324be8a8a8a1c2b67296b1e1e00be2d0
Uploaded
diff -r 000000000000 -r 731fb6cb324b pangenomeCogAnalysis.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/pangenomeCogAnalysis.xml Wed Jun 30 13:30:19 2021 +0000
@@ -0,0 +1,40 @@
+
+ Provide statistics from pangenome matrix and COG annotations
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff -r 000000000000 -r 731fb6cb324b pangenomeCogAnalysis_V1.pl
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/pangenomeCogAnalysis_V1.pl Wed Jun 30 13:30:19 2021 +0000
@@ -0,0 +1,491 @@
+#!/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[6]; # liste des espèces avec leurs orthogroupes (présence-absence)
+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 %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 = ){
+
+ $line =~s/\n//g; $line =~s/\r//g;
+ my @sp = split('\t', $line);
+ # print "$line\n";
+ # exit;
+ $hSpecies{$sp[0]} = $sp[1]; # 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 = ;
+$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: N_Id ; 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)
+
+
+while() {
+
+ my $nb_found = 0;
+ my @infos = split(/\t/,$_);
+ 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;
+
+ for (my $i=1; $i <= $#infos; $i++){ # on travaille par ligne puis dans chaque ligne (while()), 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)
+ }
+
+ else { # si jamais il n'y a rien dans la cellule...
+ $combi_abs .= "_".$i; # ... on fait la même chose mais avec $combi_abs
+ }
+
+ }
+
+ # $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";
+ 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";
+ # }
+
+ }
+ 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
+ my $first_gene = $list_of_genes[0]; # on ne prend que le premier
+ $specificgenes{$first_gene}= $orthogroup; # et pareil on crée la table de hash
+ }
+
+ else { # là c'est le génome accessoire, i.e tout le reste !
+ my @list_of_genes = split (',', $gene_random);
+ my $first_gene = $list_of_genes[0];
+ $accessorygenes{$first_gene}= $orthogroup;
+ }
+
+}
+# while (my ($k,$v) = each(%coregenes)) {
+# print "gene=$k OG=$v\n";
+# }
+# exit;
+# foreach my $oups (keys (%coregenes)) {
+# print "$oups\n";
+# }
+# exit;
+
+close (M);
+
+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";
+ }
+
+ if ($ortho_absents){
+ # open (OUT2,">$species.$combination.absents.list.txt");
+ print OUT "> $species - absent\n";
+ print OUT "$ortho_absents\n";
+ }
+
+# close(OUT2);
+}
+
+close(OUT);
+
+#//////////////////////////////////////////////////////////////////////////////////////////////////
+
+############################################### 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 ++++++++++++ #
+foreach $f (@files){
+ $hCorrespondance{$f} = $list[$l++]; # on fait correspondre pour chaque fichier de COG, un nom de souche
+}
+
+# #Affichage du hash
+# foreach $f (keys %hCorrespondance){
+# print $f."=>".$hCorrespondance{$f}."\n"
+# }
+
+# 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 %Cog_of_gene = ();
+
+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
+
+ while (my $line2 = ){
+
+ $line2 =~s/\n//g; $line2 =~s/\r//g; # on procède ligne par ligne
+ my @Genes = split('\t', $line2);
+ my $gene = $Genes[0];
+ my $first_cat = $Genes[2];
+ $Cog_of_gene{$gene} = $first_cat;
+
+ 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 ($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);
+
+ # 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;
+ }
+
+}
+# foreach my $category (sort keys (%hSpecific_Cat_Esp)) { # parcours au niveau de la 1ere clé
+
+# foreach my $especeee (keys %{$hSpecific_Cat_Esp{$category} }) { # parcours au niveau de la 2e clé pour la $category donnée
+
+# print "$category\t$especeee\t$hSpecific_Cat_Esp{$category}{$especeee}\n"; # on crée une sortie qui affiche en somme notre hash %hCount2
+# }
+# }
+# exit;
+
+# STEP 4 : AFFICHAGE DANS LE FICHIER DE SORTIE ------------------------------------------------------------------------------
+open (OUT4, ">$output4") or die $!;
+
+print OUT4 "Species"."\t"."COG categories"."\t"."Core-genome"."\t"."Accessory genome"."\t"."Specific genes"."\n";
+
+foreach my $category (sort keys (%hCount2)){ # parcours de la table %hCount2 au niveau des catégories
+ 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 @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);
+
+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 %hCore_Cat = reverse (%hCore_Cat);
+
+# ++++++++++++ parcours de 2 listes en même temps ++++++++++++ #
+foreach $g (@list_gff){
+ $hgff_order{$g} = $order_gff[$o++]; # on fait correspondre pour chaque fichier GFF, un nom de souche
+ open (G, $g);
+ while () {
+ 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=([^;]+);/){
+ my $gene = $1;
+ # print $gene."\n";
+ # exit;
+ foreach my $cog (keys (%hCore_Cat)){
+ if ($hCore_Cat{$cog} eq $gene){
+ $Cat_genes{$gene}=$cog;
+ }
+ }
+
+ $Gene_position{$gene}="$chr\t$start\t$end";
+ }
+
+
+ }
+
+ close (G);
+}
+
+open (OUT5, "> $output5");
+print OUT5 "Orthogroups\tGenes\tChromosomes\tStart\tEnd\tCOG categories\n";
+
+
+
+
+# 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;
+
+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"."\t".$Gene_position{$gene}.$cat."\n";
+
+}
+
+close (OUT5);