comparison pangenomeCogAnalysis_V1.pl @ 12:b36506d26a43 draft

Uploaded

11:9ac1a58fda89

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 $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(<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)

            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
        }

        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
        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];

        }
        # 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)) {

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 %hash_of_genes = ();


foreach $g (@list_gff){

        my $end = $table_gff[4];
        my $gene_name = $table_gff[8];
        my $type = $table_gff[2];



        if ($type && ($type eq "mRNA" or $type eq "CDS") && $gene_name =~ /ID=([^;]+);/){
            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;  
                }
            }    
            
            $Gene_position{$gene}="$chr\t$start\t$end";
        }   

        # foreach my $gene (keys (%hash_of_genes)){
        #     my $orthogrp = $hGene_OG{$gene};

    }

    close (OUT5);
}


mkdir("GroupSpecific");
foreach my $i (keys (%Genes_of_OG)){
    if (!$hCol_Annotated{$i}) { # si le fichier GFF n'existe pas 
        next;

12:b36506d26a43

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


while(<M>) {

    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;


    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 @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
        }

        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é
        
        # print "$gene_random\n";
        # print "$line\n";
        # print "$unique_col_detected\n";

        

        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 !
        my @list_of_genes = split (',', $gene_random);
        my $first_gene = $list_of_genes[0];

        }
        # 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(%strain_specie)) {
#     print "i=$k strain=$v\n";
# }
# exit;
# foreach my $oups (keys (%coregenes)) {

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){

        my $end = $table_gff[4];
        my $gene_name = $table_gff[8];
        my $type = $table_gff[2];


#or $type eq "CDS"
        if ($type && $type eq "mRNA" && $gene_name =~ /ID=([^;]+);/){
            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};

    }

    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;