#!/usr/bin/perl -w
##### Program description #######
#
# Title: convert.pl
#
# Author(s): Mattias Ohlsson
#
# Description: 
#  This program converts between the phylip and clustalw formats
#  for multiple sequence alignment. 
#  Note: This program does not handle the conservation line appearing
#        in the clustalw format.
#
# List of subroutines:
#  ReadPhylip();
#  ReadClustalW();
#  WritePhylip();
#  WriteClustalW();
#
# Overall procedure:
#  Each sequence is stored as a long string of letters (including gaps -).
#  Since the two formats have different numbers of AA's on each line
#  one cannot read and convert at the same time. The flow of the program
#  is a s follows:
#  
#  1. Read command line arguments to get the type of conversion and the
#     file to convert.
#  2. Read the file and store as an array of lines.
#  3. Parse the file in order to store each sequence as a long string of letters.
#  4. Write the alignement in the other format.
#
# Usage:
#  convert.pl {p2c|c2p} {filepath}
#
#  where p2c means conversion phylip -> clustalw and
#  c2p clustalw -> phylip. The filepath is the file to convert.
#
##################################
# CUT #
use strict;
our(%Seq, @Species);  # Global variables


#### The main program ####
open my $FH, '<', $ARGV[1] or die "Cannot open file $ARGV[1]: $!\n";
my @indata = <$FH>;
close $FH;

if( uc($ARGV[0]) eq 'P2C' ) {
    ReadPhylip(@indata);
    WriteClustalW();
} 
elsif( uc($ARGV[0]) eq 'C2P' ) {
    ReadClustalW(@indata);
    WritePhylip();
} 
else {
    die "Error: Unknown conversion mode!\n";
}

# CUT #
sub ReadPhylip {
    
    # Copy the supplied data to a local data array
    my @data = @_;

    # Use the first line to find the number of proteins
    my ($nprot, $plen) = split ' ', shift @data;
    
    # Read the first $nprot lines to get the species and the start of the seq
    my @tmp;
    for ( my $i = 0; $i < $nprot; $i++ ) {
	($Species[$i], @tmp) = split ' ', shift @data;
	chomp @tmp;
	$Seq{ $Species[$i] } = join '', @tmp;
    }
   
    # Continue to read lines of sequences
    my $np = 0;
    while ( my $line = shift @data ) {
	next if ( $line =~ /^\s*$/ );

	### TO BE COMPLETED ###

	$Seq{ $Species[$np] } .=  $line;

	$np += 1;
	if ( $np == $nprot ) {$np = 0;}
    }
}
# CUT #
sub ReadClustalW {
    
    # Copy the supplied data to a local data array
    my @data = @_;

    # Skip starting blank lines
    while ( $data[0] =~ /^\s*$/ ) {shift @data;}
	
    # The first non-blank line should contain the word CLUSTAL
    unless ( $data[0] =~ /CLUSTAL/ ) {
	die "The input file does not contain the word CLUSTAL!\n";
    }
    shift @data;

    # Store the sequenses in a hash
    my %aux;
    while ( my $line = shift @data ) {
	next unless ( $line =~ /(-|\w)/ );

	my ($name, $seqtmp, $lentmp) = split ' ', $line;
	$Seq{$name} .= $seqtmp;
	
	# Find the name in correct order
	unless ( defined $aux{$name} ) {
	    push @Species, $name;
	    $aux{$name} = 1;
	}
    }
}
# CUT #
sub WritePhylip {

    # Print the header line
    printf "%d\t%d\n", scalar @Species, length($Seq{ $Species[0] });
    
    # Now the rest of the sequences
    my $PrintName = 1;
    my $nprot = scalar @Species;
    while ( length($Seq{ $Species[0] })  > 0 ) {

	# The names
	foreach my $i (0 .. $nprot-1) {
	    my $name = $Species[$i];

	    if( $PrintName == 1 ) {
		printf "%10s\t", $name;
	    } else {
		printf "%10s\t", ' ';	
	    }

	    # The sequences
	    $Seq{$name} =~ s/^(.{1,50})//; 
	    my $aas_tmp = $1;

	    ### TO BE COMPLETED ###

	    print "$aas\n";
	}
	print "\n";
	$PrintName = 0;
    }
}
# CUT #
sub WriteClustalW {

    # Write the ClustalW header 
    print "CLUSTAL W (1.82) multiple sequence alignment\n\n\n";

    my $nprot = scalar @Species;
    my @aacount;
    while ( length($Seq{ $Species[0] }) > 0 ) {
	foreach my $i (0 .. $nprot-1) {
	    printf "%10s\t", $Species[$i];
	    $Seq{ $Species[$i] } =~ s/^(.{1,60})//; 
	    my $aas = $1;
	    print "$aas";
	    
	    # Count the number of non - in the aas
	    my $Naa = 0;
	    while ( $aas =~ /\w/g ) {
	       $Naa += 1;
	    }
	    $aacount[$i] += $Naa; 
	    if ( $aacount[$i] > 0 ) {
		print " $aacount[$i]\n";
	    }
	}
	print "\n\n";
    }
}
