# $gene is a Bio::LiveSeq::Gene object
  my $mutate = Bio::LiveSeq::Mutator->new('-gene' => $gene,
  					  '-numbering' => "coding"
  					   );
  # $mut is a Bio::LiveSeq::Mutation object
  $mutate->add_Mutation($mut);
  # $results is a Bio::Variation::SeqDiff object
  my $results=$mutate->change_gene();
  if ($results) {
      my $out = Bio::Variation::IO->new( '-format' => 'flat');
      $out->write($results);
  }

  my ($self,$value) = @_;
  if (defined $value) {
      if( ! $value->isa('Bio::LiveSeq::DNA') and ! $value->isa('Bio::LiveSeq::Transcript') ) {
	  $self->throw("Is not a Bio::LiveSeq::DNA/Transcript object but a [$value]");
	  return undef;
      }
      else {
	  $self->{'DNA'} = $value;
      }
  }
  unless (exists $self->{'DNA'}) {
      return (undef);
  } else {
      return $self->{'DNA'};
  }

  my ($self,$value) = @_;
  if (defined $value) {
      if( ! $value->isa('Bio::LiveSeq::Transcript') ) {
	  $self->throw("Is not a Bio::LiveSeq::RNA/Transcript object but a [$value]");
	  return undef;
      }
      else {
	  $self->{'RNA'} = $value;
      }
  }
  unless (exists $self->{'RNA'}) {
      return (undef);
  } else {
      return $self->{'RNA'};
  }

    my ($self, $refseq, $SeqDiff) = @_;
    my $count;
    my @bb = @{$self->{'mutations'}};
    my $cc = scalar @bb;
    foreach my $mut (@{$self->{'mutations'}}) {
#	 if ($self->numbering eq 'gene' and $mut->pos < 1) {
#	     my $tmp = $mut->pos;
#	     print STDERR "pos: ", "$tmp\n";
#	     $tmp++ if $tmp < 1;
#	     $tmp += $SeqDiff->offset;
#	     print STDERR "pos2: ", "$tmp\n";
#	     $mut->pos($tmp);
#	 }
#	elsif ($self->numbering eq 'entry') {
	if ($self->numbering eq 'entry') {
	    my $tmp = $mut->pos;
	    $tmp -= $SeqDiff->offset;
	    $tmp-- if $tmp < 1;
	    $mut->pos($tmp);
	}

	my $label = $refseq->label($mut->pos); # get the label for the position
	$mut->label($label), $count++ if $label > 0 ;
	#print STDERR "x", $mut->pos,'|' ,$mut->label, "\n";
    }
    return $count;
}

