[Bio] / FigKernelScripts / align_with_clustal_2.pl Repository:
ViewVC logotype

View of /FigKernelScripts/align_with_clustal_2.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (as text) (annotate)
Wed Nov 30 09:18:08 2005 UTC (14 years, 5 months ago) by golsen
Branch: MAIN
Changes since 1.1: +4 -5 lines
Fix the clustalw command so that it really uses the final alignment
for the tree (it seems that it is not possible to combine the -align
and -tree command flags in a single command, so it is now run as two
commands).

# -*- perl -*-

use strict;
use Carp;
use Data::Dumper;
use gjonewicklib qw(
		parse_newick_tree_str
		reroot_newick_to_approx_midpoint
		aesthetic_newick_tree
		newick_relabel_tips
		text_plot_newick
		dump_tree
		);
use FIG;

my $fig = new FIG;

$| = 1;

my $usage = "align_with_clustal_2 [-org] [-func[=user]] [-tree] [-UniProt] [-save=Dir] Id1 Id2 ... ";

(@ARGV > 1) || die "usage: $usage";

my $temp_dir = $FIG_Config::temp;
my $file     = "$temp_dir/tmp$$";
my $relabel  = {};
my $add_org = "";
my $save = "";
my $add_func = "";
my $tree = "";
my $uni_prot = "";
my $trouble = 0;
my $user = undef;

while ($ARGV[0] =~ m/^-/)
{
    if ($ARGV[0] =~ /-org/)
    {
	$add_org = shift @ARGV;
    }
    elsif ($ARGV[0] =~ /-save=(\S+)/)
    {
	$save = $1;
	shift @ARGV;
    }
    elsif ($ARGV[0] =~ /-func/)
    {
	$add_func = shift @ARGV;
	if ($add_func =~ m/=(\S+)/)  { $user = $1; }
    }
    elsif ($ARGV[0] =~ /-tree/)
    {
	$tree =	shift @ARGV;
    }
    elsif ($ARGV[0] =~ /-UniProt/i)
    {
	$uni_prot = 1;
	shift @ARGV;
    }
    else
    {
	$trouble = 1;
	print STDERR "Invalid flag $ARGV[0]\n";
	shift @ARGV;
    }
}
die "\nusage: $usage\n\n" if $trouble;

my %seen;
open(TMP,">$file.fasta") || die "could not open $file";
foreach my $id (@ARGV)
{
    next if ($seen{$id});
    $seen{$id} = 1;
    
    if (my $seq = $fig->get_translation($id))
    {
	print TMP ">$id\n$seq\n";
    }
    else
    {
	print STDERR "could not find translation for $id\n";
    }

    my $label = $id;
    if ($uni_prot && ($_ = $fig->to_alias($id,"uni")) )      { $label .= " [$_]" }
    if ($add_org  && ($_ = $fig->org_of($id)) )              { $label .= " [$_]"; }
    if ($add_func && ($_ = $fig->function_of($id, $user)) )  { $label .= " $_"; }
    
    $relabel->{$id} = $label;
}
close(TMP);
# print STDERR "$0;\n", Dumper($relabel);

system "$FIG_Config::ext_bin/clustalw -infile=$file.fasta -align -outorder=aligned > /dev/null";
if ( -s "$file.aln" && open( ALIGN, "<$file.aln" ) )
{
    print <ALIGN>;
    close( ALIGN );
}

if ( $tree && -s "$file.aln" )
{
	system "$FIG_Config::ext_bin/clustalw -infile=$file.aln -tree > /dev/null";
    my $treetext = "";
    my $tree1;
    if ( ( -s "$file.ph"  && open( TREE, "<$file.ph"  ) )
      || ( -s "$file.dnd" && open( TREE, "<$file.dnd" ) ) )
    {
        if ( $tree1 = &parse_newick_tree_str( join( "", <TREE> ) ) )
        {
            # dump_tree( $tree1 );
            my $tree1b = reroot_newick_to_approx_midpoint( $tree1 );
            my $tree2 = aesthetic_newick_tree( $tree1b );
            # my $tree2 = aesthetic_newick_tree( $tree1 );
            my $tree3 = newick_relabel_tips( $tree2, $relabel );
            #  @textlines = text_plot_newick( $node, $width, $min_dx, $dy )
            $treetext = join( "\n", text_plot_newick( $tree3, 80, 2, 2 ) );
        }
        close( TREE );
    }
    print "=======================================================================\n\n\n";
    print "$treetext\n";
}

if ($save)
{
    system "cp $file.fasta \"$save/seqs\"";
    system "cp $file.aln \"$save/aln\"";
    system "cp $file.dnd \"$save/tree\"";
    ( $tree && -s "$file.ph" ) and system "cp $file.ph \"$save/tree2\"";
}

unlink("$file.fasta","$file.aln","$file.dnd");
# -s "$file.ph" and unlink("$file.ph");

1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3