#!/usr/bin/env perl

my %hashSequences;
my %hashTrimmedPos;
main();

sub main
{
    my $Usage = "$0 dirRealignment: Trims leading and succeeding gaps and N's.";

    my ($dirRealignment) = @ARGV;

    print STDERR "\n$Usage: $dirRealignment\n";
  
    opendir(DIR, $dirRealignment) or die "ERROR (main): Could not open directory '$dirRealignment'!\n\n "; 
    my @Files = readdir(DIR);
    close DIR;

    # Need to keep track of alignment positions which have been trimmed, in order to find the positions in the database
    #my $outFileTrimmedPos = $dirRealignment;
    #$outFileTrimmedPos =~ s/[\.\/]//g;
    #$outFileTrimmedPos .= ".trimmedPos";
    #open(TRIMMED_POS, ">$outFileTrimmedPos") or die "ERROR (main): Could not create file '$outFileTrimmedPos':$!\n\n";

    foreach my $file(sort {$a cmp $b} @Files)
    {
	$file = $dirRealignment.$file;
	if($file =~ /.*\.aln$/)
	{
	    %hashSequences = ();
	    %hashTrimmedPos = ();
	    my $alignLength = read_sequences($file);
	    my $outFile = $file.".trim";
	    my $trimmedAlignLength = trim_alignment_clustalw($outFile);

	    # Remove alignment file due to space problems! - Remove it in script calling TrimRealignment.pl
 	    #system("rm $file") == 0 or die "system 'rm $file' failed: $?";

	    # Save trimmed positions
	    print ">$file\n";
	    my @KeysTrimmedPos = keys %hashTrimmedPos;
	    foreach my $kTrimmedPos (sort {$a cmp $b} @KeysTrimmedPos)
	    {
		print "SEQ|$kTrimmedPos| ";
		my @seqTrimmedPos = @{$hashTrimmedPos{$kTrimmedPos}};
		foreach my $s(@seqTrimmedPos)
		{
		    print "$s ";
		}
		print "\n";
		
	    }
	}
    }
    
    #close(TRIMMED_POS);

}# end of sub main


# PURPOSE:  Trim leading and preceeding gaps of multiple alignment and write it out in file in CLUSTALW format.
#           Trimmed alignment starts with the first contiguous block of aligned letters in all sequences and 
#           ends with such a block! Regard N's or n's also as leading or preceeding gaps.
sub trim_alignment_clustalw
{
    my ($file) = @_;

    my @Keys = keys %hashSequences;

    ## FIND POSITIONS WHICH WILL BE TRIMMED
    # Sequences have same length as they start and end with matching regions!
    my $length = length $hashSequences{$Keys[0]};
    my $startTrim = 0;
    for(my $i=0; $i<$length; $i++)
    {
	my $column = get_column($i, %hashSequences);
	# As long column contains a gap continue
	if($column =~ /.*[\-Nn]+.*/){ next; }
	else  # First match in all sequences found
	{
	    $startTrim = $i;
	    last;
	}
    }
    # Same from the end of the alignment
    my $endTrim = $length;
    for(my $i=$length-1; $i>=0; $i--)
    {
	my $column = get_column($i, %hashSequences);
	# As long column contains a gap continue
	if($column =~ /.*[\-Nn]+.*/){ next; }
	else  # First match in all sequences found
	{
	    $endTrim = $i;
	    last;
	}
    }

    ## SAVE TRIMMED POSITIONS (before trimming!)
    # At last, keep track of the trimmed positions!
    foreach my $k(@Keys)
    {
	# Here, do not reset N's or n's as they are actually part of the sequence!
	my $seqStart = substr $hashSequences{$k}, 0, $startTrim;
	$seqStart =~ s/-//g;
	my $seqEnd = substr $hashSequences{$k}, $endTrim;
	$seqEnd =~ s/-//g;

	my $cntStart = length $seqStart;
	my $cntEnd = (length $seqEnd) - 1;

	push(@{$hashTrimmedPos{$k}}, $cntStart);
	push(@{$hashTrimmedPos{$k}}, $cntEnd);
	
    }

    ## TRIM SEQUENCES
    #print "startTrim = $startTrim  endTrim = $endTrim\n";
    foreach my $k(@Keys)
    {
	$hashSequences{$k} = substr $hashSequences{$k}, $startTrim, $endTrim-$startTrim+1;
    }
    my $trimmedLength = length $hashSequences{$Keys[0]};


    # Now write trimmed sequences into input file for RNAz
    open(OUT, ">$file") or die "Could not create file $file!\n\n";
    print OUT "CLUSTAL W\n\n";
 
    my $line_length = 60;
    for(my $i=0; $i<$trimmedLength; $i += $line_length)
    {
        foreach my $k(sort {$a cmp $b} @Keys)
        {
            my $seq = $hashSequences{$k};
	    if($seq eq '') # Nachtrag 19.01.05
	    {
		die "ERROR (trim_alignment_clustalw): Could not find sequence '$k'!\n\n";
	    }
            if(length $seq != $trimmedLength)
            {
                die "ERROR (trim_alignment_clustalw): Length of sequence '$k' (", length $seq, ") != $length!\n\n";
            }
            print OUT "$k\t";
            print OUT substr($seq, $i, $line_length), "\n";
        }
        print OUT "\n\n";
    }
 
    close(OUT);

    return $trimmedLength;

}# end of trim_alignment


# PURPOSE: Gets the column of an alignment as a string
# PARAMS:  $c: Index of column [0,MAX]
sub get_column
{
    my ($c, %hashAlignSeq) = @_;
                                                                                                                                                   
    my $column = '';
                                                                                                                                                   
    my @Keys = keys %hashAlignSeq;
                                                                                                                                                   
    foreach my $k(@Keys)
    {
        $column .= substr $hashAlignSeq{$k}, $c, 1;
    }
                                                                                                                                                   
    return $column;
                                                                                                                                                   
}# end of sub get_column



sub read_sequences
{
    my ($file) = @_;

    open(IN, "$file") or die "ERROR (read sequences): Could not open file '$file'\n\n";
    my @filedata = <IN>;
    close(IN);

    my $seq = '';
    my $id = '';
    my $cnt = 0;
    foreach my $line(@filedata)
    {
	if($line =~ /^\s*CLUSTAL W\s*$/)
	{
	    next;
	}
	elsif($line =~ /^\s*$/)
	{
	    $cnt = 0;       # Alignments are spread over several lines. A blank line indicates 
	                    # new line of multiple alignment!
	    next;
	}
	elsif($line =~ /^\s*[\*\s]+\s*$/)
	{
	    next;
	}
	elsif($line =~ /^\s*(\S*)\s+(\S*)\s*$/)
	{
	    $cnt++;
	    $id = $cnt.".".$1;  # To make sure that IDs are unique as DIALIGN cuts off IDs
	    if(!exists $hashSequences{$id}){ $hashSequences{$id}='';}
	    my $subseq = $2;
	    $subseq =~ s/\s//g;
	    $hashSequences{$id} .= $subseq;
	    #print "-$id- $hashSequences{$id}\n";
	}
   
    }

    my @Keys = keys %hashSequences;
    my $l = length $hashSequences{$Keys[0]};

    #print "$Keys[0] $l\n";
    return $l;

}# end of sub read_sequences
