view release on metacpan or  search on metacpan

lib/Bio/MUST/Core/Ali.pm  view on Meta::CPAN

            # seq_id are mandatory in first block
            croak "[BMC] Error: missing id in PHILIP file at line $.; aborting!"
                unless $seq_id;

            # store first seq chunk along with seq_id
            my $new_seq = Seq->new( seq_id => $seq_id, seq => $seq );
            $ali->add_seq($new_seq);
        }

        # process remaining seq blocks
        else {
            my $curr_seq = $ali->get_seq($n);

            # elongate current with partial seq chunk mistaken as seq_id
            if ($seq_id && $seq_id ne $curr_seq->full_id) {
                $curr_seq->append_seq($seq_id);
            }

            # elongate current seq with new seq chunk
            $curr_seq->append_seq($seq);
        }

        # prepare reading of next seq or next block
        $n = 0 if ++$n == $seq_n;
    }

    my $width = $ali->width;
    croak "[BMC] Error: unexpected site number in PHYLIP file: $width;"
        . ' aborting!' if $width != $site_n;

    return $ali;
}

# PHYLIP: http://evolution.genetics.washington.edu/phylip/doc/main.html
# The information for each species follows, starting with a ten-character
# species name (which can include blanks and some punctuation marks), and
# continuing with the characters for that species. The name should be on the
# same line as the first character of the data for that species. (...) The
# name should be ten characters in length, filled out to the full ten
# characters by blanks if shorter. Any printable ASCII/ISO character is
# allowed in the name, except for parentheses ("(" and ")"), square brackets
# ("[" and "]"), colon (":"), semicolon (";") and comma (","). (...) Note that
# in these sequences we have a blank every ten sites to make them easier to
# read: any such blanks are allowed. The blank line which separates the two
# groups of lines (the ones containing sites 1-20 and ones containing sites
# 21-39) may or may not be present.

# PhyML: http://www.atgc-montpellier.fr/phyml/usersguide.php?type=command
# The input sequence file is a standard PHYLIP file of aligned DNA or
# amino-acids sequences. (...) The maximum number of characters in species
# name MUST not exceed 100. Blanks and the symbols "(),:" are not allowed
# within sequence names because the Newick tree format makes special use of
# these symbols. However, blanks (one or more) MUST appear at the end of each
# species name.

# TREE-FINDER: http://www.treefinder.de/tf-march2011-manual.pdf
# A sequence name may consist of 1 to 10 alphanumeric characters, dashes "-",
# dots ".", underscores "_", or some of "/", "?", "*", "+". No space is
# allowed inside the names. The first name character must be a letter or an
# underscore. A sequence fragment may start behind position 10 after a
# sequence name, or anywhere in a line without a name.

# TREE-PUZZLE [data/globin.a]
#   7 128
# HBB_HUMAN      HLTPEEKSAV TALWGKVNVD EVGGEALGRL LVVYPWTQRF FESFDLSMGN
# HBB_HORSE      QLSGEEKAAV LALWDKVNEE EVGGEALGRL LVVYPWTQRF FDSFDLSMGN
# HBA_HUMAN      VLSPADKTNV KAAWGKVGAG EYGAEALERM FLSFPTTKTY FPHFDLSHGS
# HBA_HORSE      VLSAADKTNV KAAWSKVGAG EYGAEALERM FLGFPTTKTY FPHFDLSHGS
# MYG_PHYCA      VLSEGEWQLV LHVWAKVEVA GHGQDILIRL FKSHPETLEK FDRFHLKKAS
# GLB5_PETMA     PLSAAEKTKI RSAWAPVYYE TSGVDILVKF FTSTPAAQEF FPKFGLTKKS
# LGB2_LUPLU     ALTESQAALV KSSWEEFNIP KHTHRFFILV LEIAPAAKDL FSFLGTSQNN

# RAXML: http://sco.h-its.org/exelixis/oldPage/RAxML-Manual.7.0.4.pdf
# The input alignment format of RAxML is relaxed interleaved or sequential
# PHYLIP. "Relaxed" means that sequence names can be of variable length
# between 1 up to 256 characters. (...) Prohibited Character(s) in taxon names
# taxon names that contain any form of whitespace character, like blanks,
# tabulators, and carriage returns, as well as one of the following prohibited
# characters: :,();[].

# PHYLOBAYES: http://megasun.bch.umontreal.ca/People/lartillot/www/phylobayes3.3e.pdf
# Taxon names may contain more than 10 characters. Avoid special characters
# such as ';' or ')', which will create problems when parsing trees. The best
# is to only use letters, digits and '_'. Sequences can be interrupted by
# space and tab, but not by return characters. Be sure that the lengths of the
# sequences are all the same, and identical to the lengths indicated in the
# header. Sequences can be interleaved, in which case the taxon names may or
# may not be repeated in each block.


sub store_phylip {
    my $self    = shift;
    my $outfile = shift;
    my $args    = shift // {};          # HashRef (should not be empty...)

    my $short = $args->{short} // 1;
    my $clean = $args->{clean} // 0;
    my $chunk = $args->{chunk} // 60;

    open my $out, '>', $outfile;

    # print data matrix dimensions
    my $height = $self->count_seqs;
    my $width  = $self->width;
    say {$out} $height . q{ } . $width;

    # setup id format (this will also affect block-like structure)
    my $format = $short ? "%-10.10s %s\n" : "%s %s\n";
    my $method = $short ? 'full_id'       : 'foreign_id';
    my $sep    = $short ? q{ }            : q{};

    # optionally disable wrapping
    $chunk = $width if $chunk < 0;

    # optionally clean seq
    my @seqs = $self->all_seqs;
       @seqs = map { $_->clone->gapify('X') } @seqs if $clean;

    # output Ali in sequential or interleaved format
    for (my $site = 0; $site < $width; $site += $chunk) {