/Bio/SimpleAlign.pm
Perl | 3213 lines | 2534 code | 545 blank | 134 comment | 294 complexity | 4e6c8816b00e40960b17aaf4d793ee81 MD5 | raw file
Large files files are truncated, but you can click here to view the full file
- # BioPerl module for SimpleAlign
- #
- # Please direct questions and support issues to <bioperl-l@bioperl.org>
- #
- # Cared for by Heikki Lehvaslaiho <heikki-at-bioperl-dot-org>
- #
- # Copyright Ewan Birney
- #
- # You may distribute this module under the same terms as perl itself
- # POD documentation - main docs before the code
- #
- # History:
- # 11/3/00 Added threshold feature to consensus and consensus_aa - PS
- # May 2001 major rewrite - Heikki Lehvaslaiho
- =head1 NAME
- Bio::SimpleAlign - Multiple alignments held as a set of sequences
- =head1 SYNOPSIS
- # Use Bio::AlignIO to read in the alignment
- $str = Bio::AlignIO->new(-file => 't/data/testaln.pfam');
- $aln = $str->next_aln();
- # Describe
- print $aln->length;
- print $aln->num_residues;
- print $aln->is_flush;
- print $aln->num_sequences;
- print $aln->score;
- print $aln->percentage_identity;
- print $aln->consensus_string(50);
- # Find the position in the alignment for a sequence location
- $pos = $aln->column_from_residue_number('1433_LYCES', 14); # = 6;
- # Extract sequences and check values for the alignment column $pos
- foreach $seq ($aln->each_seq) {
- $res = $seq->subseq($pos, $pos);
- $count{$res}++;
- }
- foreach $res (keys %count) {
- printf "Res: %s Count: %2d\n", $res, $count{$res};
- }
- # Manipulate
- $aln->remove_seq($seq);
- $mini_aln = $aln->slice(20,30); # get a block of columns
- $mini_aln = $aln->select_noncont(1,3,5,7,11); # select certain sequences
- $new_aln = $aln->remove_columns([20,30]); # remove by position
- $new_aln = $aln->remove_columns(['mismatch']); # remove by property
- # Analyze
- $str = $aln->consensus_string($threshold_percent);
- $str = $aln->match_line();
- $str = $aln->cigar_line();
- $id = $aln->percentage_identity;
- # See the module documentation for details and more methods.
- =head1 DESCRIPTION
- SimpleAlign is an object that handles a multiple sequence alignment
- (MSA). It is very permissive of types (it does not insist on sequences
- being all same length, for example). Think of it as a set of sequences
- with a whole series of built-in manipulations and methods for reading and
- writing alignments.
- SimpleAlign uses L<Bio::LocatableSeq>, a subclass of L<Bio::PrimarySeq>,
- to store its sequences. These are subsequences with a start and end
- positions in the parent reference sequence. Each sequence in the
- SimpleAlign object is a Bio::LocatableSeq.
- SimpleAlign expects the combination of name, start, and end for a
- given sequence to be unique in the alignment, and this is the key for the
- internal hashes (name, start, end are abbreviated C<nse> in the code).
- However, in some cases people do not want the name/start-end to be displayed:
- either multiple names in an alignment or names specific to the alignment
- (ROA1_HUMAN_1, ROA1_HUMAN_2 etc). These names are called
- C<displayname>, and generally is what is used to print out the
- alignment. They default to name/start-end.
- The SimpleAlign Module is derived from the Align module by Ewan Birney.
- =head1 FEEDBACK
- =head2 Mailing Lists
- User feedback is an integral part of the evolution of this and other
- Bioperl modules. Send your comments and suggestions preferably to one
- of the Bioperl mailing lists. Your participation is much appreciated.
- bioperl-l@bioperl.org - General discussion
- http://bioperl.org/wiki/Mailing_lists - About the mailing lists
- =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:
- http://bugzilla.open-bio.org/
- =head1 AUTHOR
- Ewan Birney, birney@ebi.ac.uk
- =head1 CONTRIBUTORS
- Allen Day, allenday-at-ucla.edu,
- Richard Adams, Richard.Adams-at-ed.ac.uk,
- David J. Evans, David.Evans-at-vir.gla.ac.uk,
- Heikki Lehvaslaiho, heikki-at-bioperl-dot-org,
- Allen Smith, allens-at-cpan.org,
- Jason Stajich, jason-at-bioperl.org,
- Anthony Underwood, aunderwood-at-phls.org.uk,
- Xintao Wei & Giri Narasimhan, giri-at-cs.fiu.edu
- Brian Osborne, bosborne at alum.mit.edu
- Weigang Qiu, Weigang at GENECTR-HUNTER-CUNY-EDU
- Hongyu Zhang, forward at hongyu.org
- Jay Hannah, jay at jays.net
- Alexandr Bezginov, albezg at gmail.com
- =head1 SEE ALSO
- L<Bio::LocatableSeq>
- =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::SimpleAlign;
- use vars qw(%CONSERVATION_GROUPS);
- use strict;
- use Bio::LocatableSeq; # uses Seq's as list
- use Bio::Seq;
- use Bio::SeqFeature::Generic;
- BEGIN {
- # This data should probably be in a more centralized module...
- # it is taken from Clustalw documentation.
- # These are all the positively scoring groups that occur in the
- # Gonnet Pam250 matrix. The strong and weak groups are
- # defined as strong score >0.5 and weak score =<0.5 respectively.
- %CONSERVATION_GROUPS = (
- 'strong' => [ qw(
- STA
- NEQK
- NHQK
- NDEQ
- QHRK
- MILV
- MILF
- HY
- FYW )],
- 'weak' => [ qw(
- CSA
- ATV
- SAG
- STNK
- STPA
- SGND
- SNDEQK
- NDEQHK
- NEQHRK
- FVLIM
- HFY )],);
- }
- use base qw(Bio::Root::Root Bio::Align::AlignI Bio::AnnotatableI
- Bio::FeatureHolderI);
- =head2 new
- Title : new
- Usage : my $aln = Bio::SimpleAlign->new();
- Function : Creates a new simple align object
- Returns : Bio::SimpleAlign
- Args : -source => string representing the source program
- where this alignment came from
- -annotation => Bio::AnnotationCollectionI
- -seq_annotation => Bio::AnnotationCollectionI for sequences (requires -annotation also be set)
- -seqs => array ref containing Bio::LocatableSeq or Bio::Seq::Meta
- -consensus => consensus string
- -consensus_meta => Bio::Seq::Meta object containing consensus met information (kludge)
- =cut
- sub new {
- my($class,@args) = @_;
- my $self = $class->SUPER::new(@args);
- my ($src, $score, $id, $acc, $desc, $seqs, $feats, $coll, $sa, $con, $cmeta) = $self->_rearrange([qw(
- SOURCE
- SCORE
- ID
- ACCESSION
- DESCRIPTION
- SEQS
- FEATURES
- ANNOTATION
- SEQ_ANNOTATION
- CONSENSUS
- CONSENSUS_META
- )], @args);
- $src && $self->source($src);
- defined $score && $self->score($score);
- # we need to set up internal hashs first!
- $self->{'_seq'} = {};
- $self->{'_order'} = {};
- $self->{'_start_end_lists'} = {};
- $self->{'_dis_name'} = {};
- $self->{'_id'} = 'NoName';
- # maybe we should automatically read in from args. Hmmm...
- $id && $self->id($id);
- $acc && $self->accession($acc);
- $desc && $self->description($desc);
- $coll && $self->annotation($coll);
- # sequence annotation is layered into a provided annotation collection (or dies)
- if ($sa) {
- $self->throw("Must supply an alignment-based annotation collection (-annotation) ".
- "with a sequence annotation collection")
- if !$coll;
- $coll->add_Annotation('seq_annotation', $sa);
- }
- if ($feats && ref $feats eq 'ARRAY') {
- for my $feat (@$feats) {
- $self->add_SeqFeature($feat);
- }
- }
- $con && $self->consensus($con);
- $cmeta && $self->consensus_meta($cmeta);
- # assumes these are in correct alignment order
- if ($seqs && ref($seqs) eq 'ARRAY') {
- for my $seq (@$seqs) {
- $self->add_seq($seq);
- }
- }
- return $self; # success - we hope!
- }
- =head1 Modifier methods
- These methods modify the MSA by adding, removing or shuffling complete
- sequences.
- =head2 add_seq
- Title : add_seq
- Usage : $myalign->add_seq($newseq);
- $myalign->add_seq(-SEQ=>$newseq, -ORDER=>5);
- Function : Adds another sequence to the alignment. *Does not* align
- it - just adds it to the hashes.
- If -ORDER is specified, the sequence is inserted at the
- the position spec'd by -ORDER, and existing sequences
- are pushed down the storage array.
- Returns : nothing
- Args : A Bio::LocatableSeq object
- Positive integer for the sequence position (optional)
- See L<Bio::LocatableSeq> for more information
- =cut
- sub addSeq {
- my $self = shift;
- $self->deprecated("addSeq - deprecated method. Use add_seq() instead.");
- $self->add_seq(@_);
- }
- sub add_seq {
- my $self = shift;
- my @args = @_;
- my ($seq, $order) = $self->_rearrange([qw(SEQ ORDER)], @args);
- my ($name,$id,$start,$end);
- unless ($seq) {
- $self->throw("LocatableSeq argument required");
- }
- if( ! ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
- $self->throw("Unable to process non locatable sequences [". ref($seq). "]");
- }
- !defined($order) and $order = 1 + keys %{$self->{'_seq'}}; # default
- $order--; # jay's patch (user-specified order is 1-origin)
-
- if ($order < 0) {
- $self->throw("User-specified value for ORDER must be >= 1");
- }
- $id = $seq->id() ||$seq->display_id || $seq->primary_id;
- # build the symbol list for this sequence,
- # will prune out the gap and missing/match chars
- # when actually asked for the symbol list in the
- # symbol_chars
- # map { $self->{'_symbols'}->{$_} = 1; } split(//,$seq->seq) if $seq->seq;
- $name = $seq->get_nse;
- if( $self->{'_seq'}->{$name} ) {
- $self->warn("Replacing one sequence [$name]\n") unless $self->verbose < 0;
- }
- else {
- $self->debug( "Assigning $name to $order\n");
- my $ordh = $self->{'_order'};
- if ($ordh->{$order}) {
- # make space to insert
- # $c->() returns (in reverse order) the first subsequence
- # of consecutive integers; i.e., $c->(1,2,3,5,6,7) returns
- # (3,2,1), and $c->(2,4,5) returns (2).
- my $c;
- $c = sub { return (($_[1]-$_[0] == 1) ? ($c->(@_[1..$#_]),$_[0]) : $_[0]); };
- map {
- $ordh->{$_+1} = $ordh->{$_}
- } $c->(sort {$a <=> $b} grep {$_ >= $order} keys %{$ordh});
- }
- $ordh->{$order} = $name;
- unless( exists( $self->{'_start_end_lists'}->{$id})) {
- $self->{'_start_end_lists'}->{$id} = [];
- }
- push @{$self->{'_start_end_lists'}->{$id}}, $seq;
- }
- $self->{'_seq'}->{$name} = $seq;
- }
- =head2 remove_seq
- Title : remove_seq
- Usage : $aln->remove_seq($seq);
- Function : Removes a single sequence from an alignment
- Returns :
- Argument : a Bio::LocatableSeq object
- =cut
- sub removeSeq {
- my $self = shift;
- $self->deprecated("removeSeq - deprecated method. Use remove_seq() instead.");
- $self->remove_seq(@_);
- }
- sub remove_seq {
- my $self = shift;
- my $seq = shift;
- my ($name,$id,$start,$end);
- $self->throw("Need Bio::Locatable seq argument ")
- unless ref $seq && $seq->isa( 'Bio::LocatableSeq');
- $id = $seq->id();
- $start = $seq->start();
- $end = $seq->end();
- $name = sprintf("%s/%d-%d",$id,$start,$end);
- if( !exists $self->{'_seq'}->{$name} ) {
- $self->throw("Sequence $name does not exist in the alignment to remove!");
- }
- delete $self->{'_seq'}->{$name};
- # we need to remove this seq from the start_end_lists hash
- if (exists $self->{'_start_end_lists'}->{$id}) {
- # we need to find the sequence in the array.
- my ($i, $found);;
- for ($i=0; $i < @{$self->{'_start_end_lists'}->{$id}}; $i++) {
- if (${$self->{'_start_end_lists'}->{$id}}[$i] eq $seq) {
- $found = 1;
- last;
- }
- }
- if ($found) {
- splice @{$self->{'_start_end_lists'}->{$id}}, $i, 1;
- }
- else {
- $self->throw("Could not find the sequence to remoce from the start-end list");
- }
- }
- else {
- $self->throw("There is no seq list for the name $id");
- }
- # we need to shift order hash
- my %rev_order = reverse %{$self->{'_order'}};
- my $no = $rev_order{$name};
- my $num_sequences = $self->num_sequences;
- for (; $no < $num_sequences; $no++) {
- $self->{'_order'}->{$no} = $self->{'_order'}->{$no+1};
- }
- delete $self->{'_order'}->{$no};
- return 1;
- }
- =head2 purge
- Title : purge
- Usage : $aln->purge(0.7);
- Function: Removes sequences above given sequence similarity
- This function will grind on large alignments. Beware!
- Example :
- Returns : An array of the removed sequences
- Args : float, threshold for similarity
- =cut
- sub purge {
- my ($self,$perc) = @_;
- my (%duplicate, @dups);
- my @seqs = $self->each_seq();
- for (my $i=0;$i< @seqs - 1;$i++ ) { #for each seq in alignment
- my $seq = $seqs[$i];
- #skip if already in duplicate hash
- next if exists $duplicate{$seq->display_id} ;
- my $one = $seq->seq();
- my @one = split '', $one; #split to get 1aa per array element
- for (my $j=$i+1;$j < @seqs;$j++) {
- my $seq2 = $seqs[$j];
- #skip if already in duplicate hash
- next if exists $duplicate{$seq2->display_id} ;
- my $two = $seq2->seq();
- my @two = split '', $two;
- my $count = 0;
- my $res = 0;
- for (my $k=0;$k<@one;$k++) {
- if ( $one[$k] ne '.' && $one[$k] ne '-' && defined($two[$k]) &&
- $one[$k] eq $two[$k]) {
- $count++;
- }
- if ( $one[$k] ne '.' && $one[$k] ne '-' && defined($two[$k]) &&
- $two[$k] ne '.' && $two[$k] ne '-' ) {
- $res++;
- }
- }
- my $ratio = 0;
- $ratio = $count/$res unless $res == 0;
- # if above threshold put in duplicate hash and push onto
- # duplicate array for returning to get_unique
- if ( $ratio > $perc ) {
- $self->warn("duplicate: ", $seq2->display_id) if $self->verbose > 0;
- $duplicate{$seq2->display_id} = 1;
- push @dups, $seq2;
- }
- }
- }
- foreach my $seq (@dups) {
- $self->remove_seq($seq);
- }
- return @dups;
- }
- =head2 sort_alphabetically
- Title : sort_alphabetically
- Usage : $ali->sort_alphabetically
- Function : Changes the order of the alignment to alphabetical on name
- followed by numerical by number.
- Returns :
- Argument :
- =cut
- sub sort_alphabetically {
- my $self = shift;
- my ($seq,$nse,@arr,%hash,$count);
- foreach $seq ( $self->each_seq() ) {
- $nse = $seq->get_nse;
- $hash{$nse} = $seq;
- }
- $count = 0;
- %{$self->{'_order'}} = (); # reset the hash;
- foreach $nse ( sort _alpha_startend keys %hash) {
- $self->{'_order'}->{$count} = $nse;
- $count++;
- }
- 1;
- }
- =head2 sort_by_list
- Title : sort_by_list
- Usage : $aln_ordered=$aln->sort_by_list($list_file)
- Function : Arbitrarily order sequences in an alignment
- Returns : A new Bio::SimpleAlign object
- Argument : a file listing sequence names in intended order (one name per line)
- =cut
- sub sort_by_list {
- my ($self, $list) = @_;
- my (@seq, @ids, %order);
- foreach my $seq ( $self->each_seq() ) {
- push @seq, $seq;
- push @ids, $seq->display_id;
- }
- my $ct=1;
- open(my $listfh, '<', $list) || $self->throw("can't open file for reading: $list");
- while (<$listfh>) {
- chomp;
- my $name=$_;
- $self->throw("Not found in alignment: $name") unless &_in_aln($name, \@ids);
- $order{$name}=$ct++;
- }
- close($listfh);
-
- # use the map-sort-map idiom:
- my @sorted= map { $_->[1] } sort { $a->[0] <=> $b->[0] } map { [$order{$_->id()}, $_] } @seq;
- my $aln = $self->new;
- foreach (@sorted) { $aln->add_seq($_) }
- return $aln;
- }
- =head2 set_new_reference
- Title : set_new_reference
- Usage : $aln->set_new_reference(3 or 'B31'): Select the 3rd sequence, or
- the sequence whoes name is "B31" (full, exact, and case-sensitive),
- as the reference (1st) sequence
- Function : Change/Set a new reference (i.e., the first) sequence
- Returns : a new Bio::SimpleAlign object.
- Throws an exception if designated sequence not found
- Argument : a positive integer of sequence order, or a sequence name
- in the original alignment
- =cut
- sub set_new_reference {
- my ($self, $seqid) = @_;
- my $aln = $self->new;
- my (@seq, @ids, @new_seq);
- my $is_num=0;
- foreach my $seq ( $self->each_seq() ) {
- push @seq, $seq;
- push @ids, $seq->display_id;
- }
- if ($seqid =~ /^\d+$/) { # argument is seq position
- $is_num=1;
- $self->throw("The new reference sequence number has to be a positive integer >1 and <= num_sequences ") if ($seqid <= 1 || $seqid > $self->num_sequences);
- } else { # argument is a seq name
- $self->throw("The new reference sequence not in alignment ") unless &_in_aln($seqid, \@ids);
- }
- for (my $i=0; $i<=$#seq; $i++) {
- my $pos=$i+1;
- if ( ($is_num && $pos == $seqid) || ($seqid eq $seq[$i]->display_id) ) {
- unshift @new_seq, $seq[$i];
- } else {
- push @new_seq, $seq[$i];
- }
- }
- foreach (@new_seq) { $aln->add_seq($_); }
- return $aln;
- }
- sub _in_aln { # check if input name exists in the alignment
- my ($str, $ref) = @_;
- foreach (@$ref) {
- return 1 if $str eq $_;
- }
- return 0;
- }
- =head2 uniq_seq
- Title : uniq_seq
- Usage : $aln->uniq_seq(): Remove identical sequences in
- in the alignment. Ambiguous base ("N", "n") and
- leading and ending gaps ("-") are NOT counted as
- differences.
- Function : Make a new alignment of unique sequence types (STs)
- Returns : 1a. if called in a scalar context,
- a new Bio::SimpleAlign object (all sequences renamed as "ST")
- 1b. if called in an array context,
- a new Bio::SimpleAlign object, and a hashref whose keys
- are sequence types, and whose values are arrayrefs to
- lists of sequence ids within the corresponding sequence type
- 2. if $aln->verbose > 0, ST of each sequence is sent to
- STDERR (in a tabular format)
- Argument : None
- =cut
- sub uniq_seq {
- my ($self, $seqid) = @_;
- my $aln = $self->new;
- my (%member, %order, @seq, @uniq_str, $st);
- my $order=0;
- my $len = $self->length();
- $st = {};
- foreach my $seq ( $self->each_seq() ) {
- my $str = $seq->seq();
- # it's necessary to ignore "n", "N", leading gaps and ending gaps in
- # comparing two sequence strings
- # 1st, convert "n", "N" to "?" (for DNA sequence only):
- $str =~ s/n/\?/gi if $str =~ /^[atcgn-]+$/i;
- # 2nd, convert leading and ending gaps to "?":
- $str = &_convert_leading_ending_gaps($str, '-', '?');
- # Note that '?' also can mean unknown residue.
- # I don't like making global class member changes like this, too
- # prone to errors... -- cjfields 08-11-18
- local $Bio::LocatableSeq::GAP_SYMBOLS = '-\?';
- my $new = Bio::LocatableSeq->new(
- -id => $seq->id(),
- -alphabet=> $seq->alphabet,
- -seq => $str,
- -start => $seq->start,
- -end => $seq->end
- );
- push @seq, $new;
- }
- foreach my $seq (@seq) {
- my $str = $seq->seq();
- my ($seen, $key) = &_check_uniq($str, \@uniq_str, $len);
- if ($seen) { # seen before
- my @memb = @{$member{$key}};
- push @memb, $seq;
- $member{$key} = \@memb;
- } else { # not seen
- push @uniq_str, $key;
- $order++;
- $member{$key} = [ ($seq) ];
- $order{$key} = $order;
- }
- }
- foreach my $str (sort {$order{$a} <=> $order{$b}} keys %order) { # sort by input order
- # convert leading/ending "?" back into "-" ("?" throws errors by SimpleAlign):
- my $str2 = &_convert_leading_ending_gaps($str, '?', '-');
- # convert middle "?" back into "N" ("?" throws errors by SimpleAlign):
- $str2 =~ s/\?/N/g if $str2 =~ /^[atcg\-\?]+$/i;
- my $gap='-';
- my $end= CORE::length($str2);
- $end -= CORE::length($1) while $str2 =~ m/($gap+)/g;
- my $new = Bio::LocatableSeq->new(-id =>"ST".$order{$str},
- -seq =>$str2,
- -start=>1,
- -end =>$end
- );
- $aln->add_seq($new);
- foreach (@{$member{$str}}) {
- push @{$$st{$order{$str}}}, $_->id(); # per Tristan's patch/Bug #2805
- $self->debug($_->id(), "\t", "ST", $order{$str}, "\n");
- }
- }
- return wantarray ? ($aln, $st) : $aln;
- }
- sub _check_uniq { # check if same seq exists in the alignment
- my ($str1, $ref, $length) = @_;
- my @char1=split //, $str1;
- my @array=@$ref;
- return (0, $str1) if @array==0; # not seen (1st sequence)
- foreach my $str2 (@array) {
- my $diff=0;
- my @char2=split //, $str2;
- for (my $i=0; $i<=$length-1; $i++) {
- next if $char1[$i] eq '?';
- next if $char2[$i] eq '?';
- $diff++ if $char1[$i] ne $char2[$i];
- }
- return (1, $str2) if $diff == 0; # seen before
- }
- return (0, $str1); # not seen
- }
- sub _convert_leading_ending_gaps {
- my $s=shift;
- my $sym1=shift;
- my $sym2=shift;
- my @array=split //, $s;
- # convert leading char:
- for (my $i=0; $i<=$#array; $i++) {
- ($array[$i] eq $sym1) ? ($array[$i] = $sym2):(last);
- }
- # convert ending char:
- for (my $i = $#array; $i>= 0; $i--) {
- ($array[$i] eq $sym1) ? ($array[$i] = $sym2):(last);
- }
- my $s_new=join '', @array;
- return $s_new;
- }
- =head1 Sequence selection methods
- Methods returning one or more sequences objects.
- =head2 each_seq
- Title : each_seq
- Usage : foreach $seq ( $align->each_seq() )
- Function : Gets a Seq object from the alignment
- Returns : Seq object
- Argument :
- =cut
- sub eachSeq {
- my $self = shift;
- $self->deprecated("eachSeq - deprecated method. Use each_seq() instead.");
- $self->each_seq();
- }
- sub each_seq {
- my $self = shift;
- my (@arr,$order);
- foreach $order ( sort { $a <=> $b } keys %{$self->{'_order'}} ) {
- if( exists $self->{'_seq'}->{$self->{'_order'}->{$order}} ) {
- push(@arr,$self->{'_seq'}->{$self->{'_order'}->{$order}});
- }
- }
- return @arr;
- }
- =head2 each_alphabetically
- Title : each_alphabetically
- Usage : foreach $seq ( $ali->each_alphabetically() )
- Function : Returns a sequence object, but the objects are returned
- in alphabetically sorted order.
- Does not change the order of the alignment.
- Returns : Seq object
- Argument :
- =cut
- sub each_alphabetically {
- my $self = shift;
- my ($seq,$nse,@arr,%hash,$count);
- foreach $seq ( $self->each_seq() ) {
- $nse = $seq->get_nse;
- $hash{$nse} = $seq;
- }
- foreach $nse ( sort _alpha_startend keys %hash) {
- push(@arr,$hash{$nse});
- }
- return @arr;
- }
- sub _alpha_startend {
- my ($aname,$astart,$bname,$bstart);
- ($aname,$astart) = split (/-/,$a);
- ($bname,$bstart) = split (/-/,$b);
- if( $aname eq $bname ) {
- return $astart <=> $bstart;
- }
- else {
- return $aname cmp $bname;
- }
- }
- =head2 each_seq_with_id
- Title : each_seq_with_id
- Usage : foreach $seq ( $align->each_seq_with_id() )
- Function : Gets a Seq objects from the alignment, the contents
- being those sequences with the given name (there may be
- more than one)
- Returns : Seq object
- Argument : a seq name
- =cut
- sub eachSeqWithId {
- my $self = shift;
- $self->deprecated("eachSeqWithId - deprecated method. Use each_seq_with_id() instead.");
- $self->each_seq_with_id(@_);
- }
- sub each_seq_with_id {
- my $self = shift;
- my $id = shift;
- $self->throw("Method each_seq_with_id needs a sequence name argument")
- unless defined $id;
- my (@arr, $seq);
- if (exists($self->{'_start_end_lists'}->{$id})) {
- @arr = @{$self->{'_start_end_lists'}->{$id}};
- }
- return @arr;
- }
- =head2 get_seq_by_pos
- Title : get_seq_by_pos
- Usage : $seq = $aln->get_seq_by_pos(3) # third sequence from the alignment
- Function : Gets a sequence based on its position in the alignment.
- Numbering starts from 1. Sequence positions larger than
- num_sequences() will thow an error.
- Returns : a Bio::LocatableSeq object
- Args : positive integer for the sequence position
- =cut
- sub get_seq_by_pos {
- my $self = shift;
- my ($pos) = @_;
- $self->throw("Sequence position has to be a positive integer, not [$pos]")
- unless $pos =~ /^\d+$/ and $pos > 0;
- $self->throw("No sequence at position [$pos]")
- unless $pos <= $self->num_sequences ;
- my $nse = $self->{'_order'}->{--$pos};
- return $self->{'_seq'}->{$nse};
- }
- =head2 get_seq_by_id
- Title : get_seq_by_id
- Usage : $seq = $aln->get_seq_by_id($name) # seq named $name
- Function : Gets a sequence based on its name.
- Sequences that do not exist will warn and return undef
- Returns : a Bio::LocatableSeq object
- Args : string for sequence name
- =cut
- sub get_seq_by_id {
- my ($self,$name) = @_;
- unless( defined $name ) {
- $self->warn("Must provide a sequence name");
- return;
- }
- for my $seq ( values %{$self->{'_seq'}} ) {
- if ( $seq->id eq $name) {
- return $seq;
- }
- }
- return;
- }
- =head2 seq_with_features
- Title : seq_with_features
- Usage : $seq = $aln->seq_with_features(-pos => 1,
- -consensus => 60
- -mask =>
- sub { my $consensus = shift;
- for my $i (1..5){
- my $n = 'N' x $i;
- my $q = '\?' x $i;
- while($consensus =~ /[^?]$q[^?]/){
- $consensus =~ s/([^?])$q([^?])/$1$n$2/;
- }
- }
- return $consensus;
- }
- );
- Function: produces a Bio::Seq object by first splicing gaps from -pos
- (by means of a splice_by_seq_pos() call), then creating
- features using non-? chars (by means of a consensus_string()
- call with stringency -consensus).
- Returns : a Bio::Seq object
- Args : -pos : required. sequence from which to build the Bio::Seq
- object
- -consensus : optional, defaults to consensus_string()'s
- default cutoff value
- -mask : optional, a coderef to apply to consensus_string()'s
- output before building features. this may be useful for
- closing gaps of 1 bp by masking over them with N, for
- instance
- =cut
- sub seq_with_features{
- my ($self,%arg) = @_;
- #first do the preparatory splice
- $self->throw("must provide a -pos argument") unless $arg{-pos};
- $self->splice_by_seq_pos($arg{-pos});
- my $consensus_string = $self->consensus_string($arg{-consensus});
- $consensus_string = $arg{-mask}->($consensus_string)
- if defined($arg{-mask});
- my(@bs,@es);
- push @bs, 1 if $consensus_string =~ /^[^?]/;
- while($consensus_string =~ /\?[^?]/g){
- push @bs, pos($consensus_string);
- }
- while($consensus_string =~ /[^?]\?/g){
- push @es, pos($consensus_string);
- }
- push @es, CORE::length($consensus_string) if $consensus_string =~ /[^?]$/;
- my $seq = Bio::Seq->new();
- # my $rootfeature = Bio::SeqFeature::Generic->new(
- # -source_tag => 'location',
- # -start => $self->get_seq_by_pos($arg{-pos})->start,
- # -end => $self->get_seq_by_pos($arg{-pos})->end,
- # );
- # $seq->add_SeqFeature($rootfeature);
- while(my $b = shift @bs){
- my $e = shift @es;
- $seq->add_SeqFeature(
- Bio::SeqFeature::Generic->new(
- -start => $b - 1 + $self->get_seq_by_pos($arg{-pos})->start,
- -end => $e - 1 + $self->get_seq_by_pos($arg{-pos})->start,
- -source_tag => $self->source || 'MSA',
- )
- );
- }
- return $seq;
- }
- =head1 Create new alignments
- The result of these methods are horizontal or vertical subsets of the
- current MSA.
- =head2 select
- Title : select
- Usage : $aln2 = $aln->select(1, 3) # three first sequences
- Function : Creates a new alignment from a continuous subset of
- sequences. Numbering starts from 1. Sequence positions
- larger than num_sequences() will thow an error.
- Returns : a Bio::SimpleAlign object
- Args : positive integer for the first sequence
- positive integer for the last sequence to include (optional)
- =cut
- sub select {
- my $self = shift;
- my ($start, $end) = @_;
- $self->throw("Select start has to be a positive integer, not [$start]")
- unless $start =~ /^\d+$/ and $start > 0;
- $self->throw("Select end has to be a positive integer, not [$end]")
- unless $end =~ /^\d+$/ and $end > 0;
- $self->throw("Select $start [$start] has to be smaller than or equal to end [$end]")
- unless $start <= $end;
- my $aln = $self->new;
- foreach my $pos ($start .. $end) {
- $aln->add_seq($self->get_seq_by_pos($pos));
- }
- $aln->id($self->id);
- # fix for meta, sf, ann
- return $aln;
- }
- =head2 select_noncont
- Title : select_noncont
- Usage : # 1st and 3rd sequences, sorted
- $aln2 = $aln->select_noncont(1, 3)
- # 1st and 3rd sequences, sorted (same as first)
- $aln2 = $aln->select_noncont(3, 1)
- # 1st and 3rd sequences, unsorted
- $aln2 = $aln->select_noncont('nosort',3, 1)
- Function : Creates a new alignment from a subset of sequences. Numbering
- starts from 1. Sequence positions larger than num_sequences() will
- throw an error. Sorts the order added to new alignment by default,
- to prevent sorting pass 'nosort' as the first argument in the list.
- Returns : a Bio::SimpleAlign object
- Args : array of integers for the sequences. If the string 'nosort' is
- passed as the first argument, the sequences will not be sorted
- in the new alignment but will appear in the order listed.
- =cut
- sub select_noncont {
- my $self = shift;
- my $nosort = 0;
- my (@pos) = @_;
- if ($pos[0] !~ m{^\d+$}) {
- my $sortcmd = shift @pos;
- if ($sortcmd eq 'nosort') {
- $nosort = 1;
- } else {
- $self->throw("Command not recognized: $sortcmd. Only 'nosort' implemented at this time.");
- }
- }
-
- my $end = $self->num_sequences;
- foreach ( @pos ) {
- $self->throw("position must be a positive integer, > 0 and <= $end not [$_]")
- unless( /^\d+$/ && $_ > 0 && $_ <= $end );
- }
-
- @pos = sort {$a <=> $b} @pos unless $nosort;
-
- my $aln = $self->new;
- foreach my $p (@pos) {
- $aln->add_seq($self->get_seq_by_pos($p));
- }
- $aln->id($self->id);
- # fix for meta, sf, ann
- return $aln;
- }
- =head2 slice
- Title : slice
- Usage : $aln2 = $aln->slice(20,30)
- Function : Creates a slice from the alignment inclusive of start and
- end columns, and the first column in the alignment is denoted 1.
- Sequences with no residues in the slice are excluded from the
- new alignment and a warning is printed. Slice beyond the length of
- the sequence does not do padding.
- Returns : A Bio::SimpleAlign object
- Args : Positive integer for start column, positive integer for end column,
- optional boolean which if true will keep gap-only columns in the newly
- created slice. Example:
- $aln2 = $aln->slice(20,30,1)
- =cut
- sub slice {
- my $self = shift;
- my ($start, $end, $keep_gap_only) = @_;
- $self->throw("Slice start has to be a positive integer, not [$start]")
- unless $start =~ /^\d+$/ and $start > 0;
- $self->throw("Slice end has to be a positive integer, not [$end]")
- unless $end =~ /^\d+$/ and $end > 0;
- $self->throw("Slice start [$start] has to be smaller than or equal to end [$end]")
- unless $start <= $end;
- $self->throw("This alignment has only ". $self->length . " residues. Slice start " .
- "[$start] is too big.") if $start > $self->length;
- my $cons_meta = $self->consensus_meta;
- my $aln = $self->new;
- $aln->id($self->id);
- foreach my $seq ( $self->each_seq() ) {
- my $new_seq = $seq->isa('Bio::Seq::MetaI') ?
- Bio::Seq::Meta->new
- (-id => $seq->id,
- -alphabet => $seq->alphabet,
- -strand => $seq->strand,
- -verbose => $self->verbose) :
- Bio::LocatableSeq->new
- (-id => $seq->id,
- -alphabet => $seq->alphabet,
- -strand => $seq->strand,
- -verbose => $self->verbose);
-
- # seq
- my $seq_end = $end;
- $seq_end = $seq->length if( $end > $seq->length );
- my $slice_seq = $seq->subseq($start, $seq_end);
- $new_seq->seq( $slice_seq );
- $slice_seq =~ s/\W//g;
- if ($start > 1) {
- my $pre_start_seq = $seq->subseq(1, $start - 1);
- $pre_start_seq =~ s/\W//g;
- if (!defined($seq->strand)) {
- $new_seq->start( $seq->start + CORE::length($pre_start_seq) );
- } elsif ($seq->strand < 0){
- $new_seq->start( $seq->end - CORE::length($pre_start_seq) - CORE::length($slice_seq) + 1);
- } else {
- $new_seq->start( $seq->start + CORE::length($pre_start_seq) );
- }
- } else {
- if ((defined $seq->strand)&&($seq->strand < 0)){
- $new_seq->start( $seq->end - CORE::length($slice_seq) + 1);
- } else {
- $new_seq->start( $seq->start);
- }
- }
- if ($new_seq->isa('Bio::Seq::MetaI')) {
- for my $meta_name ($seq->meta_names) {
- $new_seq->named_meta($meta_name, $seq->named_submeta($meta_name, $start, $end));
- }
- }
- $new_seq->end( $new_seq->start + CORE::length($slice_seq) - 1 );
- if ($new_seq->start and $new_seq->end >= $new_seq->start) {
- $aln->add_seq($new_seq);
- } else {
- if( $keep_gap_only ) {
- $aln->add_seq($new_seq);
- } else {
- my $nse = $seq->get_nse();
- $self->warn("Slice [$start-$end] of sequence [$nse] contains no residues.".
- " Sequence excluded from the new alignment.");
- }
- }
- }
- if ($cons_meta) {
- my $new = Bio::Seq::Meta->new();
- for my $meta_name ($cons_meta->meta_names) {
- $new->named_meta($meta_name, $cons_meta->named_submeta($meta_name, $start, $end));
- }
- $aln->consensus_meta($new);
- }
- $aln->annotation($self->annotation);
- # fix for meta, sf, ann
- return $aln;
- }
- =head2 remove_columns
- Title : remove_columns
- Usage : $aln2 = $aln->remove_columns(['mismatch','weak']) or
- $aln2 = $aln->remove_columns([0,0],[6,8])
- Function : Creates an aligment with columns removed corresponding to
- the specified type or by specifying the columns by number.
- Returns : Bio::SimpleAlign object
- Args : Array ref of types ('match'|'weak'|'strong'|'mismatch'|'gaps'|
- 'all_gaps_columns') or array ref where the referenced array
- contains a pair of integers that specify a range.
- The first column is 0
- =cut
- sub remove_columns {
- my ($self,@args) = @_;
- @args || $self->throw("Must supply column ranges or column types");
- my $aln;
- if ($args[0][0] =~ /^[a-z_]+$/i) {
- $aln = $self->_remove_columns_by_type($args[0]);
- } elsif ($args[0][0] =~ /^\d+$/) {
- $aln = $self->_remove_columns_by_num(\@args);
- } else {
- $self->throw("You must pass array references to remove_columns(), not @args");
- }
- # fix for meta, sf, ann
- $aln;
- }
- =head2 remove_gaps
- Title : remove_gaps
- Usage : $aln2 = $aln->remove_gaps
- Function : Creates an aligment with gaps removed
- Returns : a Bio::SimpleAlign object
- Args : a gap character(optional) if none specified taken
- from $self->gap_char,
- [optional] $all_gaps_columns flag (1 or 0, default is 0)
- indicates that only all-gaps columns should be deleted
- Used from method L<remove_columns> in most cases. Set gap character
- using L<gap_char()|gap_char>.
- =cut
- sub remove_gaps {
- my ($self,$gapchar,$all_gaps_columns) = @_;
- my $gap_line;
- if ($all_gaps_columns) {
- $gap_line = $self->all_gap_line($gapchar);
- } else {
- $gap_line = $self->gap_line($gapchar);
- }
- my $aln = $self->new;
- my @remove;
- my $length = 0;
- my $del_char = $gapchar || $self->gap_char;
- # Do the matching to get the segments to remove
- while ($gap_line =~ m/[$del_char]/g) {
- my $start = pos($gap_line)-1;
- $gap_line=~/\G[$del_char]+/gc;
- my $end = pos($gap_line)-1;
- #have to offset the start and end for subsequent removes
- $start-=$length;
- $end -=$length;
- $length += ($end-$start+1);
- push @remove, [$start,$end];
- }
- #remove the segments
- $aln = $#remove >= 0 ? $self->_remove_col($aln,\@remove) : $self;
- # fix for meta, sf, ann
- return $aln;
- }
- sub _remove_col {
- my ($self,$aln,$remove) = @_;
- my @new;
-
- my $gap = $self->gap_char;
-
- # splice out the segments and create new seq
- foreach my $seq($self->each_seq){
- my $new_seq = Bio::LocatableSeq->new(
- -id => $seq->id,
- -alphabet=> $seq->alphabet,
- -strand => $seq->strand,
- -verbose => $self->verbose);
- my $sequence = $seq->seq;
- foreach my $pair(@{$remove}){
- my $start = $pair->[0];
- my $end = $pair->[1];
- $sequence = $seq->seq unless $sequence;
- my $orig = $sequence;
- my $head = $start > 0 ? substr($sequence, 0, $start) : '';
- my $tail = ($end + 1) >= CORE::length($sequence) ? '' : substr($sequence, $end + 1);
- $sequence = $head.$tail;
- # start
- unless (defined $new_seq->start) {
- if ($start == 0) {
- my $start_adjust = () = substr($orig, 0, $end + 1) =~ /$gap/g;
- $new_seq->start($seq->start + $end + 1 - $start_adjust);
- }
- else {
- my $start_adjust = $orig =~ /^$gap+/;
- if ($start_adjust) {
- $start_adjust = $+[0] == $start;
- }
- $new_seq->start($seq->start + $start_adjust);
- }
- }
- # end
- if (($end + 1) >= CORE::length($orig)) {
- my $end_adjust = () = substr($orig, $start) =~ /$gap/g;
- $new_seq->end($seq->end - (CORE::length($orig) - $start) + $end_adjust);
- }
- else {
- $new_seq->end($seq->end);
- }
- }
-
- if ($new_seq->end < $new_seq->start) {
- # we removed all columns except for gaps: set to 0 to indicate no
- # sequence
- $new_seq->start(0);
- $new_seq->end(0);
- }
-
- $new_seq->seq($sequence) if $sequence;
- push @new, $new_seq;
- }
- # add the new seqs to the alignment
- foreach my $new(@new){
- $aln->add_seq($new);
- }
- # fix for meta, sf, ann
- return $aln;
- }
- sub _remove_columns_by_type {
- my ($self,$type) = @_;
- my $aln = $self->new;
- my @remove;
- my $gap = $self->gap_char if (grep { $_ eq 'gaps'} @{$type});
- my $all_gaps_columns = $self->gap_char if (grep /all_gaps_columns/,@{$type});
- my %matchchars = ( 'match' => '\*',
- 'weak' => '\.',
- 'strong' => ':',
- 'mismatch' => ' ',
- 'gaps' => '',
- 'all_gaps_columns' => ''
- );
- # get the characters to delete against
- my $del_char;
- foreach my $type (@{$type}){
- $del_char.= $matchchars{$type};
- }
- my $length = 0;
- my $match_line = $self->match_line;
- # do the matching to get the segments to remove
- if($del_char){
- while($match_line =~ m/[$del_char]/g ){
- my $start = pos($match_line)-1;
- $match_line=~/\G[$del_char]+/gc;
- my $end = pos($match_line)-1;
- #have to offset the start and end for subsequent removes
- $start-=$length;
- $end -=$length;
- $length += ($end-$start+1);
- push @remove, [$start,$end];
- }
- }
- # remove the segments
- $aln = $#remove >= 0 ? $self->_remove_col($aln,\@remove) : $self;
- $aln = $aln->remove_gaps() if $gap;
- $aln = $aln->remove_gaps('', 1) if $all_gaps_columns;
- # fix for meta, sf, ann
- $aln;
- }
- sub _remove_columns_by_num {
- my ($self,$positions) = @_;
- my $aln = $self->new;
- # sort the positions
- @$positions = sort { $a->[0] <=> $b->[0] } @$positions;
-
- my @remove;
- my $length = 0;
- foreach my $pos (@{$positions}) {
- my ($start, $end) = @{$pos};
-
- #have to offset the start and end for subsequent removes
- $start-=$length;
- $end -=$length;
- $length += ($end-$start+1);
- push @remove, [$start,$end];
- }
- #remove the segments
- $aln = $#remove >= 0 ? $self->_remove_col($aln,\@remove) : $self;
- # fix for meta, sf, ann
- $aln;
- }
- =head1 Change sequences within the MSA
- These methods affect characters in all sequences without changing the
- alignment.
- =head2 splice_by_seq_pos
- Title : splice_by_seq_pos
- Usage : $status = splice_by_seq_pos(1);
- Function: splices all aligned sequences where the specified sequence
- has gaps.
- Example :
- Returns : 1 on success
- Args : position of sequence to splice by
- =cut
- sub splice_by_seq_pos{
- my ($self,$pos) = @_;
- my $guide = $self->get_seq_by_pos($pos);
- my $guide_seq = $guide->seq;
- $guide_seq =~ s/\./\-/g;
- my @gaps = ();
- $pos = -1;
- while(($pos = index($guide_seq, '-', $pos)) > -1 ){
- unshift @gaps, $pos;
- $pos++;
- }
- foreach my $seq ($self->each_seq){
- my @bases = split '', $seq->seq;
- splice(@bases, $_, 1) foreach @gaps;
- $seq->seq(join('', @bases));
- }
- 1;
- }
- =head2 map_chars
- Title : map_chars
- Usage : $ali->map_chars('\.','-')
- Function : Does a s/$arg1/$arg2/ on the sequences. Useful for gap
- characters
- Notice that the from (arg1) is interpretted as a regex,
- so be careful about quoting meta characters (eg
- $ali->map_chars('.','-') wont do what you want)
- Returns :
- Argument : 'from' rexexp
- 'to' string
- =cut
- sub map_chars {
- my $self = shift;
- my $from = shift;
- my $to = shift;
- my ($seq,$temp);
- $self->throw("Need exactly two arguments")
- unless defined $from and defined $to;
- foreach $seq ( $self->each_seq() ) {
- $temp = $seq->seq();
- $temp =~ s/$from/$to/g;
- $seq->seq($temp);
- }
- return 1;
- }
- =head2 uppercase
- Title : uppercase()
- Usage : $ali->uppercase()
- Function : Sets all the sequences to uppercase
- Returns :
- Argument :
- =cut
- sub uppercase {
- my $self = shift;
- my $seq;
- my $temp;
- foreach $seq ( $self->each_seq() ) {
- $temp = $seq->seq();
- $temp =~ tr/[a-z]/[A-Z]/;
- $seq->seq($temp);
- }
- return 1;
- }
- =head2 cigar_line
- Title : cigar_line()
- Usage : %cigars = $align->cigar_line()
- Function : Generates a "cigar" (Compact Idiosyncratic Gapped Alignment
- Report) line for each sequence in the alignment. Examples are
- "1,60" or "5,10:12,58", where the numbers refer to conserved
- positions within the alignment. The keys of the hash are the
- NSEs (name/start/end) assigned to each sequence.
- Args : threshold (optional, defaults to 100)
- Returns : Hash of strings (cigar lines)
- =cut
- sub cigar_line {
- my $self = shift;
- my $thr=shift||100;
- my %cigars;
- my @consensus = split "",($self->consensus_string($thr));
- my $len = $self->length;
- my $gapchar = $self->gap_char;
- # create a precursor, something like (1,4,5,6,7,33,45),
- # where each number corresponds to a conserved position
- foreach my $seq ( $self->each_seq ) {
- my @seq = split "", uc ($seq->seq);
- my $pos = 1;
- for (my $x = 0 ; $x < $len ; $x++ ) {
- if ($seq[$x] eq $consensus[$x]) {
- push @{$cigars{$seq->get_nse}},$pos;
- $pos++;
- } elsif ($seq[$x] ne $gapchar) {
- $pos++;
- }
- }
- }
- # duplicate numbers - (1,4,5,6,7,33,45) becomes (1,1,4,5,6,7,33,33,45,45)
- for my $name (keys %cigars) {
- splice @{$cigars{$name}}, 1, 0, ${$cigars{$name}}[0] if
- ( ${$cigars{$name}}[0] + 1 < ${$cigars{$name}}[1] );
- push @{$cigars{$name}}, ${$cigars{$name}}[$#{$cigars{$name}}] if
- ( ${$cigars{$name}}[($#{$cigars{$name}} - 1)] + 1 <
- ${$cigars{$name}}[$#{$cigars{$name}}] );
- for ( my $x = 1 ; $x < $#{$cigars{$name}} - 1 ; $x++) {
- if (${$cigars{$name}}[$x - 1] + 1 < ${$cigars{$name}}[$x] &&
- ${$cigars{$name}}[$x + 1] > ${$cigars{$name}}[$x] + 1) {
- splice @{$cigars{$name}}, $x, 0, ${$cigars{$name}}[$x];
- }
- }
- }
- # collapse series - (1,1,4,5,6,7,33,33,45,45) becomes (1,1,4,7,33,33,45,45)
- for my $name (keys %cigars) {
- my @remove;
- for ( my $x = 0 ; $x < $#{$cigars{$name}} ; $x++) {
- if ( ${$cigars{$name}}[$x] == ${$cigars{$name}}[($x - 1)] + 1 &&
- ${$cigars{$name}}[$x] == ${$cigars{$name}}[($x + 1)] - 1 ) {
- unshift @remove,$x;
- }
- }
- for my $pos (@remove) {
- splice @{$cigars{$name}}, $pos, 1;
- }
- }
- # join and punctuate
- for my $name (keys %cigars) {
- my ($start,$end,$str) = "";
- while ( ($start,$end) = splice @{$cigars{$name}}, 0, 2 ) {
- $str .= ($start . "," . $end . ":");
- }
- $str =~ s/:$//;
- $cigars{$name} = $str;
- }
- %cigars;
- }
- =head2 match_line
- Title : match_line()
- Usage : $line = $align->match_line()
- Function : Generates a match line - much like consensus string
- except that a line indicating the '*' for a match.
- Args : (optional) Match line characters ('*' by default)
- (optional) Strong match char (':' by default)
- (optional) Weak match char ('.' by default)
- Returns : String
- =cut
- sub match_line {
- my ($self,$matchlinechar, $strong, $weak) = @_;
- my %matchchars = ('match' => $matchlinechar || '*',
- 'weak' => $weak || '.',
- 'strong' => $strong || ':',
- 'mismatch' => ' ',
- );
- my @seqchars;
- my $alphabet;
- foreach my $seq ( $self->each_seq ) {
- push @seqchars, [ split(//, uc ($seq->seq)) ];
- $alphabet = $seq->alphabet unless defined $alphabet;
- }
- my $refseq = shift @seqchars;
- # let's just march down the columns
- my $matchline;
- POS:
- foreach my $pos ( 0..$self->length ) {
- my $refchar = $refseq->[$pos];
- my $char = $matchchars{'mismatch'};
- unless( defined $refchar ) {
- last if $pos == $self->length; # short circuit on last residue
- # this in place to handle jason's soon-to-be-committed
- # intron mapping code
- goto bottom;
- }
- my %col = ($refchar => 1);
- my $dash = ($refchar eq '-' || $refchar eq '.' || $refchar eq ' ');
- foreach my $seq ( @seqchars ) {
- next if $pos >= scalar @$seq;
- $dash = 1 if( $seq->[$pos] eq '-' || $seq->[$pos] eq '.' ||
- $seq->[$pos] eq ' ' );
- $col{$seq->[$pos]}++ if defined $seq->[$pos];
- }
- my @colresidues = sort keys %col;
- # if all the values are the same
- if( $dash ) { $char = $matchchars{'mismatch'} }
- elsif( @colresidues == 1 ) { $char = $matchchars{'match'} }
- elsif( $alphabet eq 'protein' ) { # only try to do weak/strong
- # matches for protein seqs
- TYPE:
- foreach my $type ( qw(strong weak) ) {
- # iterate through categories
- my %groups;
- # iterate through each of the aa in the col
- # look to see which groups it is in
- foreach my $c ( @colresidues ) {
- foreach my $f ( grep { index($_,$c) >= 0 } @{$CONSERVATION_GROUPS{$type}} ) {
- push @{$groups{$f}},$c;
- }
- }
- GRP:
- foreach my $cols ( values %groups ) {
- @$cols = sort @$cols;
- # now we are just testing to see if two arrays
- # are identical w/o changing either one
- # have to be same len
- next if( scalar @$cols != scalar @colresidues );
- # walk down the length and check each slot
- for($_=0;$_ < (scalar @$cols);$_++ ) {
- next GRP if( $cols->[$_] ne $colresidues[$_] );
- }
- $char = $matchchars{$type};
- last TYPE;
- }
- }
- }
- bottom:
- $matchline .= $char;
- }
- return $matchline;
- }
- =head2 gap_line
- Title : gap_line()
- Usage : $line = $align->gap_line()
- Function : Generates a gap line - much like consensus string
- except that a line where '-' represents gap
- Args : (optional) gap line characters ('-' by default)
- Returns : string
- =cut
- sub gap_line {
- my ($self,$gapchar) = @_;
- $gapchar = $gapchar || $self->gap_char;
- my %gap_hsh; # column gaps vector
- foreach my $seq ( $self->each_seq ) {
- my $i = 0;
- map {$gap_hsh{$_->[0]} = undef} grep {$_->[1] eq $gapchar}
- map {[$i++, $_]} split(//, uc ($seq->seq));
- }
- my $gap_line;
- foreach my $pos ( 0..$self->length-1 ) {
- $gap_line .= (exists $gap_hsh{$pos}) ? $gapchar:'.';
- }
- return $gap_line;
- }
- =head2 all_gap_line
- Title : all_gap_line()
- Usage : $line = $align->all_gap_line()
- Function : Generates a gap line - much like consensus string
- except that a line where '-' represents all-gap column
- Args : (optional) gap line characters ('-' by default)
- Returns : string
- =cut
- sub all_gap_line {
- my ($self,$gapchar) = @_;
- $gapchar = $gapchar || $self->gap_char;
- my %gap_hsh; # column gaps counter hash
- my @seqs = $self->each_seq;
- foreach my $seq ( @seqs ) {
- my $i = 0;
- map {$gap_hsh{$_->[0]}++} grep {$_->[1] eq $gapchar}
- map {[$i++, $_]} split(//, uc ($seq->seq));
- }
- my …
Large files files are truncated, but you can click here to view the full file