[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.1 - (download) (as text) (annotate)
Wed Nov 30 08:44:26 2005 UTC (14 years, 6 months ago) by golsen
Branch: MAIN
Introduce a different tree printing routine to fid_checked.  This required
a new version of align_with_clustal.pl (= align_with_clustal_2.pl) and a
new tree utilities library (gjonewicklib.pm).

The new tree printer does not require deleting the extra lines (they are
not included in the first place).  This would probably allow some other
restructuring of the code in fid_checked, but I did the minimum.

The other reason for the changes (actually the real reason) is that it is
now possible to introduce midpoint rooting to the tree, and to reorder
the taxa so that the tree is MUCH easier to read (aesthetic tree order).

Also I'm changing it to use a second-generation tree from the alignment,
not just show the guide tree.

# -*- 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 "
     . ( $tree ? "-tree " : "" )
     . "-outorder=aligned > /dev/null";
if ( -s "$file.aln" && open( ALIGN, "<$file.aln" ) )
{
    print <ALIGN>;
    close( ALIGN );
}

if ( $tree )
{
    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