Bio::LiveSeq Mutator
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvs
Summary
Bio::LiveSeq::Mutator - Package mutating LiveSequences
Package variables
No package variables defined.
Included modules
Bio::LiveSeq::Mutation
Bio::Variation::AAChange
Bio::Variation::Allele
Bio::Variation::DNAMutation
Bio::Variation::RNAChange
Bio::Variation::SeqDiff
Inherit
Bio::Root::Root
Synopsis
  # $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);
}
Description
This class mutates Bio::LiveSeq::Gene objects and returns a
Bio::Variation::SeqDiff object. Mutations are described as
Bio::LiveSeq::Mutation objects. See Bio::LiveSeq::Gene,
Bio::Variation::SeqDiff, and Bio::LiveSeq::Mutation for details.
Methods
new
No description
Code
geneDescriptionCode
numberingDescriptionCode
add_MutationDescriptionCode
each_MutationDescriptionCode
mutationDescriptionCode
DNADescriptionCode
RNADescriptionCode
dnamutDescriptionCode
rnachangeDescriptionCode
aachangeDescriptionCode
exonsDescriptionCode
change_gene_with_alignmentDescriptionCode
create_mutationDescriptionCode
change_geneDescriptionCode
_mutationpos2labelDescriptionCode
_set_DNAMutationDescriptionCode
_rnaAffected
No description
Code
_set_effectsDescriptionCode
_untranslatedDescriptionCode
_post_mutation
No description
Code
Methods description
genecode    nextTop
 Title   : gene
Usage : $mutobj = $obj->gene;
: $mutobj = $obj->gene($objref);
Function:
Returns or sets the link-reference to a Bio::LiveSeq::Gene object. If no value has ben set, it will return undef Returns : an object reference or undef Args : a Bio::LiveSeq::Gene
See Bio::LiveSeq::Gene for more information.
numberingcodeprevnextTop
 Title   : numbering
Usage : $obj->numbering();
Function:
Sets and returns coordinate system used in positioning the mutations. See change_gene for details.
Example : Returns : string Args : string (coding [transcript number] | gene | entry)
add_MutationcodeprevnextTop
 Title   : add_Mutation
Usage : $self->add_Mutation($ref)
Function: adds a Bio::LiveSeq::Mutation object
Example :
Returns :
Args : a Bio::LiveSeq::Mutation
See Bio::LiveSeq::Mutation for more information.
each_MutationcodeprevnextTop
 Title   : each_Mutation
Usage : foreach $ref ( $a->each_Mutation )
Function: gets an array of Bio::LiveSeq::Mutation objects
Example :
Returns : array of Mutations
Args :
See Bio::LiveSeq::Mutation for more information.
mutationcodeprevnextTop
 Title   : mutation
Usage : $mutobj = $obj->mutation;
: $mutobj = $obj->mutation($objref);
Function:
Returns or sets the link-reference to the current mutation object. If the value is not set, it will return undef. Internal method. Returns : an object reference or undef
DNAcodeprevnextTop
 Title   : DNA
Usage : $mutobj = $obj->DNA;
: $mutobj = $obj->DNA($objref);
Function:
Returns or sets the reference to the LiveSeq object holding the reference sequence. If there is no link, it will return undef. Internal method. Returns : an object reference or undef
RNAcodeprevnextTop
 Title   : RNA
Usage : $mutobj = $obj->RNA;
: $mutobj = $obj->RNA($objref);
Function:
Returns or sets the reference to the LiveSeq object holding the reference sequence. If the value is not set, it will return undef. Internal method. Returns : an object reference or undef
dnamutcodeprevnextTop
 Title   : dnamut
Usage : $mutobj = $obj->dnamut;
: $mutobj = $obj->dnamut($objref);
Function:
Returns or sets the reference to the current DNAMutation object. If the value is not set, it will return undef. Internal method. Returns : a Bio::Variation::DNAMutation object or undef
See Bio::Variation::DNAMutation for more information.
rnachangecodeprevnextTop
 Title   : rnachange
Usage : $mutobj = $obj->rnachange;
: $mutobj = $obj->rnachange($objref);
Function:
Returns or sets the reference to the current RNAChange object. If the value is not set, it will return undef. Internal method. Returns : a Bio::Variation::RNAChange object or undef
See Bio::Variation::RNAChange for more information.
aachangecodeprevnextTop
 Title   : aachange
Usage : $mutobj = $obj->aachange;
: $mutobj = $obj->aachange($objref);
Function:
Returns or sets the reference to the current AAChange object. If the value is not set, it will return undef. Internal method. Returns : a Bio::Variation::AAChange object or undef
See Bio::Variation::AAChange for more information.
exonscodeprevnextTop
 Title   : exons
Usage : $mutobj = $obj->exons;
: $mutobj = $obj->exons($objref);
Function:
Returns or sets the reference to a current array of Exons. If the value is not set, it will return undef. Internal method. Returns : an array of Bio::LiveSeq::Exon objects or undef
See Bio::LiveSeq::Exon for more information.
change_gene_with_alignmentcodeprevnextTop
 Title   : change_gene_with_alignment
