Raw content of Bio::EnsEMBL::Compara::RunnableDB::MCoffee
<![CDATA[#
# You may distribute this module under the same terms as perl itself
#
# POD documentation - main docs before the code
=pod
=head1 NAME
Bio::EnsEMBL::Compara::RunnableDB::MCoffee
=cut
=head1 SYNOPSIS
my $db = Bio::EnsEMBL::Compara::DBAdaptor->new($locator);
my $mcoffee = Bio::EnsEMBL::Compara::RunnableDB::Mcoffee->new (
-db => $db,
-input_id => $input_id,
-analysis => $analysis );
$mcoffee->fetch_input(); #reads from DB
$mcoffee->run();
$mcoffee->output();
$mcoffee->write_output(); #writes to DB
=cut
=head1 DESCRIPTION
This Analysis/RunnableDB is designed to take a protein_tree cluster as input
Run an MCOFFEE multiple alignment on it, and store the resulting alignment
back into the protein_tree_member and protein_tree_member_score table.
input_id/parameters format eg: "{'protein_tree_id'=>726093, 'clusterset_id'=>1}"
protein_tree_id : use family_id to run multiple alignment on its members
options : commandline options to pass to the 'mcoffee' program
=cut
=head1 CONTACT
Contact Albert Vilella on module implemetation/design detail: avilella@ebi.ac.uk
Contact Javier Herrero on EnsEMBL/Compara: jherrero@ebi.ac.uk
Contact Ewan Birney on EnsEMBL in general: birney@sanger.ac.uk
=cut
=head1 APPENDIX
The rest of the documentation details each of the object methods.
Internal methods are usually preceded with a _
=cut
package Bio::EnsEMBL::Compara::RunnableDB::MCoffee;
use strict;
use Getopt::Long;
use IO::File;
use File::Basename;
use File::Path;
use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::BaseAlignFeature;
use Bio::EnsEMBL::Compara::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Compara::DBSQL::PeptideAlignFeatureAdaptor;
use Bio::EnsEMBL::Compara::Member;
use Time::HiRes qw(time gettimeofday tv_interval);
# use POSIX qw(ceil floor);
use Bio::EnsEMBL::Hive;
our @ISA = qw(Bio::EnsEMBL::Hive::Process);
=head2 fetch_input
Title : fetch_input
Usage : $self->fetch_input
Function: Fetches input data for mcoffee from the database
Returns : none
Args : none
=cut
sub fetch_input {
my( $self) = @_;
$self->{'comparaDBA'} = Bio::EnsEMBL::Compara::DBSQL::DBAdaptor->new(-DBCONN=>$self->db->dbc);
### DEFAULT PARAMETERS ###
my $p = '';
$p = 'use_exon_boundaries';
$self->{$p} = 0;
$p = 'method';
$self->{$p} = 'fmcoffee';
$p = 'output_table';
$self->{$p} = 'protein_tree_member';
$p = 'max_gene_count';
$self->{$p} = 1500;
$p = 'options';
$self->{$p} = '';
#########################
# Fetch parameters from the two possible locations. Input_id takes precedence!
$self->get_params($self->parameters);
$self->get_params($self->input_id);
$self->check_if_exit_cleanly;
$self->{treeDBA} = $self->{'comparaDBA'}->get_ProteinTreeAdaptor;
#
# A little logic, depending on the input params.
#
# Protein Tree input.
if (defined $self->{'protein_tree_id'}) {
my $id = $self->{'protein_tree_id'};
$self->{'protein_tree'} = $self->{'comparaDBA'}->get_ProteinTreeAdaptor->fetch_node_by_node_id($id);
$self->{'protein_tree'}->flatten_tree; # This makes retries safer
$self->{'input_fasta'} = $self->dumpProteinTreeToWorkdir($self->{'protein_tree'});
if ($self->{'use_exon_boundaries'}) {
$self->{'input_fasta_exons'} = $self->dumpProteinTreeToWorkDir($self->{'protein_tree'},1); # The extra "1" at the end adds the exon markers.
}
}
if (defined($self->{'redo'}) && $self->{'method'} eq 'unalign') {
# Redo - take previously existing alignment - post-process it
$self->{redo_sa} = $self->{'protein_tree'}->get_SimpleAlign(-id_type => 'MEMBER');
$self->{redo_sa}->set_displayname_flat(1);
$self->{redo_alnname} = $self->worker_temp_directory . $self->{protein_tree}->node_id . ".fasta";
my $alignout = Bio::AlignIO->new(-file => ">".$self->{redo_alnname},
-format => "fasta");
$alignout->write_aln($self->{redo_sa});
}
# Auto-switch to fmcoffee if gene count is too big.
if ($self->{'method'} eq 'cmcoffee') {
if ($self->{'protein_tree'}->get_tagvalue('gene_count') > 100) {
$self->{'method'} = 'mafft';
# $self->{'method'} = 'fmcoffee';
print "MCoffee, auto-switch method to mafft because gene count > 100 \n";
}
}
print "RETRY COUNT:".$self->input_job->retry_count()."\n";
# if ($self->input_job->retry_count >= 1) {
# if ($self->{'protein_tree'}->get_tagvalue('gene_count') > 200) {
# $self->{'method'} = 'mafft';
# # # HIGHMEM ensembl-hive code still experimental
# # unless (defined($self->worker->{HIGHMEM})) {
# # $self->input_job->update_status('HIGHMEM');
# # $self->DESTROY;
# # throw("Mcoffee job too big: try something else and FAIL it");
# # }
# }
# }
# Auto-switch to fmcoffee on two failures.
if ($self->input_job->retry_count >= 2) {
$self->{'method'} = 'fmcoffee';
}
# Auto-switch to muscle on a third failure.
if ($self->input_job->retry_count >= 3) {
$self->{'method'} = 'mafft';
# actually, we are going to run mafft directly here, not through mcoffee
# maybe in the future we want to use this option in tcoffee:
# t_coffee ..... -dp_mode myers_miller_pair_wise
}
print "MCoffee alignment method: ".$self->{'method'}."\n";
#
# Ways to fail the job before running.
#
# No input specified.
if (!defined($self->{'protein_tree'})) {
$self->DESTROY;
throw("MCoffee job no input protein_tree");
}
# Error writing input Fasta file.
if (!$self->{'input_fasta'}) {
$self->DESTROY;
throw("MCoffee job, error writing input Fasta");
}
# Gene count too big.
if ($self->{'protein_tree'}->get_tagvalue('gene_count') > $self->{'max_gene_count'}) {
$self->dataflow_output_id($self->input_id, 2);
$self->input_job->update_status('FAILED');
$self->DESTROY;
throw("Mcoffee job too big: try something else and FAIL it");
}
# Retry count >= 3.
if ($self->input_job->retry_count >= 5) {
$self->dataflow_output_id($self->input_id, 2);
$self->input_job->update_status('FAILED');
$self->DESTROY;
throw("Mcoffee job failed >=3 times: try something else and FAIL it");
}
return 1;
}
=head2 run
Title : run
Usage : $self->run
Function: runs MCOFFEE
Returns : none
Args : none
=cut
sub run
{
my $self = shift;
$self->check_if_exit_cleanly;
$self->{'mcoffee_starttime'} = time()*1000;
$self->run_mcoffee;
}
=head2 write_output
Title : write_output
Usage : $self->write_output
Function: parse mcoffee output and update protein_tree_member tables
Returns : none
Args : none
=cut
sub write_output {
my $self = shift;
$self->check_if_exit_cleanly;
$self->parse_and_store_alignment_into_proteintree;
#
# Store various alignment tags.
#
$self->_store_aln_tags unless ($self->{redo});
}
sub DESTROY {
my $self = shift;
if($self->{'protein_tree'}) {
$self->{'protein_tree'}->release_tree;
$self->{'protein_tree'} = undef;
}
# Cleanup temp files and stuff.
unlink ($self->{'input_params'}) if($self->{'input_params'});
unlink ($self->{'input_fasta'}) if($self->{'input_fasta'});
if($self->{'mcoffee_output'}) {
unlink ($self->{'mcoffee_output'});
unlink ($self->{'mcoffee_output'} . ".log");
}
$self->SUPER::DESTROY if $self->can("SUPER::DESTROY");
}
##########################################
#
# internal methods
#
##########################################
sub get_params {
my $self = shift;
my $param_string = shift;
return unless($param_string);
print("parsing parameter string : ",$param_string,"\n") if($self->debug);
my $params = eval($param_string);
return unless($params);
if($self->debug) {
foreach my $key (keys %$params) {
print(" $key\t=>\t", $params->{$key}, "\n");
}
}
my $p;
# METHOD: The style of MCoffee to be run for this alignment.
# Could be: fmcoffee or cmcoffee
$p = 'method';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
# cutoff: for filtering
$p = 'cutoff';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
# OUTPUT_TABLE: self-explanatory.
$p = 'output_table';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
# Loads the protein tree if we have a protein_tree_id
$p = 'protein_tree_id';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
# clusterset_id
$p = 'clusterset_id';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
# Extra command-line options.
$p = 'options';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
# Set a limit on the number of members to align.
$p = 'max_gene_count';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
$p = 'use_exon_boundaries';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
# This is analysis, not production: 'redo' e.g. '1:1000000' from clusterset_id 1 to a different clusterset_id 10000000
$p = 'redo';
$self->{$p} = $params->{$p} if (defined($params->{$p}));
return;
}
sub run_mcoffee
{
my $self = shift;
my $input_fasta = $self->{'input_fasta'};
my $mcoffee_output = $self->worker_temp_directory . "output.mfa";
$mcoffee_output =~ s/\/\//\//g; $self->{'mcoffee_output'} = $mcoffee_output;
# (Note: t_coffee automatically uses the .mfa output as the basename for the score output)
my $mcoffee_scores = $mcoffee_output . ".score_ascii";
$mcoffee_scores =~ s/\/\//\//g; $self->{'mcoffee_scores'} = $mcoffee_scores;
my $tree_temp = $self->worker_temp_directory . "tree_temp.dnd";
$tree_temp =~ s/\/\//\//g; $self->{'mcoffee_tree'} = $tree_temp;
my $mcoffee_executable = $self->analysis->program_file;
unless (-e $mcoffee_executable) {
print "Using default T-Coffee executable!\n";
$mcoffee_executable = "/nfs/acari/avilella/src/tcoffee/compara/t_coffee";
}
throw("can't find a M-Coffee executable to run\n") unless(-e $mcoffee_executable);
#
# Make the t_coffee temp dir.
#
my $tempdir = $self->worker_temp_directory;
print "TEMP DIR: $tempdir\n" if ($self->debug);
#
# Output the params file.
#
$self->{'input_params'} = $self->worker_temp_directory. "temp.params";
my $paramsfile = $self->{'input_params'};
$paramsfile =~ s/\/\//\//g; # converts any // in path to /
open(OUTPARAMS, ">$paramsfile")
or throw("Error opening $paramsfile for write");
my $method_string = '-method=';
if ($self->{'method'} && $self->{'method'} eq 'cmcoffee') {
# CMCoffee, slow, comprehensive multiple alignments.
$method_string .= "mafftgins_msa, muscle_msa, kalign_msa, probcons_msa"; #, t_coffee_msa";
} elsif ($self->{'method'} eq 'fmcoffee') {
# FMCoffee, fast but accurate alignments.
$method_string .= "mafft_msa, muscle_msa, clustalw_msa, kalign_msa";
} elsif ($self->{'method'} eq 'muscle') {
# MUSCLE: quick, kind of crappy alignments.
$method_string .= "muscle_msa";
} elsif ($self->{'method'} eq 'mafft') {
# MAFFT FAST: very quick alignments.
$method_string .= "mafft_msa";
} elsif ($self->{'method'} eq 'prank') {
# PRANK: phylogeny-aware alignment.
$method_string .= "prank_msa";
} elsif (defined($self->{'redo'}) && $self->{'method'} eq 'unalign') {
my $cutoff = $self->{'cutoff'} || 2;
# Unalign module
$method_string = " -other_pg seq_reformat -in " . $self->{redo_alnname} ." -action +aln2overaln lower $cutoff 1>$mcoffee_output";
} else {
throw ("Improper method parameter: ".$self->{'method'});
}
if ($self->{'use_exon_boundaries'}) {
$method_string .= ", exon_pair";
my $exon_file = $self->{'input_fasta_exons'};
print OUTPARAMS "-template_file=$exon_file\n";
}
$method_string .= "\n";
print OUTPARAMS $method_string;
print OUTPARAMS "-mode=mcoffee\n";
print OUTPARAMS "-output=fasta_aln,score_ascii\n";
print OUTPARAMS "-outfile=$mcoffee_output\n";
print OUTPARAMS "-newtree=$tree_temp\n";
close(OUTPARAMS);
my $t_env_filename = $tempdir . "t_coffee_env";
open(TCOFFEE_ENV, ">$t_env_filename")
or throw("Error opening $t_env_filename for write");
print TCOFFEE_ENV "http_proxy_4_TCOFFEE=\n";
print TCOFFEE_ENV "EMAIL_4_TCOFFEE=cedric.notredame\@europe.com\n";
close TCOFFEE_ENV;
# Commandline
my $cmd = $mcoffee_executable;
$cmd .= " ".$input_fasta unless ($self->{redo});
$cmd .= " ". $self->{'options'};
if (defined($self->{'redo'}) && $self->{'method'} eq 'unalign') {
$self->{'mcoffee_scores'} = undef; #these wont have scores
$cmd .= " ". $method_string;
} else {
$cmd .= " -parameters=$paramsfile";
}
print("$cmd\n") if($self->debug);
#
# Output some environment variables for tcoffee
#
my $prefix = "export HOME_4_TCOFFEE=\"$tempdir\";";
$prefix .= "export DIR_4_TCOFFEE=\"$tempdir\";";
$prefix .= "export TMP_4_TCOFFEE=\"$tempdir\";";
$prefix .= "export CACHE_4_TCOFFEE=\"$tempdir\";";
$prefix .= "export NO_ERROR_REPORT_4_TCOFFEE=1;";
$prefix .= "export MAFFT_BINARIES=/nfs/acari/gj1/bin/mafft-bins/binaries;";
print $prefix.$cmd."\n" if ($self->debug);
#
# Run the command.
#
$self->{'comparaDBA'}->dbc->disconnect_when_inactive(1);
my $rc;
if ($self->{'method'} eq 'mafft') {
my $mafft_executable = "/nfs/acari/avilella/src/mafft/mafft-6.522/scripts/mafft";
BEGIN {$ENV{MAFFT_BINARIES} = '/nfs/acari/avilella/src/mafft/mafft-6.522/binaries'; }
print STDERR "### $mafft_executable --auto $input_fasta > $mcoffee_output\n";
$rc = system("$mafft_executable --auto $input_fasta > $mcoffee_output");
$self->{'mcoffee_scores'} = undef; #these wont have scores
} else {
$DB::single=1;1;
$rc = system($prefix.$cmd);
}
$self->{'comparaDBA'}->dbc->disconnect_when_inactive(0);
unless($rc == 0) {
$self->DESTROY;
throw("MCoffee job, error running executable: $\n");
}
}
########################################################
#
# ProteinTree input/output section
#
########################################################
sub update_single_peptide_tree
{
my $self = shift;
my $tree = shift;
foreach my $member (@{$tree->get_all_leaves}) {
next unless($member->isa('Bio::EnsEMBL::Compara::AlignedMember'));
next unless($member->sequence);
$member->cigar_line(length($member->sequence)."M");
$self->{'comparaDBA'}->get_ProteinTreeAdaptor->store($member);
printf("single_pepide_tree %s : %s\n", $member->stable_id, $member->cigar_line) if($self->debug);
}
}
sub dumpProteinTreeToWorkdir
{
my $self = shift;
my $tree = shift;
my $use_exon_boundaries = shift;
my $fastafile;
if ($use_exon_boundaries) {
$fastafile = $self->worker_temp_directory. "proteintree_exon_". $tree->node_id. ".fasta";
} else {
my $node_id = $tree->node_id;
$fastafile = $self->worker_temp_directory. "proteintree_". $node_id. ".fasta";
}
$fastafile =~ s/\/\//\//g; # converts any // in path to /
return $fastafile if(-e $fastafile);
print("fastafile = '$fastafile'\n") if ($self->debug);
open(OUTSEQ, ">$fastafile")
or $self->throw("Error opening $fastafile for write!");
my $seq_id_hash = {};
my $residues = 0;
my $member_list = $tree->get_all_leaves;
$DB::single=1;1;
$self->{'tag_gene_count'} = scalar(@{$member_list});
foreach my $member (@{$member_list}) {
# Double-check we are only using longest
my $longest_member = $member->gene_member->get_longest_peptide_Member;
unless ($longest_member->member_id eq $member->member_id) {
$member->disavow_parent;
$self->{treeDBA}->delete_flattened_leaf($member);
my $updated_gene_count = $tree->get_tagvalue('gene_count');
$updated_gene_count--;
$tree->adaptor->delete_tag($tree->node_id,'gene_count');
$tree->store_tag('gene_count', $updated_gene_count);
next;
}
####
return undef unless ($member->isa("Bio::EnsEMBL::Compara::AlignedMember"));
next if($seq_id_hash->{$member->sequence_id});
$seq_id_hash->{$member->sequence_id} = 1;
my $seq = '';
if ($use_exon_boundaries) {
$seq = $member->get_exon_bounded_sequence;
} else {
$seq = $member->sequence;
}
$residues += $member->seq_length;
$seq =~ s/(.{72})/$1\n/g;
chomp $seq;
print OUTSEQ ">". $member->sequence_id. "\n$seq\n";
}
close OUTSEQ;
if(scalar keys (%{$seq_id_hash}) <= 1) {
$self->update_single_peptide_tree($tree);
}
$self->{'tag_residue_count'} = $residues;
return $fastafile;
}
sub parse_and_store_alignment_into_proteintree
{
my $self = shift;
my $mcoffee_output = $self->{'mcoffee_output'};
my $mcoffee_scores = $self->{'mcoffee_scores'};
my $tree = $self->{'protein_tree'};
return unless($mcoffee_output and -e $mcoffee_output);
#
# Read in the alignment using Bioperl.
#
use Bio::AlignIO;
my $alignio = Bio::AlignIO->new(-file => "$mcoffee_output",
-format => "fasta");
my $aln = $alignio->next_aln();
my %align_hash;
foreach my $seq ($aln->each_seq) {
my $id = $seq->display_id;
$align_hash{$id} = $seq->seq;
}
#
# Read in the scores file manually.
#
my %score_hash;
if (defined $mcoffee_scores) {
my $FH = IO::File->new();
$FH->open($mcoffee_scores) || throw("Could not open alignment scores file [$mcoffee_scores]");
<$FH>; #skip header
my $i=0;
while(<$FH>) {
$i++;
next if ($i < 7); # skip first 7 lines.
next if($_ =~ /^\s+/); #skip lines that start with space
if ($_ =~ /:/) {
my ($id,$overall_score) = split(/:/,$_);
$id =~ s/^\s+|\s+$//g;
$overall_score =~ s/^\s+|\s+$//g;
print "___".$id."___".$overall_score."___\n";
next;
}
chomp;
my ($id, $align) = split;
$score_hash{$id} ||= '';
$score_hash{$id} .= $align;
}
$FH->close;
}
#
# Convert alignment strings into cigar_lines
#
my $alignment_length;
foreach my $id (keys %align_hash) {
next if ($id eq 'cons');
my $alignment_string = $align_hash{$id};
unless (defined $alignment_length) {
$alignment_length = length($alignment_string);
} else {
if ($alignment_length != length($alignment_string)) {
throw("While parsing the alignment, some id did not return the expected alignment length\n");
}
}
# Call the method to do the actual conversion
$align_hash{$id} = $self->_to_cigar_line(uc($alignment_string));
}
if (defined($self->{redo}) && $self->{'output_table'} eq 'protein_tree_member') {
# We clone the tree, attach it to the new clusterset_id, then store it.
# protein_tree_member is now linked to the new one
my ($from_clusterset_id, $to_clusterset_id) = split(":",$self->{'redo'});
throw("malformed redo option: ". $self->{'redo'}." should be like 1:1000000")
unless (defined($from_clusterset_id) && defined($to_clusterset_id));
my $clone_tree = $self->{protein_tree}->copy;
my $clusterset = $self->{treeDBA}->fetch_node_by_node_id($to_clusterset_id);
$clusterset->add_child($clone_tree);
$self->{treeDBA}->store($clone_tree);
# Maybe rerun indexes - restore
# $self->{treeDBA}->sync_tree_leftright_index($clone_tree);
$self->_store_aln_tags($clone_tree);
# Point $tree object to the new tree from now on
$tree->release_tree; $tree = $clone_tree;
}
#
# Align cigar_lines to members and store
#
foreach my $member (@{$tree->get_all_leaves}) {
# Redo alignment is member_id based, new alignment is sequence_id based
if ($align_hash{$member->sequence_id} eq "" && $align_hash{$member->member_id} eq "") {
throw("mcoffee produced an empty cigar_line for ".$member->stable_id."\n");
}
# Redo alignment is member_id based, new alignment is sequence_id based
$member->cigar_line($align_hash{$member->sequence_id} || $align_hash{$member->member_id});
## Check that the cigar length (Ms) matches the sequence length
# Take the M lengths into an array
my @cigar_match_lengths = map { if ($_ eq '') {$_ = 1} else {$_ = $_;} } map { $_ =~ /^(\d*)/ } ( $member->cigar_line =~ /(\d*[M])/g );
# Sum up the M lengths
my $seq_cigar_length; map { $seq_cigar_length += $_ } @cigar_match_lengths;
my $member_sequence = $member->sequence; $member_sequence =~ s/\*//g;
if ($seq_cigar_length != length($member_sequence)) {
print $member_sequence."\n".$member->cigar_line."\n" if ($self->debug);
throw("While storing the cigar line, the returned cigar length did not match the sequence length\n");
}
if ($self->{'output_table'} eq 'protein_tree_member') {
#
# We can use the default store method for the $member.
$self->{'comparaDBA'}->get_ProteinTreeAdaptor->store($member);
} else {
#
# Do a manual insert into the correct output table.
#
my $table_name = $self->{'output_table'};
printf("Updating $table_name %s : %s\n",$member->stable_id,$member->cigar_line) if ($self->debug);
my $sth = $self->{treeDBA}->prepare("INSERT ignore INTO $table_name
(node_id,member_id,method_link_species_set_id,cigar_line) VALUES (?,?,?,?)");
$sth->execute($member->node_id,$member->member_id,$member->method_link_species_set_id,$member->cigar_line);
$sth->finish;
}
if (defined $self->{'mcoffee_scores'}) {
#
# Do a manual insert of the *scores* into the correct score output table.
#
my $table_name = $self->{'output_table'} . "_score";
my $sth = $self->{treeDBA}->prepare("INSERT ignore INTO $table_name
(node_id,member_id,method_link_species_set_id,cigar_line) VALUES (?,?,?,?)");
my $score_string = $score_hash{$member->sequence_id} || '';
$score_string =~ s/[^\d-]/9/g; # Convert non-digits and non-dashes into 9s. This is necessary because t_coffee leaves some leftover letters.
printf("Updating $table_name %s : %s\n",$member->stable_id,$score_string) if ($self->debug);
$sth->execute($member->node_id,$member->member_id,$member->method_link_species_set_id,$score_string);
$sth->finish;
}
}
}
# Converts the given alignment string to a cigar_line format.
sub _to_cigar_line {
my $self = shift;
my $alignment_string = shift;
$alignment_string =~ s/\-([A-Z])/\- $1/g;
$alignment_string =~ s/([A-Z])\-/$1 \-/g;
my @cigar_segments = split " ",$alignment_string;
my $cigar_line = "";
foreach my $segment (@cigar_segments) {
my $seglength = length($segment);
$seglength = "" if ($seglength == 1);
if ($segment =~ /^\-+$/) {
$cigar_line .= $seglength . "D";
} else {
$cigar_line .= $seglength . "M";
}
}
return $cigar_line;
}
sub _store_aln_tags {
my $self = shift;
my $tree = shift || $self->{'protein_tree'};
my $output_table = $self->{'output_table'};
my $pta = $self->{'comparaDBA'}->get_ProteinTreeAdaptor;
print "Storing Alignment tags...\n";
#
# Retrieve a tree with the "correct" cigar lines.
#
if ($output_table ne "protein_tree_member") {
$tree = $self->_get_alternate_alignment_tree($pta,$tree->node_id,$output_table);
}
my $sa = $tree->get_SimpleAlign;
# Alignment percent identity.
my $aln_pi = $sa->average_percentage_identity;
$tree->store_tag("aln_percent_identity",$aln_pi);
# Alignment length.
my $aln_length = $sa->length;
$tree->store_tag("aln_length",$aln_length);
# Alignment runtime.
my $aln_runtime = int(time()*1000-$self->{'mcoffee_starttime'});
$tree->store_tag("aln_runtime",$aln_runtime);
# Alignment method.
my $aln_method = $self->{'method'};
$tree->store_tag("aln_method",$aln_method);
# Alignment residue count.
my $aln_num_residues = $sa->no_residues;
$tree->store_tag("aln_num_residues",$aln_num_residues);
# Alignment redo mapping.
my ($from_clusterset_id, $to_clusterset_id) = split(":",$self->{'redo'});
my $redo_tag = "MCoffee_redo_".$from_clusterset_id."_".$to_clusterset_id;
$tree->store_tag("$redo_tag",$self->{'protein_tree_id'}) if ($self->{'redo'});
}
sub _get_alternate_alignment_tree {
my $self = shift;
my $pta = shift;
my $node_id = shift;
my $table = shift;
my $tree = $pta->fetch_node_by_node_id($node_id);
foreach my $leaf (@{$tree->get_all_leaves}) {
# "Release" the stored / cached values for the alignment strings.
undef $leaf->{'cdna_alignment_string'};
undef $leaf->{'alignment_string'};
# Grab the correct cigar line for each leaf node.
my $id = $leaf->member_id;
my $cmd = "SELECT cigar_line FROM $table where member_id=$id;";
my $sth = $pta->prepare($cmd);
$sth->execute();
my $data = $sth->fetchrow_hashref();
$sth->finish();
my $cigar = $data->{'cigar_line'};
die "No cigar line for member $id!\n" unless ($cigar);
$leaf->cigar_line($cigar);
}
return $tree;
}
1;
]]>