#
# Calculate labels around mutated nucleotide
#

    my ($self, $seqDiff) = @_;

    if ($self->rnachange and $self->rnachange->region eq 'coding') {

	#$seqDiff->add_Variant($self->rnachange);

	my $aachange=$self->aachange;
	my ($AAobj,$aa_start_prelabel,$aa_start,$mut_translation);
	$AAobj=$self->RNA->get_Translation;
	$aa_start_prelabel=$aachange->start;
	$aa_start=$AAobj->position($self->RNA->label(2,$aa_start_prelabel));
	$aachange->start($aa_start);
	$mut_translation=$AAobj->seq;

	# this now takes in account possible preinsertions
	my $aa_m = Bio::Variation::Allele->new;
	$aa_m->seq(substr($mut_translation,$aa_start-1)) if substr($mut_translation,$aa_start-1);
	$aachange->allele_mut($aa_m);
	$aachange->add_Allele($aa_m);
	#$aachange->allele_mut(substr($mut_translation,$aa_start-1));
	#$aachange->allele_mut($mut_translation);
	my ($rlenori, $rlenmut);
	$rlenori = CORE::length($aachange->RNAChange->allele_ori->seq);
	$rlenmut = CORE::length($aachange->RNAChange->allele_mut->seq);
	#point mutation

	if ($rlenori == 1 and $rlenmut == 1 and $aachange->allele_ori->seq ne '*') {
	     my $alleleseq;
	     if ($aachange->allele_mut->seq) {
		 $alleleseq = substr($aachange->allele_mut->seq, 0, 1);
		 $aachange->allele_mut->seq($alleleseq);
	     }
	     $aachange->end($aachange->start);
	     $aachange->length(1);
	 }
	elsif ( $rlenori == $rlenmut and 
		$aachange->allele_ori->seq ne '*' ) { #complex inframe mutation
	    $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 
				       0, 
				       length($aachange->allele_ori->seq));
	}
	#inframe mutation
	elsif ((int($rlenori-$rlenmut))%3 == 0) {
	    if ($aachange->RNAChange->allele_mut->seq  and
		$aachange->RNAChange->allele_ori->seq ) {
		# complex
		my $rna_len = length ($aachange->RNAChange->allele_mut->seq);
		my $len = $rna_len/3;
		$len++ unless $rna_len%3 == 0;
		$aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0, $len );
	    }
	    elsif  ($aachange->RNAChange->codon_pos == 1){
		 # deletion
		if ($aachange->RNAChange->allele_mut->seq eq  '') {
		    $aachange->allele_mut->seq('');
		    $aachange->end($aachange->start + $aachange->length - 1 );
		}
		 # insertion
		 elsif ($aachange->RNAChange->allele_ori->seq eq '' ) {
		     $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0,
					   length ($aachange->RNAChange->allele_mut->seq) / 3);
		     $aachange->allele_ori->seq('');
		     $aachange->end($aachange->start + $aachange->length - 1 );
		     $aachange->length(0);
		 }
	    } else {
		#elsif  ($aachange->RNAChange->codon_pos == 2){
		 # deletion
		 if (not $aachange->RNAChange->allele_mut->seq ) {
		     $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0, 1);
		 }
		 # insertion
		 elsif (not $aachange->RNAChange->allele_ori->seq) {
		     $aachange->allele_mut->seq(substr $aachange->allele_mut->seq, 0,
						length ($aachange->RNAChange->allele_mut->seq) / 3 +1);
		 }
	     }
	 } else {
	     #frameshift
	     #my $pos = index $aachange->allele_mut
	     #$aachange->allele_mut(substr($aachange->allele_mut, 0, 1));
	     $aachange->length(CORE::length($aachange->allele_ori->seq));
	     my $aaend = $aachange->start + $aachange->length -1;
	     $aachange->end($aachange->start);
	 }

	 # splicing site deletion check
	 my @beforeexons=@{$self->exons};
	 my @afterexons=$self->RNA->all_Exons;
	 my $i;
	 if (scalar(@beforeexons) ne scalar(@afterexons)) {
	     my $mut_number = $self->mutation->issue;
	     $self->warn("Exons have been modified at mutation n.$mut_number!");
	     $self->rnachange->exons_modified(1);
	 } else {
	   EXONCHECK:
	     foreach $i (0..$#beforeexons) {
		 if ($beforeexons[$i] ne $afterexons[$i]) {
	     my $mut_number = $self->mutation->issue;
		     $self->warn("Exons have been modified at mutation n.$mut_number!");
		     $self->rnachange->exons_modified(1);
		     last EXONCHECK;
		 }
	     }
	 }
     } else {
	 #$seqDiff->rnachange(undef);
	 #print "getting here?";
     }
    return 1;
}