Usage : $results=$mutate->change_gene_with_alignment($aln);
Function: Returns a Bio::Variation::SeqDiff object containing the results of the changes in the alignment. The alignment has to be pairwise and have one sequence named 'QUERY', the other one is assumed to be a part of the sequence from $gene. This method offers a shortcut to change_gene and automates the creation of Bio::LiveSeq::Mutation objects. Use it with almost identical sequnces, e.g. to locate a SNP. Args : Bio::SimpleAlign object representing a short local alignment Returns : Bio::Variation::SeqDiff object or 0 on error
See Bio::LiveSeq::Mutation, Bio::SimpleAlign, and
Bio::Variation::SeqDiff for more information.
create_mutationcodeprevnextTop
 Title   : create_mutation
Usage :
Function:
Formats sequence differences from two sequences into Bio::LiveSeq::Mutation objects which can be applied to a gene. To keep it generic, sequence arguments need not to be Bio::LocatableSeq. Coordinate change to parent sequence numbering needs to be done by the calling code. Called from change_gene_with_alignment Args : Bio::PrimarySeqI inheriting object for the reference sequence Bio::PrimarySeqI inheriting object for the query sequence integer for the start position of the local sequence difference integer for the length of the sequence difference Returns : Bio::LiveSeq::Mutation object
change_genecodeprevnextTop
 Title   : change_gene
Usage : my $mutate = Bio::LiveSeq::Mutator->new(-gene => $gene,
numbering => "coding"
);
# $mut is Bio::LiveSeq::Mutation object
$mutate->add_Mutation($mut);
my $results=$mutate->change_gene();
Function: Returns a Bio::Variation::SeqDiff object containing the results of the changes performed according to the instructions present in Mutation(s). The -numbering argument decides what molecule is being changed and what numbering scheme being used: -numbering => "entry" determines the DNA level, using the numbering from the beginning of the sequence -numbering => "coding" determines the RNA level, using the numbering from the beginning of the 1st transcript Alternative transcripts can be used by specifying "coding 2" or "coding 3" ... -numbering => "gene" determines the DNA level, using the numbering from the beginning of the 1st transcript and inluding introns. The meaning equals 'coding' if the reference molecule is cDNA. Args : Bio::LiveSeq::Gene object Bio::LiveSeq::Mutation object(s) string specifying a numbering scheme (defaults to 'coding') Returns : Bio::Variation::SeqDiff object or 0 on error
_mutationpos2labelcodeprevnextTop
 Title   : _mutationpos2label
Usage :
Function: converts mutation positions into labels
Example :
Returns : number of valid mutations
Args : LiveSeq sequence object
_set_DNAMutationcodeprevnextTop
 Title   : _set_DNAMutation
Usage :
Function:
Stores DNA level mutation attributes before mutation into Bio::Variation::DNAMutation object. Links it to SeqDiff object. Example : Returns : Bio::Variation::DNAMutation object Args : Bio::Variation::SeqDiff object
See Bio::Variation::DNAMutation and Bio::Variation::SeqDiff.
_set_effectscodeprevnextTop
 Title   : _set_effects
Usage :
Function:
Stores RNA and AA level mutation attributes before mutation into Bio::Variation::RNAChange and Bio::Variation::AAChange objects. Links them to SeqDiff object. Example : Returns : Args : Bio::Variation::SeqDiff object Bio::Variation::DNAMutation object
See Bio::Variation::RNAChange, Bio::Variation::RNAChange,
Bio::Variation::SeqDiff, and Bio::Variation::DNAMutation.
_untranslatedcodeprevnextTop
 Title   : _untranslated
