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

BioPerl documentation Contained in the BioPerl distribution.

Index

Code Index:

NAME

Top

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

SYNOPSIS

Top

Do not use this module directly. Use it via the 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
    }

DESCRIPTION

Top

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

ACKNOWLEDGEMENTS

Top

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

FEEDBACK

Top

Support

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

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.

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/

AUTHORS - Heikki Lehvaslaiho

Top

Email: heikki-at-bioperl-dot-org

APPENDIX

Top

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

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]

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    :

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

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)
BioPerl documentation Contained in the BioPerl distribution. 
#
# BioPerl module for Bio::AlignIO::nexus
#
# Copyright Heikki Lehvaslaiho
#


# 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
}


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);
    }
}




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;
}


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;
}


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

1;