1;

    my ($self) = @_;
    my @exons=$self->RNA->all_Exons;
    my $RNAstart=$self->RNA->start;
    my $RNAend=$self->RNA->end;
    my ($firstexon,$before,$after,$i);
    my ($rnaAffected) = 0;

    # check for inserted labels (that require follows instead of <,>)
    my $DNAend=$self->RNA->{'seq'}->end;
    if ($self->mutation->prelabel > $DNAend or $self->mutation->postlabel > $DNAend) {
	 #this means one of the two labels is an inserted one
	 #(coming from a previous mutation. This would falsify all <,>
	 #checks, so the follow() has to be used
	 $self->warn("Attention, workaround not fully tested yet! Expect unpredictable results.\n");
	 if (($self->mutation->postlabel==$RNAstart) or (follows($self->mutation->postlabel,$RNAstart))) {
	     $self->warn("RNA not affected because change occurs before RNAstart");
	 }
	 elsif (($RNAend==$self->mutation->prelabel) or (follows($RNAend,$self->mutation->prelabel))) {
	     $self->warn("RNA not affected because change occurs after RNAend");
	 }
	 elsif (scalar @exons == 1) {
	     #no introns, just one exon
	     $rnaAffected = 1; # then RNA is affected!
	 } else {
	     # otherwise check for change occurring inside an intron
	     $firstexon=shift(@exons);
	     $before=$firstexon->end;
	
	     foreach $i (0..$#exons) {
		 $after=$exons[$i]->start;
		 if (follows($self->mutation->prelabel,$before) or
			($after==$self->mutation->prelabel) or
			follows($after,$self->mutation->prelabel) or
			follows($after,$self->mutation->postlabel)) {

		     $rnaAffected = 1;
		     # $i is number of exon and can be used for proximity check
		 }
		 $before=$exons[$i]->end;
	     }
	
	 }
    } else {
	my $strand = $exons[0]->strand;
	if (($strand == 1 and $self->mutation->postlabel <= $RNAstart) or
	    ($strand != 1 and $self->mutation->postlabel >= $RNAstart)) {
	    #$self->warn("RNA not affected because change occurs before RNAstart");
	    $rnaAffected = 0;
	}
	elsif (($strand == 1 and $self->mutation->prelabel >= $RNAend) or
		($strand != 1 and $self->mutation->prelabel <= $RNAend)) {
	     #$self->warn("RNA not affected because change occurs after RNAend");
	     $rnaAffected = 0;
	     my $dist;
	     if ($strand == 1){
		 $dist = $self->mutation->prelabel - $RNAend;
	     } else {
		 $dist = $RNAend - $self->mutation->prelabel;
	     }
	     $self->dnamut->region_dist($dist);
	 }
	 elsif (scalar @exons == 1) {
	     #if just one exon -> no introns, 
	     $rnaAffected = 1; # then RNA is affected!
	 } else {	
	     # otherwise check for mutation occurring inside an intron
	     $firstexon=shift(@exons);
	     $before=$firstexon->end;
	     if ( ($strand == 1 and $self->mutation->prelabel < $before) or 
		  ($strand == -1 and $self->mutation->prelabel > $before) 
		  ) {
		 $rnaAffected = 1 ;

		 #print "Exon 1 : ", $firstexon->start, " - ", $firstexon->end, "<br>\n";
		 my $afterdist = $self->mutation->prelabel - $firstexon->start;
		 my $beforedist =  $firstexon->end - $self->mutation->postlabel;
		 my $exonvalue = $i + 1;
		 $self->dnamut->region('exon');
		 $self->dnamut->region_value($exonvalue);
		 if ($afterdist < $beforedist) {
		     $afterdist++; 		  
		     $afterdist++;
		     $self->dnamut->region_dist($afterdist);
		     #print "splice site $afterdist nt upstream!<br>";
		 } else {
		     $self->dnamut->region_dist($beforedist);
		     #print "splice site $beforedist nt downstream!<br>";
		 }
	     } else {
		 #print "first exon  : ", $firstexon->start, " - ", $firstexon->end, "<br>\n";
		 foreach $i (0..$#exons) {
		     $after=$exons[$i]->start;
		     #proximity test for intronic mutations
		     if ( ($strand == 1 and 
			   $self->mutation->prelabel >=  $before and 
			   $self->mutation->postlabel <= $after) 
			  or
			  ($strand == -1 and 
			   $self->mutation->prelabel <=  $before and 
			   $self->mutation->postlabel >= $after)  ) {
			 $self->dnamut->region('intron');
			 #$self->dnamut->region_value($i);
			 my $afterdist = $self->mutation->prelabel - $before;
			 my $beforedist =  $after - $self->mutation->postlabel;
			 my $intronvalue = $i + 1;
			 if ($afterdist < $beforedist) {
			     $afterdist++;
			     $self->dnamut->region_value($intronvalue);
			     $self->dnamut->region_dist($afterdist);
			     #print "splice site $afterdist nt upstream!<br>";
			 } else {
			     $self->dnamut->region_value($intronvalue);
			     $self->dnamut->region_dist($beforedist * -1);
			     #print "splice site $beforedist nt downstream!<br>";
			 }
			 $self->rnachange(undef);
			 last;
		     } 
		     #proximity test for exon mutations
		     elsif ( ( $strand == 1 and 
			       $exons[$i]->start <= $self->mutation->prelabel and 
			       $exons[$i]->end >= $self->mutation->postlabel) or 
			     ( $strand == -1 and 
			       $exons[$i]->start >= $self->mutation->prelabel and 
			       $exons[$i]->end <= $self->mutation->postlabel) ) {
			 $rnaAffected = 1;

			 my $afterdist = $self->mutation->prelabel - $exons[$i]->start;
			 my $beforedist =  $exons[$i]->end - $self->mutation->postlabel;
			 my $exonvalue = $i + 1;
			 $self->dnamut->region('exon');
			 if ($afterdist < $beforedist) {
			     $afterdist++;
			     $self->dnamut->region_value($exonvalue+1);
			     $self->dnamut->region_dist($afterdist);
			     #print "splice site $afterdist nt upstream!<br>";
			 } else {
			     #$beforedist;
			     $self->dnamut->region_value($exonvalue+1);
			     $self->dnamut->region_dist($beforedist * -1);
			     #print "splice site $beforedist nt downstream!<br>";
			 }
			 last;
		     }
		     $before=$exons[$i]->end;
		 }
	     }
	 }
     }
    #$self->warn("RNA not affected because change occurs inside an intron");
    #return(0); # if still not returned, then not affected, return 0
    return $rnaAffected;
}

#
# ### Creation of RNA and AA variation objects
#

    my ($self, $seqDiff) = @_;

    my $dnamut_start = $self->mutation->label - $seqDiff->offset;
    # if negative DNA positions (before ATG)
    $dnamut_start-- if $dnamut_start <= 0;
    my $dnamut_end;
    ($self->mutation->len == 0 or $self->mutation->len == 1) ?
	($dnamut_end = $dnamut_start) :
	($dnamut_end = $dnamut_start+$self->mutation->len);
    #print "start:$dnamut_start, end:$dnamut_end\n";
    my $dnamut = Bio::Variation::DNAMutation->new(-start => $dnamut_start,
						  -end => $dnamut_end,
						  );
    $dnamut->mut_number($self->mutation->issue);
    $dnamut->isMutation(1);
    my $da_m = Bio::Variation::Allele->new;
    $da_m->seq($self->mutation->seq) if $self->mutation->seq;
    $dnamut->allele_mut($da_m);
    $dnamut->add_Allele($da_m);
    # allele_ori
    my $allele_ori = $self->DNA->labelsubseq($self->mutation->prelabel,
					     undef,
					     $self->mutation->postlabel); # get seq
    chop $allele_ori; # chop the postlabel nucleotide
    $allele_ori=substr($allele_ori,1); # away the prelabel nucleotide
    my $da_o = Bio::Variation::Allele->new;
    $da_o->seq($allele_ori) if $allele_ori;
    $dnamut->allele_ori($da_o);
    ($self->mutation->len == 0) ?
	($dnamut->length($self->mutation->len)) : ($dnamut->length(CORE::length $allele_ori));
    #print " --------------- $dnamut_start -$len-  $dnamut_end -\n";
    $seqDiff->add_Variant($dnamut);
    $self->dnamut($dnamut);
    $dnamut->mut_number($self->mutation->issue);
    # setting proof
    if ($seqDiff->numbering eq "entry" or $seqDiff->numbering eq "gene") {
	 $dnamut->proof('experimental');
    } else {
	 $dnamut->proof('computed');
    }
    # how many nucleotides to store upstream and downstream of the change
    my $flanklen = $self->{'flanklen'};
    #print  `date`, " flanking sequences start\n";
    my $uplabel = $self->DNA->label(1-$flanklen,$self->mutation->prelabel); # this could be unavailable!

    my $upstreamseq;
    if ($uplabel > 0) {
	 $upstreamseq =
	     $self->DNA->labelsubseq($uplabel, undef, $self->mutation->prelabel);
    } else { # from start (less than $flanklen nucleotides)
	 $upstreamseq =
	     $self->DNA->labelsubseq($self->DNA->start, undef, $self->mutation->prelabel);
    }
    $dnamut->upStreamSeq($upstreamseq);
    my $dnstreamseq = $self->DNA->labelsubseq($self->mutation->postlabel, $flanklen);
    $dnamut->dnStreamSeq($dnstreamseq); # $flanklen or less nucleotides
    return $dnamut;
}


#
### Check if mutation propagates to RNA (and AA) level
#  
# side effect: sets intron/exon information
# returns a boolean value
#

    my ($self, $seqDiff, $dnamut) = @_;
    my ($rnapos_end, $upstreamseq, $dnstreamseq);
    my $flanklen = $self->{'flanklen'};

    ($self->mutation->len == 0) ?
	($rnapos_end = $self->mutation->transpos) :
	($rnapos_end = $self->mutation->transpos + $self->mutation->len -1);
    my $rnachange = Bio::Variation::RNAChange->new(-start => $self->mutation->transpos,
						    -end =>  $rnapos_end
						    );
    $rnachange->isMutation(1);

    # setting proof
    if ($seqDiff->numbering eq "coding") {
	 $rnachange->proof('experimental');
    } else {
	 $rnachange->proof('computed');
    }

    $seqDiff->add_Variant($rnachange);
    $self->rnachange($rnachange);
    $rnachange->DNAMutation($dnamut);
    $dnamut->RNAChange($rnachange);
    $rnachange->mut_number($self->mutation->issue);

    # setting the codon_position of the "start" nucleotide of the change
    $rnachange->codon_pos(($self->RNA->frame($self->mutation->label))+1); # codon_pos=frame+1

    my @exons=$self->RNA->all_Exons;
    $self->exons(\@exons);
    #print  `date`, " before flank, after exons. RNAObj query\n";
    # if cannot retrieve from Transcript, Transcript::upstream_seq will be used
    # before "fac7 g 65" bug discovered
    # $uplabel=$self->RNA->label(1-$flanklen,$prelabel);
    my $RNAprelabel=$self->RNA->label(-1,$self->mutation->label); # to fix fac7g65 bug
    # for the fix, all prelabel used in the next block have been changed to RNAprelabel
    my $uplabel=$self->RNA->label(1-$flanklen,$RNAprelabel);
    if ($self->RNA->valid($uplabel)) {
	 $upstreamseq = $self->RNA->labelsubseq($uplabel, undef, $RNAprelabel);
    } else {
	$upstreamseq = $self->RNA->labelsubseq($self->RNA->start, undef, $RNAprelabel)
	    if $self->RNA->valid($RNAprelabel);
	my $lacking=$flanklen-length($upstreamseq); # how many missing
	my $upstream_atg=$exons[0]->subseq(-$lacking,-1);
	$upstreamseq=$upstream_atg . $upstreamseq;
    }

    $rnachange->upStreamSeq($upstreamseq);

    # won't work OK if postlabel NOT in Transcript
    # now added RNApostlabel but this has to be /fully tested/
    # for the fix, all postlabel used in the next block have been changed to RNApostlabel
    my $RNApostlabel; # to fix fac7g64 bug
    if ($self->mutation->len == 0) {
      $RNApostlabel=$self->mutation->label;
    } else {
      my $mutlen = 1 + $self->mutation->len;
      $RNApostlabel=$self->RNA->label($mutlen,$self->mutation->label);
    }
    $dnstreamseq=$self->RNA->labelsubseq($RNApostlabel, $flanklen);
    if ($dnstreamseq eq '-1') { # if out of transcript was requested
	 my $lastexon=$exons[-1];
	 my $lastexonlength=$lastexon->length;
	 $dnstreamseq=$self->RNA->labelsubseq($RNApostlabel); # retrieves till RNAend
	 my $lacking=$flanklen-length($dnstreamseq); # how many missing
	 my $downstream_stop=$lastexon->subseq($lastexonlength+1,undef,$lacking);
	 $dnstreamseq .= $downstream_stop;
    } else {
	 $rnachange->dnStreamSeq($dnstreamseq);
    }
    # AAChange creation
    my $AAobj=$self->RNA->get_Translation;
    # storage of prelabel here, to be used in create_mut_objs_after
    my $aachange = Bio::Variation::AAChange->new(-start => $RNAprelabel
						  );
    $aachange->isMutation(1);
    $aachange->proof('computed');

    $seqDiff->add_Variant($aachange);
    $self->aachange($aachange);
    $rnachange->AAChange($aachange);
    $aachange->RNAChange($rnachange);

    $aachange->mut_number($self->mutation->issue);
#    $before_mutation{aachange}=$aachange;

    my $ra_o = Bio::Variation::Allele->new;
    $ra_o->seq($dnamut->allele_ori->seq) if $dnamut->allele_ori->seq;
    $rnachange->allele_ori($ra_o);

    $rnachange->length(CORE::length $rnachange->allele_ori->seq);

    my $ra_m = Bio::Variation::Allele->new;
    $ra_m->seq($self->mutation->seq) if $self->mutation->seq;
    $rnachange->allele_mut($ra_m);
    $rnachange->add_Allele($ra_m);

    #$rnachange->allele_mut($seq);
    $rnachange->end($rnachange->start) if $rnachange->length == 0;

    # this holds the aminoacid sequence that will be affected by the mutation
    my $aa_allele_ori=$AAobj->labelsubseq($self->mutation->label,undef, 
					  $self->mutation->lastlabel);

    my $aa_o = Bio::Variation::Allele->new;
    $aa_o->seq($aa_allele_ori) if $aa_allele_ori;
    $aachange->allele_ori($aa_o);
    #$aachange->allele_ori($aa_allele_ori);

    my $aa_length_ori = length($aa_allele_ori);
    $aachange->length($aa_length_ori); #print "==========$aa_length_ori\n";
    $aachange->end($aachange->start + $aa_length_ori - 1 );

    my ($self, $seqDiff, $dnamut) = @_;
    my $rnapos_end;
    ($self->mutation->len == 0) ?
	($rnapos_end = $self->mutation->transpos) :
	($rnapos_end = $self->mutation->transpos + $self->mutation->len -1);
    my $rnachange = Bio::Variation::RNAChange->new(-start => $self->mutation->transpos,
						    -end =>  $rnapos_end
						    );
    #my $rnachange = Bio::Variation::RNAChange->new;

    $rnachange->isMutation(1);
    my $ra_o = Bio::Variation::Allele->new;
    $ra_o->seq($dnamut->allele_ori->seq) if $dnamut->allele_ori->seq;
    $rnachange->allele_ori($ra_o);
    my $ra_m = Bio::Variation::Allele->new;
    $ra_m->seq($dnamut->allele_mut->seq) if $dnamut->allele_mut->seq;
    $rnachange->allele_mut($ra_m);
    $rnachange->add_Allele($ra_m);
    $rnachange->upStreamSeq($dnamut->upStreamSeq);
    $rnachange->dnStreamSeq($dnamut->dnStreamSeq);
    $rnachange->length($dnamut->length);
    $rnachange->mut_number($dnamut->mut_number);
    # setting proof
    if ($seqDiff->numbering eq "coding") {
	$rnachange->proof('experimental');
    } else {
	$rnachange->proof('computed');
    }

    my $dist; 
    if ($rnachange->end < 0) {
	$rnachange->region('5\'UTR');
	$dnamut->region('5\'UTR');
	my $dist = $dnamut->end ;
	$dnamut->region_dist($dist);
	$dist = $seqDiff->offset - $self->gene->maxtranscript->start + 1 + $dist;
	$rnachange->region_dist($dist);
	return if $dist < 1; # if mutation is not in mRNA 
    } else {
	$rnachange->region('3\'UTR');
	$dnamut->region('3\'UTR');
	my $dist = $dnamut->start - $seqDiff->cds_end + $seqDiff->offset;
	$dnamut->region_dist($dist);
	$dist = $seqDiff->cds_end - $self->gene->maxtranscript->end -1 + $dist;
	$rnachange->region_dist($dist);
	return if $dist > 0; # if mutation is not in mRNA 
    }
    $seqDiff->add_Variant($rnachange);
    $self->rnachange($rnachange);
    $rnachange->DNAMutation($dnamut);
    $dnamut->RNAChange($rnachange);
}

# args: reference to label changearray, reference to position changearray
# Function: take care of the creation of mutation objects, with
# information AFTER the change takes place

  my ($self,$value) = @_;
  if (defined $value) {
      if( ! $value->isa('Bio::Variation::AAChange') ) {
	  $self->throw("Is not a Bio::Variation::AAChange object but a [$value]");
	  return undef;
      }
      else {
	  $self->{'aachange'} = $value;
      }
  }
  unless (exists $self->{'aachange'}) {
      return (undef);
  } else {
      return $self->{'aachange'};
  }

    my ($self,$value) = @_;
    if( $value->isa('Bio::Liveseq::Mutation') ) {
	my $com = ref $value;
	$self->throw("Is not a Mutation object but a [$com]" );
	return undef;
    }
    if (! $value->pos) {
	$self->warn("No value for mutation position in the sequence!");
	return undef;
    }
    if (! $value->seq && ! $value->len) {
	$self->warn("Either mutated sequence or length of the deletion must be given!");
	return undef;
    }
    push(@{$self->{'mutations'}},$value);

    my ($self) = @_;

    #
    # Sanity check
    #
    unless ($self->gene) {
	$self->warn("Input object Bio::LiveSeq::Gene is not given");
	return 0;
    }
    #
    # Setting the reference sequence based on -numbering
    #
    my @transcripts=@{$self->gene->get_Transcripts};
    my $refseq; # will hold Bio::LiveSeq:Transcript object or Bio::LiveSeq::DNA

    # 'gene' eq 'coding' if reference sequence is cDNA
    $self->numbering ('coding') if $self->gene->get_DNA->alphabet eq 'rna' and $self->numbering eq 'gene';

    if ($self->numbering =~ /(coding)( )?(\d+)?/ ) {
	$self->numbering($1);
	my $transnumber = $3;
	$transnumber-- if $3; # 1 -> 0, 2 -> 1
	if ($transnumber && $transnumber >= 0 && $transnumber <= $#transcripts) {
	    $refseq=$transcripts[$transnumber];
	} else {
	    $transnumber && $self->warn("The alternative transcript number ". $transnumber+1 .
	    "- does not exist. Reverting to the 1st transcript\n");
	    $refseq=$transcripts[0];
	}
    } else {
	$refseq=$transcripts[0]->{'seq'};
    }
    #
    # Recording the state: SeqDiff object creation  ?? transcript no.??
    #
    my $seqDiff = Bio::Variation::SeqDiff->new();
    $seqDiff->alphabet($self->gene->get_DNA->alphabet);
    $seqDiff->numbering($self->numbering);
    my ($DNAobj, $RNAobj);
    if ($refseq->isa("Bio::LiveSeq::Transcript")) {
	$self->RNA($refseq);
	$self->DNA($refseq->{'seq'});
	$seqDiff->rna_ori($refseq->seq );
	$seqDiff->aa_ori($refseq->get_Translation->seq);
    } else {
	$self->DNA($refseq);
	$self->RNA($transcripts[0]);
	$seqDiff->rna_ori($self->RNA->seq);
	$seqDiff->aa_ori($self->RNA->get_Translation->seq);
    }
    $seqDiff->dna_ori($self->DNA->seq);
    # put the accession number into the SeqDiff object ID
    $seqDiff->id($self->DNA->accession_number);

    # the atg_offset takes in account that DNA object could be a subset of the
    # whole entry (via the light_weight loader)
    my $atg_label=$self->RNA->start;
    my $atg_offset=$self->DNA->position($atg_label)+($self->DNA->start)-1;
    $seqDiff->offset($atg_offset - 1);
    $self->DNA->coordinate_start($atg_label);

    my @exons = $self->RNA->all_Exons;
    $seqDiff->cds_end($exons[$#exons]->end);

    #
    # Converting mutation positions to labels
    #
    $self->warn("no mutations"), return 0 
	unless $self->_mutationpos2label($refseq, $seqDiff);

    # need to add more than one rna & aa
    #foreach $transcript (@transcripts) {
    #  $seqDiff{"ori_transcript_${i}_seq"}=$transcript->seq;
    #  $seqDiff{"ori_translation_${i}_seq"}=$transcript->get_Translation->seq;
    #}

    # do changes
    my $k;
    foreach my $mutation ($self->each_Mutation) {
	next unless $mutation->label > 0;
	$self->mutation($mutation);

	$mutation->issue(++$k);
	#
	# current position on the transcript
	#
	if ($self->numbering =~ /coding/) {
	    $mutation->transpos($mutation->pos); # transpos given by user
	} else {
	    #transpos of label / It will be 0 if mutating an intron, negative if upstream of ATG
	    $mutation->transpos($self->RNA->position($mutation->label,$atg_label));
	}
	#
	# Calculate adjacent labels based on the position on the current sequence
	#
	$mutation->prelabel($self->DNA->label(-1, $mutation->label)); # 1 before label
	if ($mutation->len == 0) {
	    $mutation->postlabel($mutation->label);
	    $mutation->lastlabel($mutation->label);
	} elsif ($mutation->len == 1) {
	    $mutation->lastlabel($mutation->label); # last nucleotide affected
	    $mutation->postlabel($self->DNA->label(2,$mutation->lastlabel)); # $len after label
	} else {
	    $mutation->lastlabel($self->DNA->label($mutation->len,$mutation->label));
	    $mutation->postlabel($self->DNA->label(2,$mutation->lastlabel));
	}
	my $dnamut = $self->_set_DNAMutation($seqDiff);
	#
	#
	#
	if ($self->_rnaAffected) {
	    $self->_set_effects($seqDiff, $dnamut);
	}
	elsif ($seqDiff->offset != 0 and $dnamut->region ne 'intron') {
	    $self->_untranslated ($seqDiff, $dnamut);
	} else {
	    #$self->warn("Mutation starts outside coding region, RNAChange object not created");
	}

	#########################################################################
	# Mutations are done here!                                              #
	$refseq->labelchange($mutation->seq, $mutation->label, $mutation->len); #
	#########################################################################

	$self->_post_mutation ($seqDiff);

	$self->dnamut(undef);
	$self->rnachange(undef);
	$self->aachange(undef);
	$self->exons(undef);
    }
    # record the final state of all three sequences
    $seqDiff->dna_mut($self->DNA->seq);
    $seqDiff->rna_mut($self->RNA->seq);
    if ($refseq->isa("Bio::LiveSeq::Transcript")) {
	$seqDiff->aa_mut($refseq->get_Translation->seq);
    } else {
	$seqDiff->aa_mut($self->RNA->get_Translation->seq);
    }

    #$seqDiff{mut_dna_seq}=$gene->get_DNA->seq;
    #my $i=1;
    #foreach $transcript (@transcripts) {
    #  $seqDiff{"mut_transcript_${i}_seq"}=$transcript->seq;
    #  $seqDiff{"mut_translation_${i}_seq"}=$transcript->get_Translation->seq;
    #}
    return $seqDiff;

    my ($self, $aln) = @_;

    #
    # Sanity checks
    #

    $self->throw("Argument is not a Bio::SimpleAlign object but a [$aln]")
	unless $aln->isa('Bio::SimpleAlign');
    $self->throw("'Pairwise alignments only, please") 
	if $aln->no_sequences != 2;

    # find out the order the two sequences are given
    my $queryseq_pos = 1; #default
    my $refseq_pos = 2;
    unless ($aln->get_seq_by_pos(1)->id eq 'QUERY') {
	carp('Query sequence has to be named QUERY') 
	    if $aln->get_seq_by_pos(2)->id ne 'QUERY';
	$queryseq_pos = 2; # alternative
	$refseq_pos = 1;
    }

    # trim the alignment
    my $start =  $aln->column_from_residue_number('QUERY', 1);
    my $end =  $aln->column_from_residue_number('QUERY', 
						$aln->get_seq_by_pos($queryseq_pos)->end );
    
    my $aln2 = $aln->slice($start, $end);

    #
    # extracting mutations
    #

    my $cs = $aln2->consensus_string(51);
    my $queryseq = $aln2->get_seq_by_pos($queryseq_pos);
    my $refseq = $aln2->get_seq_by_pos($refseq_pos);

    while ( $cs =~ /(\?+)/g) {
	# pos in local coordinates
	my $pos = pos($cs) - length($1) + 1;
	my $mutation = create_mutation($self, 
				       $refseq, 
				       $queryseq, 
				       $pos, 
				       CORE::length($1)
				       );
	# reset pos to refseq coordinates
	$pos +=  $refseq->start - 1;
	$mutation->pos($pos);

        $self->add_Mutation($mutation);
    }
    return $self->change_gene();

    my ($self, $refseq, $queryseq, $pos, $len) = @_;
    
    $self->throw("Is not a Bio::PrimarySeqI object but a [$refseq]") 
	unless $refseq->isa('Bio::PrimarySeqI');
    $self->throw("Is not a Bio::PrimarySeqI object but a [$queryseq]") 
	unless $queryseq->isa('Bio::PrimarySeqI');
    $self->throw("Position is not a positive integer but [$pos]") 
	unless $pos =~ /^\+?\d+$/;
    $self->throw("Length is not a positive integer but [$len]") 
	unless $len =~ /^\+?\d+$/;

    my $mutation;
    my $refstring = $refseq->subseq($pos, $pos + $len - 1);
    my $varstring = $queryseq->subseq($pos, $pos + $len - 1);
    
    if ($len == 1 and $refstring =~ /[^\.\-\*\?]/ and 
	$varstring  =~ /[^\.\-\*\?]/ ) { #point
	$mutation = new Bio::LiveSeq::Mutation (-seq => $varstring,
						-pos => $pos,
						);
    }
    elsif ( $refstring =~ /^[^\.\-\*\?]+$/ and 
	    $varstring  !~ /^[^\.\-\*\?]+$/ ) { # deletion
	$mutation = new Bio::LiveSeq::Mutation (-pos => $pos,
						-len => $len
						);
    }
    elsif ( $refstring !~ /^[^\.\-\*\?]+$/ and 
	    $varstring  =~ /^[^\.\-\*\?]+$/ ) { # insertion
	$mutation = new Bio::LiveSeq::Mutation (-seq => $varstring,
						-pos => $pos,
						-len => 0
						);
    } else { # complex
	$mutation = new Bio::LiveSeq::Mutation (-seq => $varstring,
						-pos => $pos,
						-len => $len
						);
    }
    
    return $mutation;

  my ($self,$value) = @_;
  if (defined $value) {
      if( ! $value->isa('Bio::Variation::DNAMutation') ) {
	  $self->throw("Is not a Bio::Variation::DNAMutation object but a [$value]");
	  return undef;
      }
      else {
	  $self->{'dnamut'} = $value;
      }
  }
  unless (exists $self->{'dnamut'}) {
      return (undef);
  } else {
      return $self->{'dnamut'};
  }

   my ($self) = @_;
   return @{$self->{'mutations'}};

  my ($self,$value) = @_;
  if (defined $value) {
      $self->{'exons'} = $value;
  }
  unless (exists $self->{'exons'}) {
      return (undef);
  } else {
      return $self->{'exons'};
  }

  my ($self,$value) = @_;
  if (defined $value) {
      if( ! $value->isa('Bio::LiveSeq::Gene') ) {
	  $self->throw("Is not a Bio::LiveSeq::Gene object but a [$value]");
	  return undef;
      }
      else {
	  $self->{'gene'} = $value;
      }
  }
  unless (exists $self->{'gene'}) {
      return (undef);
  } else {
      return $self->{'gene'};
  }

  my ($self,$value) = @_;
  if (defined $value) {
      if( ! $value->isa('Bio::LiveSeq::Mutation') ) {
	  $self->throw("Is not a Bio::LiveSeq::Mutation object but a [$value]");
	  return undef;
      }
      else {
	  $self->{'mutation'} = $value;
      }
  }
  unless (exists $self->{'mutation'}) {
      return (undef);
  } else {
      return $self->{'mutation'};
  }

    my($class,@args) = @_;
    my $self;
    $self = {};
    bless $self, $class;

    my ($gene, $numbering) =
	    $self->_rearrange([qw(GENE
				  NUMBERING
				  )],
			      @args);

    $self->{ 'mutations' } = [];

    $gene && $self->gene($gene);
    $numbering && $self->numbering($numbering);

    #class constant;
    $self->{'flanklen'} = 25;
    return $self; # success - we hope!

    my ($self,$value) = @_;
    if( defined $value) {
	if ($value =~ /(coding)( )?(\d+)?/ or $value eq 'entry' or $value eq 'gene') {
	    $self->{'numbering'} = $value;
	} else { # defaulting to 'coding'
	    $self->{'numbering'} = 'coding';
	}
    }
    unless (exists $self->{'numbering'}) {
	return 'coding';
    } else {
	return $self->{'numbering'};
    }

  my ($self,$value) = @_;
  if (defined $value) {
      if( ! $value->isa('Bio::Variation::RNAChange') ) {
	  $self->throw("Is not a Bio::Variation::RNAChange object but a [$value]");
	  return undef;
      }
      else {
	  $self->{'rnachange'} = $value;
      }
  }
  unless (exists $self->{'rnachange'}) {
      return (undef);
  } else {
      return $self->{'rnachange'};
  }

  bioperl-l@bioperl.org                         - General discussion
  http://bio.perl.org/MailList.html             - About the mailing lists

  bioperl-bugs@bio.perl.org
  http://bugzilla.bioperl.org/

  Email:  heikki@ebi.ac.uk
          insana@ebi.ac.uk, jinsana@gmx.net

  Address:

     EMBL Outstation, European Bioinformatics Institute
     Wellcome Trust Genome Campus, Hinxton
     Cambs. CB10 1SD, United Kingdom

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