Usage :
Function:
Stores RNA change attributes before mutation into Bio::Variation::RNAChange object. Links it to SeqDiff object. Example : Returns : Args : Bio::Variation::SeqDiff object Bio::Variation::DNAMutation object
See Bio::Variation::RNAChange, Bio::Variation::SeqDiff and
Bio::Variation::DNAMutation for details.
Methods code
newdescriptionprevnextTop
sub new {
    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!
}
genedescriptionprevnextTop
sub gene {
  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;
      }
      else {
	  $self->{'gene'} = $value;
      }
  }
  unless (exists $self->{'gene'}) {
      return;
  } else {
      return $self->{'gene'};
  }
}
numberingdescriptionprevnextTop
sub numbering {
    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'}; }
}
add_MutationdescriptionprevnextTop
sub add_Mutation {
    my ($self,$value) = @_;
    if( $value->isa('Bio::Liveseq::Mutation') ) {
	my $com = ref $value;
	$self->throw("Is not a Mutation object but a [$com]" );
	return;
    }
    if (! $value->pos) {
	$self->warn("No value for mutation position in the sequence!");
	return;
    }
    if (! $value->seq && ! $value->len) {
	$self->warn("Either mutated sequence or length of the deletion must be given!");
	return;
    }
    push(@{$self->{'mutations'}},$value);
}
each_MutationdescriptionprevnextTop
sub each_Mutation {
   my ($self) = @_;
   return @{$self->{'mutations'}};
}
mutationdescriptionprevnextTop
sub mutation {
  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;
      }
      else {
	  $self->{'mutation'} = $value;
      }
  }
  unless (exists $self->{'mutation'}) {
      return;
  } else {
      return $self->{'mutation'};
  }
}
DNAdescriptionprevnextTop
sub DNA {
  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;
      }
      else {
	  $self->{'DNA'} = $value;
      }
  }
  unless (exists $self->{'DNA'}) {
      return;
  } else {
      return $self->{'DNA'};
  }
}
RNAdescriptionprevnextTop
sub RNA {
  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;
      }
      else {
	  $self->{'RNA'} = $value;
      }
  }
  unless (exists $self->{'RNA'}) {
      return;
  } else {
      return $self->{'RNA'};
  }
}
dnamutdescriptionprevnextTop
sub dnamut {
  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;
      }
      else {
	  $self->{'dnamut'} = $value;
      }
  }
  unless (exists $self->{'dnamut'}) {
      return;
  } else {
      return $self->{'dnamut'};
  }
}
rnachangedescriptionprevnextTop
sub rnachange {
  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;
      }
      else {
	  $self->{'rnachange'} = $value;
      }
  }
  unless (exists $self->{'rnachange'}) {
      return;
  } else {
      return $self->{'rnachange'};
  }
}
aachangedescriptionprevnextTop
sub aachange {
  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;
      }
      else {
	  $self->{'aachange'} = $value;
      }
  }
  unless (exists $self->{'aachange'}) {
      return;
  } else {
      return $self->{'aachange'};
  }
}
exonsdescriptionprevnextTop
sub exons {
  my ($self,$value) = @_;
  if (defined $value) {
      $self->{'exons'} = $value;
  }
  unless (exists $self->{'exons'}) {
      return;
  } else {
      return $self->{'exons'};
  }
}
change_gene_with_alignmentdescriptionprevnextTop
sub change_gene_with_alignment {
    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();
}
create_mutationdescriptionprevnextTop
sub create_mutation {
    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 = Bio::LiveSeq::Mutation->new(-seq => $varstring, -pos => $pos, ); } elsif ( $refstring =~ /^[^\.\-\*\?]+$/ and $varstring !~ /^[^\.\-\*\?]+$/ ) { # deletion
$mutation = Bio::LiveSeq::Mutation->new(-pos => $pos, -len => $len ); } elsif ( $refstring !~ /^[^\.\-\*\?]+$/ and $varstring =~ /^[^\.\-\*\?]+$/ ) { # insertion
$mutation = Bio::LiveSeq::Mutation->new(-seq => $varstring, -pos => $pos, -len => 0 ); } else { # complex
$mutation = Bio::LiveSeq::Mutation->new(-seq => $varstring, -pos => $pos, -len => $len ); } return $mutation;
}
change_genedescriptionprevnextTop
sub change_gene {
    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(-verbose => $self->verbose); $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;
}
_mutationpos2labeldescriptionprevnextTop
sub _mutationpos2label {
    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
#
}
_set_DNAMutationdescriptionprevnextTop
sub _set_DNAMutation {
    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
#
}
_rnaAffecteddescriptionprevnextTop
sub _rnaAffected {
    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
#proximity test for exon mutations
elsif ( ( $strand == 1 and $exons[$i]->start < $self->mutation->prelabel and $exons[$i]->end > $self->mutation->prelabel) or ( $strand == 1 and $exons[$i]->start < $self->mutation->postlabel and $exons[$i]->end > $self->mutation->postlabel) or ( $strand == -1 and $exons[$i]->start > $self->mutation->prelabel and $exons[$i]->end < $self->mutation->prelabel) or ( $strand == -1 and $exons[$i]->start > $self->mutation->postlabel 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
#
}
_set_effectsdescriptionprevnextTop
sub _set_effects {
    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 );
}
_untranslateddescriptionprevnextTop
sub _untranslated {
    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
}
_post_mutationdescriptionprevnextTop
sub _post_mutation {
    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;
}
General documentation
FEEDBACKTop
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://bioperl.org/wiki/Mailing_lists - About the mailing lists
Support Top
Please direct usage questions or support issues to the mailing list:
bioperl-l@bioperl.org
rather than to the module maintainer directly. Many experienced and
reponsive experts will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.
Reporting BugsTop
report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution. Bug reports can be submitted via the
web:
  https://redmine.open-bio.org/projects/bioperl/
AUTHOR - Heikki Lehvaslaiho & Joseph A.L. InsanaTop
  Email:  heikki-at-bioperl-dot-org
insana@ebi.ac.uk, jinsana@gmx.net
APPENDIXTop
  The rest of the documentation details each of the object
methods. Internal methods are usually preceded with a _