diff variant_effect_predictor/Bio/Variation/SeqDiff.pm @ 0:1f6dce3d34e0

Uploaded
author mahtabm
date Thu, 11 Apr 2013 02:01:53 -0400
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/variant_effect_predictor/Bio/Variation/SeqDiff.pm	Thu Apr 11 02:01:53 2013 -0400
@@ -0,0 +1,1146 @@
+# $Id: SeqDiff.pm,v 1.16 2002/10/22 07:38:49 lapp Exp $
+# bioperl module for Bio::Variation::SeqDiff
+#
+# Cared for by Heikki Lehvaslaiho <heikki@ebi.ac.uk>
+#
+# Copyright Heikki Lehvaslaiho
+#
+# You may distribute this module under the same terms as perl itself
+#
+# POD documentation - main docs before the code
+
+# cds_end definition?
+
+=head1 NAME
+
+Bio::Variation::SeqDiff - Container class for mutation/variant descriptions
+
+=head1 SYNOPSIS
+
+  $seqDiff = Bio::Variation::SeqDiff->new (
+                                           -id => $M20132,
+					   -alphabet => 'rna',
+                                           -gene_symbol => 'AR'
+                                           -chromosome => 'X',
+                                           -numbering => 'coding'
+                                           );
+  # get a DNAMutation object somehow
+  $seqDiff->add_Variant($dnamut);
+  print  $seqDiff->sys_name(), "\n"; 
+
+=head1 DESCRIPTION
+
+SeqDiff stores Bio::Variation::VariantI object references and
+descriptive information common to all changes in a sequence. Mutations
+are understood to be any kind of sequence markers and are expected to
+occur in the same chromosome. See L<Bio::Variation::VariantI> for details.
+
+The methods of SeqDiff are geared towards describing mutations in
+human genes using gene-based coordinate system where 'A' of the
+initiator codon has number 1 and the one before it -1. This is
+according to conventions of human genetics.
+
+There will be class Bio::Variation::Genotype to describe markers in
+different chromosomes and diploid genototypes.
+
+Classes implementing Bio::Variation::VariantI interface are 
+Bio::Variation::DNAMutation, Bio::Variation::RNAChange, and
+Bio::Variation::AAChange. See L<Bio::Variation::VariantI>,
+L<Bio::Variation::DNAMutation>, L<Bio::Variation::RNAChange>, and
+L<Bio::Variation::AAChange> for more information.
+
+Variant objects can be added using two ways: an array passed to the
+constructor or as individual Variant objects with add_Variant
+method.
+
+
+=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 the 
+Bioperl mailing lists  Your participation is much appreciated.
+
+  bioperl-l@bioperl.org                         - General discussion
+  http://bio.perl.org/MailList.html             - About the mailing lists
+
+=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
+ email or the web:
+
+  bioperl-bugs@bio.perl.org
+  http://bugzilla.bioperl.org/
+
+=head1 AUTHOR - Heikki Lehvaslaiho
+
+Email:  heikki@ebi.ac.uk
+Address: 
+
+     EMBL Outstation, European Bioinformatics Institute
+     Wellcome Trust Genome Campus, Hinxton
+     Cambs. CB10 1SD, United Kingdom 
+
+=head1 CONTRIBUTORS
+
+Eckhard Lehmann, ecky@e-lehmann.de
+
+=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::Variation::SeqDiff;
+my $VERSION=1.0;
+
+use strict;
+use vars qw($VERSION @ISA);
+use Bio::Root::Root;
+use Bio::Tools::CodonTable;
+use Bio::PrimarySeq;
+
+@ISA = qw( Bio::Root::Root );
+
+
+=head2 new
+
+  Title   : new
+  Usage   : $seqDiff = Bio::Variation::SeqDiff->new;
+  Function: generates a new Bio::Variation::SeqDiff
+  Returns : reference to a new object of class SeqDiff
+  Args    : 
+
+=cut
+
+sub new {
+    my($class,@args) = @_;
+    my $self = $class->SUPER::new(@args);
+
+    my($id, $sysname, $trivname, $chr, $gene_symbol, 
+       $desc, $alphabet, $numbering, $offset, $rna_offset, $rna_id, $cds_end,
+       $dna_ori, $dna_mut, $rna_ori, $rna_mut, $aa_ori, $aa_mut
+       #@variants, @genes
+       ) =
+	   $self->_rearrange([qw(ID
+				 SYSNAME
+				 TRIVNAME
+				 CHR
+				 GENE_SYMBOL
+				 DESC
+				 ALPHABET
+				 NUMBERING
+				 OFFSET
+				 RNA_OFFSET
+				 RNA_ID
+				 CDS_END
+				 DNA_ORI
+				 DNA_MUT
+				 RNA_ORI
+				 AA_ORI
+				 AA_MUT
+				 )],
+			    @args);
+    
+    #my $make = $self->SUPER::_initialize(@args);
+    
+    $id        && $self->id($id);           
+    $sysname   && $self->sysname($sysname); 
+    $trivname  && $self->trivname($trivname);
+    $chr       && $self->chromosome($chr);  
+    $gene_symbol && $self->gene_symbol($chr);
+    $desc      && $self->description($desc);
+    $alphabet   && $self->alphabet($alphabet);
+    $numbering && $self->numbering($numbering);
+    $offset    && $self->offset($offset);   
+    $rna_offset && $self->rna_offset($rna_offset);   
+    $rna_id    && $self->rna_id($rna_id);   
+    $cds_end   && $self->cds_end($cds_end);   
+
+    $dna_ori   && $self->dna_ori($dna_ori); 
+    $dna_mut   && $self->dna_mut($dna_mut); 
+    $rna_ori   && $self->rna_ori($rna_ori); 
+    $rna_mut   && $self->rna_mut($rna_mut); 
+    $aa_ori    && $self->aa_ori ($aa_ori);  
+    $aa_mut    && $self->aa_mut ($aa_mut);  
+
+    $self->{ 'variants' } = [];
+    #@variants && push(@{$self->{'variants'}},@variants);
+
+    $self->{ 'genes' } = [];
+    #@genes && push(@{$self->{'genes'}},@genes);
+
+    return $self; # success - we hope!
+}
+
+
+=head2 id
+
+ Title   : id
+ Usage   : $obj->id(H0001); $id = $obj->id();
+ Function: 
+
+           Sets or returns the id of the seqDiff.
+           Should be used to give the collection of variants a UID
+           without semantic associations.
+
+ Example : 
+ Returns : value of id, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub id {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'id'} = $value;
+  }
+#  unless (exists $self->{'id'}) {
+#    return "undefined";
+#  }
+ else {
+    return $self->{'id'};
+  }
+}
+
+
+=head2 sysname
+
+ Title   : sysname
+ Usage   : $obj->sysname('5C>G'); $sysname = $obj->sysname();
+ Function: 
+
+           Sets or returns the systematic name of the seqDiff.  The
+           name should follow the HUGO Mutation Database Initiative
+           approved nomenclature. If called without first setting the
+           value, will generate it from L<Bio::Variation::DNAMutation>
+           objects attached.
+
+ Example : 
+ Returns : value of sysname, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub sysname {
+    my ($self,$value) = @_;
+    if (defined $value) {
+	$self->{'sysname'} = $value;
+    }
+    elsif (not defined $self->{'sysname'}) {
+
+	my $sysname = ''; 
+	my $c = 0;
+	foreach my $mut ($self->each_Variant) {
+	    if( $mut->isa('Bio::Variation::DNAMutation') ) {
+		$c++;
+		if ($c == 1 ) {
+		    $sysname = $mut->sysname ;
+		}
+		else {
+		    $sysname .= ";". $mut->sysname;
+		}
+	    }
+	}
+	$sysname  = "[". $sysname. "]" if $c > 1;
+	$self->{'sysname'} = $sysname;
+    }
+    return $self->{'sysname'};
+}
+
+
+=head2 trivname
+
+ Title   : trivname
+ Usage   : $obj->trivname('[A2G;T56G]'); $trivname = $obj->trivname();
+ Function: 
+
+           Sets or returns the trivial name of the seqDiff.
+           The name should follow the HUGO Mutation Database Initiative
+           approved nomenclature. If called without first setting the
+           value, will generate it from L<Bio::Variation::AAChange>
+           objects attached.
+
+ Example : 
+ Returns : value of trivname, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub trivname {
+    my ($self,$value) = @_;
+    if (defined $value) {
+	$self->{'trivname'} = $value;
+    }
+    elsif (not defined $self->{'trivname'}) {
+	
+	my $trivname = ''; 
+	my $c = 0;
+	foreach my $mut ($self->each_Variant) {
+	    if( $mut->isa('Bio::Variation::AAChange') ) {
+		$c++;
+		if ($c == 1 ) {
+		    $trivname = $mut->trivname ;
+		}
+		else {
+		    $trivname .= ";". $mut->trivname;
+		}
+	    }
+	}
+	$trivname  = "[". $trivname. "]" if $c > 1;
+	$self->{'trivname'} = $trivname;
+    }
+
+  else {
+      return $self->{'trivname'};
+  }
+}
+
+
+=head2 chromosome
+
+ Title   : chromosome
+ Usage   : $obj->chromosome('X'); $chromosome = $obj->chromosome();
+ Function: 
+
+           Sets or returns the chromosome ("linkage group") of the seqDiff.
+
+ Example : 
+ Returns : value of chromosome, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub chromosome {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'chromosome'} = $value;
+  }
+  else {
+      return $self->{'chromosome'};
+  }
+}
+
+
+=head2 gene_symbol
+
+ Title   : gene_symbol
+ Usage   : $obj->gene_symbol('FOS'); $gene_symbol = $obj->gene_symbol;
+ Function: 
+
+           Sets or returns the gene symbol for the studied CDS.
+
+ Example : 
+ Returns : value of gene_symbol, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub gene_symbol {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'gene_symbol'} = $value;
+  }
+  else {
+      return $self->{'gene_symbol'};
+  }
+}
+
+
+
+=head2 description
+
+ Title   : description
+ Usage   : $obj->description('short description'); $descr = $obj->description();
+ Function: 
+
+           Sets or returns the short description of the seqDiff.
+
+ Example : 
+ Returns : value of description, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub description {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'description'} = $value;
+  }
+  else {
+      return $self->{'description'};
+  }
+}
+
+
+=head2 alphabet
+
+ Title   : alphabet
+ Usage   : if( $obj->alphabet eq 'dna' ) { /Do Something/ }
+ Function: Returns the type of primary reference sequence being one of 
+           'dna', 'rna' or 'protein'. This is case sensitive.
+
+ Returns : a string either 'dna','rna','protein'. 
+ Args    : none
+
+
+=cut
+
+sub alphabet {
+   my ($self,$value) = @_;
+   my %type = (dna => 1,
+	       rna => 1,
+	       protein => 1);
+   if( defined $value ) {
+       if ($type{$value}) {
+	   $self->{'alphabet'} = $value;
+       } else {
+	   $self->throw("$value is not valid alphabet value!");
+       }
+   }
+   return $self->{'alphabet'};
+}
+
+
+=head2 numbering
+
+ Title   : numbering
+ Usage   : $obj->numbering('coding'); $numbering = $obj->numbering();
+ Function: 
+
+           Sets or returns the string giving the numbering schema used
+           to describe the variants.
+
+ Example : 
+ Returns : value of numbering, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+
+sub numbering {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'numbering'} = $value;
+  }
+  else {
+      return $self->{'numbering'};
+  }
+}
+
+
+=head2 offset
+
+ Title   : offset
+ Usage   : $obj->offset(124); $offset = $obj->offset();
+ Function: 
+
+           Sets or returns the offset from the beginning of the DNA sequence 
+           to the coordinate start used to describe variants. Typically
+           the beginning of the coding region of the gene. 
+           The cds_start should be 1 + offset.
+
+ Example : 
+ Returns : value of offset, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+
+sub offset {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'offset'} = $value;
+  }
+  elsif (not defined $self->{'offset'} ) {
+      return $self->{'offset'} = 0;
+  }
+  else {
+      return $self->{'offset'};
+  }
+}
+
+
+=head2 cds_start
+
+ Title   : cds_start
+ Usage   : $obj->cds_start(123); $cds_start = $obj->cds_start();
+ Function: 
+
+           Sets or returns the cds_start from the beginning of the DNA
+           sequence to the coordinate start used to describe
+           variants. Typically the beginning of the coding region of
+           the gene. Needs to be and is implemented as 1 + offset.
+
+ Example : 
+ Returns : value of cds_start, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+
+sub cds_start {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'offset'} = $value - 1;
+  }
+  else {
+      return $self->{'offset'} + 1;
+  }
+}
+
+
+=head2 cds_end
+
+ Title   : cds_end
+ Usage   : $obj->cds_end(321); $cds_end = $obj->cds_end();
+ Function: 
+
+           Sets or returns the position of the last nucleotitide of the
+           termination codon. The coordinate system starts from cds_start.
+
+ Example : 
+ Returns : value of cds_end, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+
+sub cds_end {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'cds_end'} = $value;
+  }
+  else {
+      return $self->{'cds_end'};
+      #$self->{'cds_end'} = CORE::length($self->SeqDiff->rna_ori)/3;
+  }
+}
+
+
+=head2 rna_offset
+
+ Title   : rna_offset
+ Usage   : $obj->rna_offset(124); $rna_offset = $obj->rna_offset();
+ Function: 
+
+           Sets or returns the rna_offset from the beginning of the RNA sequence 
+           to the coordinate start used to describe variants. Typically
+           the beginning of the coding region of the gene. 
+
+ Example : 
+ Returns : value of rna_offset, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+
+sub rna_offset {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'rna_offset'} = $value;
+  }
+  elsif (not defined $self->{'rna_offset'} ) {
+      return $self->{'rna_offset'} = 0;
+  }
+  else {
+      return $self->{'rna_offset'};
+  }
+}
+
+
+=head2 rna_id
+
+ Title   : rna_id
+ Usage   : $obj->rna_id('transcript#3'); $rna_id = $obj->rna_id();
+ Function: 
+
+	    Sets or returns the ID for original RNA sequence of the seqDiff.
+
+ Example : 
+ Returns : value of rna_id, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub rna_id {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'rna_id'} = $value;
+  }
+  else {
+      return $self->{'rna_id'};
+  }
+}
+
+
+
+=head2 add_Variant
+
+ Title   : add_Variant
+ Usage   : $obj->add_Variant($variant)
+ Function: 
+
+           Pushes one Bio::Variation::Variant into the list of variants.
+           At the same time, creates a link from the Variant to SeqDiff
+           using its SeqDiff method.
+
+ Example : 
+ Returns : 1 when succeeds, 0 for failure.
+ Args    : Variant object
+
+=cut
+
+sub add_Variant {
+  my ($self,$value) = @_;
+  if (defined $value) {
+      if( ! $value->isa('Bio::Variation::VariantI') ) {
+	  $self->throw("Is not a VariantI complying  object but a [$self]");
+	  return 0;
+      }
+      else {
+	  push(@{$self->{'variants'}},$value);
+	  $value->SeqDiff($self);
+	  return 1;
+      }
+  }
+  else {
+      return 0;
+  }
+}
+
+
+=head2 each_Variant
+
+ Title   : each_Variant
+ Usage   : $obj->each_Variant();
+ Function: 
+
+            Returns a list of Variants.
+
+ Example : 
+ Returns : list of Variants
+ Args    : none
+
+=cut
+
+sub each_Variant{
+   my ($self,@args) = @_;
+   
+   return @{$self->{'variants'}}; 
+}
+
+
+
+=head2 add_Gene
+
+ Title   : add_Gene
+ Usage   : $obj->add_Gene($gene)
+ Function: 
+
+           Pushes one L<Bio::LiveSeq::Gene> into the list of genes.
+
+ Example : 
+ Returns : 1 when succeeds, 0 for failure.
+ Args    : Bio::LiveSeq::Gene object
+
+See L<Bio::LiveSeq::Gene> for more information.
+
+=cut
+
+
+sub add_Gene {
+  my ($self,$value) = @_;
+  if (defined $value) {
+      if( ! $value->isa('Bio::LiveSeq::Gene') ) {
+	  $value->throw("Is not a Bio::LiveSeq::Gene object but a  [$value]");
+	  return 0;
+      }
+      else {
+	  push(@{$self->{'genes'}},$value);
+	  return 1;
+      }
+  }
+  else {
+      return 0;
+  }
+}
+
+
+=head2 each_Gene
+
+ Title   : each_Gene
+ Usage   : $obj->each_Gene();
+ Function: 
+
+            Returns a list of L<Bio::LiveSeq::Gene>s.
+
+ Example : 
+ Returns : list of Genes
+ Args    : none
+
+=cut
+
+sub each_Gene{
+   my ($self,@args) = @_;
+
+   return @{$self->{'genes'}}; 
+}
+
+
+=head2 dna_ori
+
+ Title   : dna_ori
+ Usage   : $obj->dna_ori('atgctgctgctgct'); $dna_ori = $obj->dna_ori();
+ Function: 
+
+	    Sets or returns the original DNA sequence string of the seqDiff.
+
+ Example : 
+ Returns : value of dna_ori, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub dna_ori {
+  my ($self,$value) = @_;
+  if (defined $value) {
+      $self->{'dna_ori'} = $value;
+  }
+  else {
+      return $self->{'dna_ori'};
+  }
+}
+
+
+=head2 dna_mut
+
+ Title   : dna_mut
+ Usage   : $obj->dna_mut('atgctggtgctgct'); $dna_mut = $obj->dna_mut();
+ Function: 
+
+	    Sets or returns the mutated DNA sequence of the seqDiff.
+            If sequence has not been set generates it from the
+            original sequence and DNA mutations.
+
+ Example : 
+ Returns : value of dna_mut, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub dna_mut {
+  my ($self,$value) = @_;
+  if (defined $value) {
+      $self->{'dna_mut'} = $value;
+  }
+  else {
+      $self->_set_dnamut() unless $self->{'dna_mut'};
+      return $self->{'dna_mut'};
+  }
+}
+
+sub _set_dnamut {
+    my $self = shift;
+
+    return undef unless $self->{'dna_ori'}  && $self->each_Variant;
+
+    $self->{'dna_mut'} = $self->{'dna_ori'};
+    foreach ($self->each_Variant) {
+	next unless $_->isa('Bio::Variation::DNAMutation');
+	next unless $_->isMutation;
+
+	my ($s, $la, $le);
+	#lies the mutation less than 25 bases after the start of sequence?
+	if ($_->start < 25) {
+	    $s = 0; $la = $_->start - 1;
+	} else {
+	    $s = $_->start - 25; $la = 25;
+	}
+
+	#is the mutation an insertion?
+	$_->end($_->start) unless $_->allele_ori->seq;
+
+	#does the mutation end greater than 25 bases before the end of
+	#sequence?
+	if (($_->end + 25) > length($self->{'dna_mut'})) {
+	    $le = length($self->{'dna_mut'}) - $_->end;
+	} else {
+	    $le = 25;
+	}
+
+	$_->dnStreamSeq(substr($self->{'dna_mut'}, $s, $la));
+	$_->upStreamSeq(substr($self->{'dna_mut'}, $_->end, $le));
+
+	my $s_ori = $_->dnStreamSeq . $_->allele_ori->seq . $_->upStreamSeq;
+	my $s_mut = $_->dnStreamSeq . $_->allele_mut->seq . $_->upStreamSeq;
+
+	(my $str = $self->{'dna_mut'}) =~ s/$s_ori/$s_mut/;
+	$self->{'dna_mut'} = $str;
+    }
+}
+
+
+=head2 rna_ori
+
+ Title   : rna_ori
+ Usage   : $obj->rna_ori('atgctgctgctgct'); $rna_ori = $obj->rna_ori();
+ Function: 
+
+	    Sets or returns the original RNA sequence of the seqDiff.
+
+ Example : 
+ Returns : value of rna_ori, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub rna_ori {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'rna_ori'} = $value;
+  }
+  else {
+      return $self->{'rna_ori'};
+  }
+}
+
+
+=head2 rna_mut
+
+ Title   : rna_mut
+ Usage   : $obj->rna_mut('atgctggtgctgct'); $rna_mut = $obj->rna_mut();
+ Function: 
+
+	    Sets or returns the mutated RNA sequence of the seqDiff.
+
+ Example : 
+ Returns : value of rna_mut, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub rna_mut {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'rna_mut'} = $value;
+  }
+  else {
+      return $self->{'rna_mut'};
+  }
+}
+
+
+=head2 aa_ori
+
+ Title   : aa_ori
+ Usage   : $obj->aa_ori('MAGVLL*'); $aa_ori = $obj->aa_ori();
+ Function: 
+
+	    Sets or returns the original protein sequence of the seqDiff.
+
+ Example : 
+ Returns : value of aa_ori, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub aa_ori {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'aa_ori'} = $value;
+  }
+  else {
+      return $self->{'aa_ori'};
+  }
+}
+
+
+=head2 aa_mut
+
+ Title   : aa_mut
+ Usage   : $obj->aa_mut('MA*'); $aa_mut = $obj->aa_mut();
+ Function: 
+
+	    Sets or returns the mutated protein sequence of the seqDiff.
+
+ Example : 
+ Returns : value of aa_mut, a scalar
+ Args    : newvalue (optional)
+
+=cut
+
+
+sub aa_mut {
+  my ($self,$value) = @_;
+  if (defined $value) {
+    $self->{'aa_mut'} = $value;
+  }
+  else {
+      return $self->{'aa_mut'};
+  }
+}
+
+
+=head2 seqobj
+
+ Title   : seqobj
+ Usage   : $dnaobj = $obj->seqobj('dna_mut');
+ Function: 
+
+	    Returns the any original or mutated sequences as a
+	    Bio::PrimarySeq object.
+
+ Example : 
+ Returns : Bio::PrimarySeq object for the requested sequence
+ Args    : string, method name for the sequence requested
+
+See L<Bio::PrimarySeq> for more information.
+
+=cut
+
+sub seqobj {
+  my ($self,$value) = @_;
+  my $out;
+  my %valid_obj = 
+      map {$_, 1} qw(dna_ori rna_ori aa_ori dna_mut rna_mut aa_mut);
+  $valid_obj{$value} ||
+      $self->throw("Sequence type '$value' is not a valid type (".
+                  join(',', map "'$_'", sort keys %valid_obj) .") lowercase");
+  my ($alphabet) = $value =~ /([^_]+)/;
+  my $id =  $self->id;
+  $id =  $self->rna_id if $self->rna_id;
+  $alphabet = 'protein' if $alphabet eq 'aa';
+  $out = Bio::PrimarySeq->new
+      ( '-seq' => $self->{$value},
+	'-display_id'  => $id,
+	'-accession_number' => $self->id,
+	'-alphabet' => $alphabet
+	) if   $self->{$value} ;
+  return $out;
+}
+
+=head2 alignment
+
+ Title   : alignment
+ Usage   : $obj->alignment
+ Function: 
+
+           Returns a pretty RNA/AA sequence alignment from linked
+           objects.  Under construction: Only simple coding region
+           point mutations work.
+
+ Example : 
+ Returns : 
+ Args    : none
+
+=cut
+
+
+sub alignment {
+    my $self = shift;
+    my (@entry, $text);
+
+    my $maxflanklen = 12;
+
+    foreach my $mut ($self->each_Variant) {
+	if( $mut->isa('Bio::Variation::RNAChange') ) {
+
+	    my $upflank = $mut->upStreamSeq;
+	    my $dnflank = $mut->dnStreamSeq;
+	    my $cposd = $mut->codon_pos;
+	    my $rori = $mut->allele_ori->seq;
+	    my $rmut =  $mut->allele_mut->seq;
+	    my $rseqoriu = '';
+	    my $rseqmutu = '';
+	    my $rseqorid = '';
+	    my $rseqmutd = '';
+	    my $aaseqmutu = '';
+	    my (@rseqori, @rseqmut );
+
+	    #  point
+	    if ($mut->DNAMutation->label =~ /point/) {
+		if ($cposd == 1 ) {
+		    my $nt2d = substr($dnflank, 0, 2);
+		    push @rseqori, $rori. $nt2d;
+		    push @rseqmut, uc ($rmut). $nt2d;
+		    $dnflank = substr($dnflank, 2);
+		}
+		elsif ($cposd == 2) { 
+		    my $ntu = chop $upflank;
+		    my $ntd = substr($dnflank, 0, 1);
+		    push @rseqori, $ntu. $rori. $ntd;
+		    push @rseqmut,  $ntu. uc ($rmut). $ntd;
+		    $dnflank =  substr($dnflank, 1);
+		}
+		elsif ($cposd == 3) {
+		    my $ntu1 = chop $upflank;
+		    my $ntu2 = chop $upflank;
+		    push (@rseqori, $ntu2. $ntu1. $rori);
+		    push (@rseqmut, $ntu2. $ntu1. uc $rmut);
+		}		
+	    }
+	    #deletion
+	    elsif ($mut->DNAMutation->label =~ /deletion/) {
+		if ($cposd == 2 ) {
+		    $rseqorid = chop $upflank;
+		    $rseqmutd = $rseqorid;
+		}
+		for (my $i=1; $i<=$mut->length; $i++) {
+		    my $ntd .= substr($mut->allele_ori, $i-1, 1);
+		    $rseqorid .= $ntd;
+		    if  (length($rseqorid) == 3 ) {
+			push (@rseqori, $rseqorid);
+			push (@rseqmut, "   ");
+			$rseqorid = '';
+		    }		    
+		}
+
+		if ($rseqorid) {
+		    $rseqorid .= substr($dnflank, 0, 3-$rseqorid);
+		    push (@rseqori, $rseqorid);
+		    push (@rseqmut, "   ");
+		    $dnflank = substr($dnflank,3-$rseqorid);
+		} 
+	    }
+	    $upflank = reverse $upflank;
+	    # loop throught the flanks
+	    for (my $i=1; $i<=length($dnflank); $i++) {
+		
+		last if  $i > $maxflanklen;
+
+		my $ntd .= substr($dnflank, $i-1, 1);
+		my $ntu .= substr($upflank, $i-1, 1);
+
+		$rseqmutd .= $ntd;
+		$rseqorid .= $ntd;
+		$rseqmutu = $ntu. $rseqmutu;
+		$rseqoriu = $ntu. $rseqoriu;
+		
+		if  (length($rseqorid) == 3  and length($rseqorid) == 3) {
+		    push (@rseqori, $rseqorid);
+		    push (@rseqmut, $rseqmutd);
+		    $rseqorid =  $rseqmutd ='';
+		}
+		if  (length($rseqoriu) == 3  and length($rseqoriu) == 3) {
+		    unshift (@rseqori, $rseqoriu);
+		    unshift (@rseqmut, $rseqmutu);
+		    $rseqoriu =  $rseqmutu ='';
+		}
+
+		#print "|i=$i,  $cposd, $rseqmutd, $rseqorid\n";
+		#print "|i=$i,  $cposu, $rseqmutu, $rseqoriu\n\n";
+
+	    }
+
+	    push (@rseqori, $rseqorid);
+	    unshift (@rseqori, $rseqoriu);
+	    push (@rseqmut, $rseqmutd);
+	    unshift (@rseqmut, $rseqmutu);
+	    
+	    return unless $mut->AAChange;
+	    #translate
+	    my $tr = new Bio::Tools::CodonTable ('-id' => $mut->codon_table);
+	    my $apos =  $mut->AAChange->start;
+	    my $aposmax = CORE::length($self->aa_ori); #terminator codon no 
+	    my $rseqori;
+	    my $rseqmut;
+	    my $aaseqori;
+	    my $aaseqmut = "";
+	    for (my $i = 0; $i <= $#rseqori; $i++) {
+		 my $a = '';
+
+		 $a =  $tr->translate($rseqori[$i]) if length($rseqori[$i]) == 3;
+		 
+		 if (length($a) != 1 or 
+		     $apos - ( $maxflanklen/2 -1) + $i < 1 or 
+		     $apos - ( $maxflanklen/2 -1) + $i > $aposmax ) {
+		     $aaseqori .= "    ";
+		 } else {
+		     $aaseqori .= " ". $a. "  ";
+		 }
+		 my $b = '';
+		 if (length($rseqmut[$i]) == 3) {
+		     if ($rseqmut[$i] eq '   ') {
+			 $b = "_";
+		     } else {
+			 $b = $tr->translate($rseqmut[$i]);
+		     }
+		 }
+		 if (( $b ne $a and
+		       length($b) == 1 and 
+		       $apos - ( $maxflanklen/2 -1) + $i >= 1 ) or
+		     ( $apos - ( $maxflanklen/2 -1) + $i >= $aposmax and 
+		       $mut->label =~ 'termination')
+		     ) {
+		     $aaseqmut .= " ". $b. "  ";
+		 } else {
+		     $aaseqmut .= "    ";
+		 }
+		 
+		 if ($i == 0 and length($rseqori[$i]) != 3) {
+		     my $l = 3 - length($rseqori[$i]);
+		     $rseqori[$i] = (" " x $l). $rseqori[$i];
+		     $rseqmut[$i] = (" " x $l). $rseqmut[$i];
+		 }
+		 $rseqori .= $rseqori[$i]. " " if $rseqori[$i] ne '';
+		 $rseqmut .= $rseqmut[$i]. " " if $rseqmut[$i] ne '';
+	     }
+	    
+	    # collect the results
+	    push (@entry, 
+		  "\n"
+		  );   	    
+	    $text = "           ". $aaseqmut; 
+	    push (@entry, 
+		  $text
+		  );   	    
+	    $text = "Variant  : ". $rseqmut;
+	    push (@entry, 
+		  $text
+		  );   	    
+	    $text = "Reference: ". $rseqori;
+	    push (@entry, 
+		  $text
+		  );   	    
+	    $text = "           ". $aaseqori;
+	    push (@entry, 
+		  $text
+		  );   
+	    push (@entry, 
+		  "\n"
+		  );   
+	}
+
+    }
+
+    my $res;
+    foreach my $line (@entry) {
+       $res .=  "$line\n";
+    }
+    return $res;
+}
+
+1;