Commit 14f7e449 authored by Sebastien's avatar Sebastien

Selection de la proteine representant le groupe a partir matrice de distance

parent 76a7c1ef
......@@ -7,7 +7,7 @@ RUN apt-get -q update && \
hardinfo wget patch vim gawk unzip apache2 zlib1g-dev libpng-dev libgd2-noxpm-dev build-essential && \
apt-get autoremove -y && apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
RUN apt-get -q update && apt-get install -yq cron at libtree-simple-perl libbio-perl-perl ncbi-blast+ sqlite liberror-perl libswitch-perl libxml-simple-perl libjson-perl libapache-session-perl libnet-ldap-perl libapache-htpasswd-perl ssmtp fasttree mafft biosquid python2.7 python-scipy python-sklearn python-numpy python-fastcluster mcl libarray-utils-perl liburi-encode-perl cd-hit
RUN apt-get -q update && apt-get install -yq cron at libtree-simple-perl libbio-perl-perl ncbi-blast+ sqlite liberror-perl libswitch-perl libxml-simple-perl libjson-perl libapache-session-perl libnet-ldap-perl libapache-htpasswd-perl ssmtp fasttree mafft biosquid python2.7 python-scipy python-sklearn python-numpy python-fastcluster mcl libarray-utils-perl liburi-encode-perl embassy-phylip
WORKDIR /var/www/
......
......@@ -27,7 +27,8 @@ sub CreatePanProteome
foreach my $code (keys %{$rh_analysis_params->{'proteomes'}})
{
if ($rh_analysis_params->{'proteomes'}->{$code}->{'reference'} eq 'true')
if ( $rh_analysis_params->{'proteomes'}->{$code}->{'reference'} eq 'true'
or $rh_analysis_params->{'proteomes'}->{$code}->{'reference'} == &TRUE)
{
&Debug("INFO - $code is a reference proteome");
push(@a_references, $code);
......@@ -38,36 +39,18 @@ sub CreatePanProteome
foreach my $groupid (keys %{$rh_homologygroups})
{
$refexist = 0;
@a_allseqs1 = ();
foreach my $code (@a_references)
foreach my $code (keys %{$rh_homologygroups->{$groupid}})
{
if (exists($rh_homologygroups->{$groupid}->{$code}))
foreach my $seqid (@{$rh_homologygroups->{$groupid}->{$code}})
{
$refexist = 1;
foreach my $seqid (@{$rh_homologygroups->{$groupid}->{$code}})
{
@a_seqs1 = ();
push(@a_seqs1, $code, $seqid);
push(@a_allseqs1, [@a_seqs1]);
}
@a_seqs1 = ();
push(@a_seqs1, $code, $seqid);
push(@a_allseqs1, [@a_seqs1]);
}
}
if ($refexist == 0)
{
foreach my $code (keys %{$rh_homologygroups->{$groupid}})
{
foreach my $seqid (@{$rh_homologygroups->{$groupid}->{$code}})
{
@a_seqs1 = ();
push(@a_seqs1, $code, $seqid);
push(@a_allseqs1, [@a_seqs1]);
}
}
}
$h_list_group_select_seqid_code{$groupid} = &SelectSequenceRepGroupe($rh_species, \@a_allseqs1);
$h_list_group_select_seqid_code{$groupid} = &SelectSequenceRepGroupe($rh_species, \@a_allseqs1, \@a_references);
}
&PrintPanProteome($outdir, $rh_specific, \%h_list_group_select_seqid_code,
......@@ -90,40 +73,62 @@ sub CreatePanProteome
sub SelectSequenceRepGroupe
{
my ($rh_species, $ra_species_seqid_list) = @_;
my ($rh_species, $ra_species_seqid_list, $ra_reference) = @_;
my $selecSeqid = "";
my $selecCode = "";
my ($fh_tmp, $tmpfile) = File::Temp::tempfile(UNLINK => 1);
foreach my $ra_species_seqid (@{$ra_species_seqid_list})
my $id = 0;
my ($fh_tmp, $tmpfile) = File::Temp::tempfile(UNLINK => 1);
foreach my $ra_species_seqid (@{$ra_species_seqid_list})
{
my $species_code = $ra_species_seqid->[0];
my $seqid = $ra_species_seqid->[1];
my $sequence = $rh_species->{$species_code}->{'sequences'}->{$seqid}->{'sequence'};
print $fh_tmp <<END;
>$species_code|$seqid acc=$seqid sp=$species_code
my $species_code = $ra_species_seqid->[0];
my $seqid = $ra_species_seqid->[1];
my $sequence = $rh_species->{$species_code}->{'sequences'}->{$seqid}->{'sequence'};
print $fh_tmp <<END;
>$id $species_code|$seqid acc=$seqid sp=$species_code
$sequence
END
}
$fh_tmp->close();
my $cdhit_cmd = "cdhit -i $tmpfile -o $tmpfile.cdhit -c 0.5 -n 3";
&System("($cdhit_cmd > /dev/null ) 2> /dev/null");
my %h_cdhit = ();
&FastaToHash("$tmpfile.cdhit",\%h_cdhit, 'noseq');
my $longest = 0;
my $representative;
foreach my $cdhit_seq (keys %h_cdhit)
{
my $len = $h_cdhit{$cdhit_seq}->{len};
if ($len > $longest)
{
$longest = $len;
$representative = $cdhit_seq;
}
}
return [$h_cdhit{$representative}->{'acc'}, $h_cdhit{$representative}->{'sp'}];
$id++;
}
$fh_tmp->close();
#Calcul de la matrice de distance sur alignement du groupe
my $distance_cmd = "mafft --auto --quiet $tmpfile | fprotdist -sequence stdin -out stdout -auto -noprogress";
&Log("CMD - $distance_cmd");
my @a_distance_lines = `$distance_cmd`;
chomp @a_distance_lines;
#Calcul de la somme par ligne de la matrice de distance (score distance globale de chaque proteine avec les autres)
foreach my $line (@a_distance_lines)
{
if ($line =~ /^\d+\s+\d\.\d+/)
{
my ($id, @a_distances) = split(/\s+/, $line);
my $score = eval join '+', @a_distances;
push(@{$ra_species_seqid_list->[$id]}, $score);
}
}
#Tri sur ce score: la proteine avec le plus petit score (1ere du tri) ou la premiere protein issue d'une proteome de reference est choisie
my @a_representative = ();
foreach my $ra_species_seqid (sort {$a->[2] <=> $b->[2]} @{$ra_species_seqid_list})
{
my $species_code = $ra_species_seqid->[0];
my $seqid = $ra_species_seqid->[1];
@a_representative = ($seqid, $species_code) if (scalar @a_representative == 0);
my $score = $ra_species_seqid->[2];
if (grep (/^$species_code$/, @{$ra_reference}))
{
@a_representative = ($seqid, $species_code);
last;
}
}
return \@a_representative;
}
=head3 function PrintPanProteome
......@@ -143,10 +148,9 @@ sub PrintPanProteome
my ($outdir, $rh_specific, $rh_h_list_group_select_seqid_code, $rh_species, $rh_nb_repre_group_by_proteome) = @_;
my $subdirectory = "pan_proteome";
my $subdirectory_description = "Pan proteome analysis";
my $rh_description = {$subdirectory => $subdirectory_description};
my $rh_description = {$subdirectory => $subdirectory_description};
mkdir "$outdir/$subdirectory";
$rh_description->{"$subdirectory/PanProteome.fasta"} = 'Pan proteome multifasta file';
my $fh_panproteome = &GetStreamOut("$outdir/$subdirectory/PanProteome.fasta");
......@@ -172,7 +176,7 @@ sub PrintPanProteome
}
$fh_panproteome->close;
$rh_description->{"$subdirectory/PanProteomeDescription.xls"} = 'Pan proteome description table';
$rh_description->{"$subdirectory/PanProteomeDescription.xls"} = 'Pan proteome description table';
my $fh_panproteome_desc = &GetStreamOut("$outdir/$subdirectory/PanProteomeDescription.xls");
print $fh_panproteome_desc
"#Species\tNb seq in panproteome\tNb seq rep group in panproteome\tNb seq uniq in panproteome\n";
......@@ -200,8 +204,8 @@ For each species
Add specific singlecopy protein
END
$fh_readme->close;
push(@A_OUTPUT_DATA, {'directory' => $subdirectory, 'description' => $rh_description});
push(@A_OUTPUT_DATA, {'directory' => $subdirectory, 'description' => $rh_description});
return;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment