libbio-perl-perl 1.6.923-1 / usr / share / perl5 / Bio / AlignIO / nexus.pm

#
# BioPerl module for Bio::AlignIO::nexus
#
# Copyright Heikki Lehvaslaiho
#

=head1 NAME

Bio::AlignIO::nexus - NEXUS format sequence input/output stream

=head1 SYNOPSIS

Do not use this module directly.  Use it via the L<Bio::AlignIO> class.

    use Bio::AlignIO;

    my $in = Bio::AlignIO->new(-format => 'nexus',
                              -file   => 'aln.nexus');
    while( my $aln = $in->next_aln ) {
        # do something with the alignment
    }

=head1 DESCRIPTION

This object can transform L<Bio::Align::AlignI> objects to and from NEXUS
data blocks. See method documentation for supported NEXUS features.

=head1 ACKNOWLEDGEMENTS

Will Fisher has written an excellent standalone NEXUS format parser in
Perl, readnexus. A number of tricks were adapted from it.

=head1 FEEDBACK

=head2 Support 

Please direct usage questions or support issues to the mailing list:

I<bioperl-l@bioperl.org>

rather than to the module maintainer directly. Many experienced and 
reponsive experts will be able look at the problem and quickly 
address it. Please include a thorough description of the problem 
with code and data examples if at all possible.

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution.  Bug reports can be submitted via the
web:

  https://redmine.open-bio.org/projects/bioperl/

=head1 AUTHORS - Heikki Lehvaslaiho

Email: heikki-at-bioperl-dot-org

=head1 APPENDIX

The rest of the documentation details each of the object
methods. Internal methods are usually preceded with a _

=cut

# Let the code begin...

package Bio::AlignIO::nexus;
use vars qw(%valid_type);
use strict;
no strict "refs";


use base qw(Bio::AlignIO);

BEGIN {
    %valid_type = map {$_, 1} qw( dna rna protein standard );
    # standard throws error: inherited from Bio::PrimarySeq
}

=head2 new

 Title   : new
 Usage   : $alignio = Bio::AlignIO->new(-format => 'nexus', -file => 'filename');
 Function: returns a new Bio::AlignIO object to handle clustalw files
 Returns : Bio::AlignIO::clustalw object
 Args    : -verbose => verbosity setting (-1,0,1,2)
           -file    => name of file to read in or with ">" - writeout
           -fh      => alternative to -file param - provide a filehandle 
                       to read from/write to 
           -format  => type of Alignment Format to process or produce

           Customization of nexus flavor output

           -show_symbols => print the symbols="ATGC" in the data definition
                            (MrBayes does not like this)
                            boolean [default is 1] 
           -show_endblock => print an 'endblock;' at the end of the data
                            (MyBayes does not like this)
                            boolean [default is 1] 

=cut

sub _initialize {
    my ($self, @args) = @_;
    $self->SUPER::_initialize(@args);
    my ($show_symbols, $endblock) = 
	$self->_rearrange([qw(SHOW_SYMBOLS SHOW_ENDBLOCK)], @args);
    my @names = qw(symbols endblock);
    for my $v ( $show_symbols, $endblock ) {
	$v = 1 unless defined $v; # default value is 1
	my $n = shift @names;
	$self->flag($n, $v);
    }
}


=head2 next_aln

 Title   : next_aln
 Usage   : $aln = $stream->next_aln()
 Function: Returns the next alignment in the stream.

           Supports the following NEXUS format features:
           - The file has to start with '#NEXUS'
           - Reads in the name of the alignment from a comment
             (anything after 'TITLE: ') .
           - Sequence names can be given in a taxa block, too.
           - If matchchar notation is used, converts
             them back to sequence characters.
           - Does character conversions specified in the
             NEXUS equate command.
           - Sequence names of type 'Homo sapiens' and
             Homo_sapiens are treated identically.

 Returns : L<Bio::Align::AlignI> object
 Args    :


=cut


sub next_aln {
    my $self = shift;
    my $entry;
    my ($aln_name, $seqcount, $residuecount, %hash, $alphabet,
	$match, $gap, $missing, $equate, $interleave,
	$name,$str,@names,$seqname,$start,$end,$count,$seq);
    local $Bio::LocatableSeq::OTHER_SYMBOLS = '\*\?\.';
    local $Bio::LocatableSeq::GAP_SYMBOLS = '\-';
    my $aln =  Bio::SimpleAlign->new(-source => 'nexus');

    # file starts with '#NEXUS' but we allow white space only lines before it
    $entry = $self->_readline;
    $entry = $self->_readline while defined $entry && $entry =~ /^\s+$/;

    return unless $entry;
    $self->throw("Not a valid interleaved NEXUS file! [#NEXUS] not starting the file\n$entry")
	unless ($entry && $entry =~ /^#NEXUS/i);

    # skip anything before either the taxa or data block
    # but read in the optional title in a comment
    while (defined($entry = $self->_readline)) {
	local ($_) = $entry;
	/\[TITLE. *([^\]]+)]\s+/i and $aln_name = $1;
	last if /^begin +data/i || /^begin +taxa/i;
    }
    $aln_name =~ s/\s/_/g and $aln->id($aln_name) if $aln_name;

    # data and taxa blocks
    my $incomment;
    while (defined ($entry = $self->_readline)) {
	local ($_) =  $entry;
	next if s/\[[^\]]+\]//g; # remove comments
	if( s/\[[^\]]+$// ) {
	    $incomment = 1;
		 # skip line if it is now empty or contains only whitespace
	    next if /^\s*$/;
	} elsif($incomment) {
	    if( s/^[^\]]*\]// ) {
			 $incomment = 0;
	    } else {
			 next;
	    }
	} elsif( /taxlabels/i ) {
	    # doesn't deal with taxlabels adequately and can mess things up!
	    # @names = $self->_read_taxlabels;
	} else {

	    /ntax\s*=\s*(\d+)/i        and $seqcount = $1;
	    /nchar\s*=\s*(\d+)/i       and $residuecount = $1;
	    /matchchar\s*=\s*(.)/i     and $match = $1;
	    /gap\s*=\s*(.)/i           and $gap = $1;
	    /missing\s*=\s*(.)/i       and $missing = $1;
	    /equate\s*=\s*\"([^\"]+)/i and $equate = $1;  # "e.g. equate="T=C G=A";
	    /datatype\s*=\s*(\w+)/i    and $alphabet = lc $1;
	    /interleave/i              and $interleave = 1 ;
	    last if /matrix/io;
	}
    }
    $self->throw("Not a valid NEXUS sequence file. Datatype not specified.")
	unless $alphabet;
    $self->throw("Not a valid NEXUS sequence file. Datatype should not be [$alphabet]")
	unless $valid_type{$alphabet};
    $self->throw("\"$gap\" is not a valid gap character. For compatability, gap char can not be one of: ()[]{}/\,;:=*'`\"<>^")
    	if $gap && $gap =~ /[\(\)\[\]\{\}\/\\\,\;\:\=\*\'\`\<\>\^]/;
    $self->throw("\"$missing\" is not a valid missing character. For compatability, missing char can not be one of: ()[]{}/\,;:=*'`\"<>^")
    	if $missing && $missing =~ /[\(\)\[\]\{\}\/\\\,\;\:\=\*\'\`\<\>\^]/;

    $aln->gap_char($gap);
    $aln->missing_char($missing);

    #
    # if data is not right after the matrix line
    #  read the empty lines out
    #
    while ($entry = $self->_readline) {
	unless ($entry =~ /^\s+$/) {
	    $self->_pushback($entry);
	    last;
	}
    }

    #
    # matrix command
    #
    # first alignment section
    if (@names == 0) {		# taxa block did not exist
	while ($entry = $self->_readline) {
		local ($_) =  $entry;
		if( s/\[[^\]]+\]//g ) { #]  remove comments
			next if /^\s*$/; 
			# skip line if it is now empty or contains only whitespace
		}
		if ($interleave && defined$count && ($count <= $seqcount)) {
			/^\s+$/ and last;
		} else {
			/^\s+$/ and next;
		}
		/^\s*;/ and last;	# stop if colon at end of matrix is on it's own line
		#/^\s*;\s*$/ and last;
		if ( /^\s*([\"\'](.+?)[\"\']|(\S+))\s+(.*)\s*$/ ) {	
			# get single and double quoted names, or all the first 
         # nonwhite word as the name, and remained is seq
			#if (/^\s*('([^']*?)'|([^']\S*))\s+(.*)$/) { #'
			$name = ($2 || $3);
			if  ($4) {
				# seq is on same line as name
				# this is the usual NEXUS format
				$str = $4;
			} else {
				# otherwise get seq from following lines. No comments allowed
				# a less common matrix format, usually used for very long seqs
				$str='';
				while (local ($_) = $self->_readline) {
					my $str_tmp = $_;
					$str_tmp =~ s/[\s;]//g;
					$str .= $str_tmp;
					last if length$str == $residuecount;
				}
			}
			$name =~ s/ /_/g;
			push @names, $name;

			$str =~ s/[\s;]//g;
			$count =  @names;
			$hash{$count} = $str;
		}
		$self->throw("Not a valid interleaved NEXUS file! seqcount [$count] > predeclared [$seqcount] in the first section") if $count > $seqcount;
		/;/ and last;	# stop if colon at end of matrix is on the same line as the last seq
	}
}

    # interleaved sections
    $count = 0;
    if ( $interleave ) {	# only read next section if file is interleaved
	while( $entry = $self->_readline) {
	    local ($_) =  $entry;
	    if( s/\[[^\]]+\]//g ) { #]  remove comments
		next if /^\s*$/; # skip line if it is now empty or contains only whitespace
	    }
	    /^\s*;/ and last;		# stop if colon at end of matrix is on it's own line
	    $count = 0, next if $entry =~ /^\s*$/;
	    if (/^\s*(\'([^\']*?)\'|([^\']\S*))\s+(.*)$/) { 
		$str = $4;
		$str =~ s/[\s;]//g;
		$count++;
		$hash{$count} .= $str;
	    };
	    $self->throw("Not a valid interleaved NEXUS file!
    		seqcount [$count] > predeclared [$seqcount] ") if $count > $seqcount;
	    /;/ and last;	# stop if colon at end of matrix is on the same line as the last seq
	}
    }
    
    return 0 if @names < 1;
    
    # sequence creation
    $count = 0;
    foreach $name ( @names ) {
	$count++;
	if( $name =~ /(\S+)\/(\d+)-(\d+)/ ) {
	    ($seqname,$start,$end) = ($1,$2,$3);
	} else {
	    ($seqname,$start,$str) = ($name,1,$hash{$count});
	    $str =~ s/[$Bio::LocatableSeq::GAP_SYMBOLS]//g;
	    $end = length($str);
	}
	
	# consistency test
	$self->throw("Length of sequence [$seqname] is not [$residuecount]; got".CORE::length($hash{$count}))
	    unless CORE::length($hash{$count}) == $residuecount;
	
	$seq = Bio::LocatableSeq->new('-seq'        => $hash{$count},
				      '-display_id' => $seqname,
				      '-start'      => $start,
				      '-end'        => $end,
				     '-alphabet'    => $alphabet
				      );
	$aln->add_seq($seq);
    }

    # if matchchar is used
    $aln->unmatch($match) if $match;

    # if equate ( e.g. equate="T=C G=A") is used
    if ($equate) {
	$aln->map_chars($1, $2) while $equate =~ /(\S)=(\S)/g;
    }

    while  (defined $entry &&
	    $entry !~ /endblock/i) {
        $entry = $self->_readline;
    }

    return $aln if $aln->num_sequences;
	return;
}

sub _read_taxlabels {
    my ($self) = @_;
    my ($name, @names);
    while (my $entry = $self->_readline) {
	last if $entry =~ m/^\s*(END)?;/i;
	if( $entry =~ m/\s*(\S+)\s+/ ) {
	    ($name) = ($1);
	    $name =~ s/\[[^\[]+\]//g;
	    $name =~ s/\W/_/g;
	    push @names, $name;
	}
    }
    return @names;
}

=head2 write_aln

 Title   : write_aln
 Usage   : $stream->write_aln(@aln)
 Function: Writes the $aln object into the stream in interleaved NEXUS
           format. Everything is written into a data block.
           SimpleAlign methods match_char, missing_char and gap_char must be set
           if you want to see them in the output.
 Returns : 1 for success and 0 for error
 Args    : L<Bio::Align::AlignI> object

=cut

sub write_aln {
    my ($self,@aln) = @_;
    my $count = 0;
    my $wrapped = 0;
    my $maxname;
    my ($length,$date,$name,$seq,$miss,$pad,%hash,@arr,$tempcount,$index );
    my ($match, $missing, $gap,$symbols) = ('', '', '','');

    foreach my $aln (@aln) {
	if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) {
	    $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
	    next;
	}
	$self->throw("All sequences in the alignment must be the same length")
	    unless $aln->is_flush($self->verbose);

	$length  = $aln->length();

	$self->_print (sprintf("#NEXUS\n[TITLE: %s]\n\nbegin data;\ndimensions ntax=%s nchar=%s;\n",
			       $aln->id, $aln->num_sequences, $length));
	$match = "match=". $aln->match_char if $aln->match_char;
	$missing = "missing=". $aln->missing_char if $aln->missing_char;
	$gap = "gap=". $aln->gap_char if $aln->gap_char;

	$symbols = 'symbols="'.join('',$aln->symbol_chars). '"' 
	    if( $self->flag('symbols') && $aln->symbol_chars);
	$self->_print 
	    (sprintf("format interleave datatype=%s %s %s %s %s;\n\nmatrix\n",
		     $aln->get_seq_by_pos(1)->alphabet, $match, 
		     $missing, $gap, $symbols));

                     # account for single quotes round names
	my $indent = $aln->maxdisplayname_length+2;

	$aln->set_displayname_flat();
	foreach $seq ( $aln->each_seq() ) {
	    my $nmid = $aln->displayname($seq->get_nse());
	    if( $nmid =~ /[^\w\d\.]/ ) {
              # put name in single quotes incase it contains any of
              # the following chars: ()[]{}/\,;:=*'"`+-<> that are not
              # allowed in PAUP* and possible other software

		$name = sprintf("%-${indent}s", "\'" . $nmid . "\'");
	    } else { 
		$name = sprintf("%-${indent}s", $nmid);
	    }
	    $hash{$name} = $seq->seq;
	    push(@arr,$name);
	}

	while( $count < $length ) {
	    # there is another block to go!
	    foreach $name ( @arr ) {
		my $dispname = $name;
#		 $dispname = '' if $wrapped;
		$self->_print (sprintf("%${indent}s  ",$dispname));
		$tempcount = $count;
		$index = 0;
		while( ($tempcount + 10 < $length) && ($index < 5)  ) {
		    $self->_print (sprintf("%s ",substr($hash{$name},$tempcount,10)));
		    $tempcount += 10;
		    $index++;
		}
		# last
		if( $index < 5) {
		    # space to print!
		    $self->_print (sprintf("%s ",substr($hash{$name},$tempcount)));
		    $tempcount += 10;
		}
		$self->_print ("\n");
	    }
	    $self->_print ("\n\n");
	    $count = $tempcount;
	    $wrapped = 1;
	}
	if( $self->flag('endblock') ) {
	    $self->_print (";\n\nendblock;\n");
	} else { 
	    $self->_print (";\n\nend;\n");
	}
    }
    $self->flush if $self->_flush_on_write && defined $self->_fh;
    return 1;
}

=head2 flag

 Title   : flag
 Usage   : $obj->flag($name,$value)
 Function: Get/Set a flag value
 Returns : value of flag (a scalar)
 Args    : on set, new value (a scalar or undef, optional)


=cut

sub flag{
    my ($self,$name,$val) = @_;
    return $self->{'flag'}->{$name} = $val if defined $val;
    return $self->{'flag'}->{$name};
}

1;