Toolbar

Summary

CodingSNP - DESCRIPTION of Object

Package variables

Included modules

Inherit

Bio::EnsEMBL::Root

Synopsis

Give standard usage here

Description

Describe the object here

Methods

Methods description

 Title   : exon_struct
 Usage   : $obj->exon_struct($newval)
 Function: 
 Returns : value of exon_struct
 Args    : newvalue (optional)

 Title   : gene
 Usage   : $geneobj = $obj->gene;
 Function: Returns or sets the link-reference to a Gene object.
           If there is no link, it will return undef
 Returns : an obj_ref or undef

 Title   : snp2transcript
 Usage   :
 Function:
 Example :
 Returns : 
 Args    :

 Title   : snp_struct
 Usage   : $obj->snp_struct($newval)
 Function: 
 Returns : value of snp_struct
 Args    : newvalue (optional)

 Title   : transcript
 Usage   : $transcriptobj = $obj->transcript;
 Function: Returns or sets the link-reference to a Transcript object.
           If there is no link, it will return undef
 Returns : an obj_ref or undef

Methods code

sub _calculate_gene_coordinates {

    my ($self, $snp, $rna, $aa, $aaseq) = @_;
    my $seqDiff = undef;
    $snp || $self->throw("Give a Bio::EnsEMBL::ExternalData::Variation object as an argument");
    $snp->isa('Bio::EnsEMBL::ExternalData::Variation') ||
        $self->throw("Is not a Bio::EnsEMBL::ExternalData::Variation object but a [$snp]");

    $rna->isa('Bio::EnsEMBL::Transcript') ||
        $self->throw("Is not a Bio::EnsEMBL::Transcript object but a [$rna]");
    $aa->isa('Bio::EnsEMBL::Translation') ||
        $self->throw("Is not a Bio::EnsEMBL::Translation object but a [$aaseq]");
    $aaseq->isa('Bio::PrimarySeqI') ||
        $self->throw("Is not a Bio::PrimarySeqI object but a [$aaseq]");

    #reject the SNP if it is not withing 5kb of the coding region
    return undef if $snp->start < ($aa->start - 5000);
    return undef if $snp->end > ($aa->end + 5000);
    #only real SNPs taken into account!
    return undef if $snp->start != $snp->start; 

    #
    #first create the container object
    #
    $seqDiff = Bio::Variation::SeqDiff->new();
    $seqDiff->moltype('rna');
    $seqDiff->numbering('entry');
    #$seqDiff->dna_ori(); leave undefined, 
    #we do not want to burden the object with long contig sequnce
    $seqDiff->rna_ori($rna->dna_seq->seq);
    $seqDiff->aa_ori($rna->translate); #slow opeation?
    $seqDiff->id($rna->id);
    $seqDiff->rna_id($rna->id);
    $seqDiff->offset($aa->start -1);
    $seqDiff->cds_end($aa->end - $aa->start + 1);
    
    print STDERR "AA: ".$aa->start."\t".$aa->end."\n";


##################################
#Do some kind of dirty stuff...don't want to know about it.  
##################################
    
    my $snp_ref = $self->snp_struct;
    my %snps = %$snp_ref;

    my $ex_str = $self->exon_struct;
    my %exons = %$ex_str;

    my @ex_obj = $rna->each_Exon;
    
    my @ex_array;

    foreach my $e(@ex_obj) {
	push(@ex_array,$exons{$e->id});
    }

#Here we get the position of the SNP in what we could call DNA coordinates
    my $start_pos = $snps{$snp->id}->{'start'};
    
    my $ex_st = $ex_array[0]->{'start'};

    my $dna_pos = abs($start_pos - $ex_st + 1);
    
#Lets get the position of the SNP in transcript coordinates
    my $snp_ex = $snps{$snp->id}->{'exon'};
    my $ex_rank = $exons{$snp_ex}->{'rank'};
    my $ex_array = $ex_rank -1;
    
    
    my $count = 0;
    my $tr_length = 0;
    while ($count < ($ex_rank-1)) {
	my $length = abs($ex_array[$count]->{'end'} - $ex_array[$count]->{'start'} + 1);
	$tr_length = $tr_length + $length;
	$count++;
    }
    my $length1 =  abs($ex_array[$ex_array]->{'start'} - $start_pos + 1);
    $tr_length = $tr_length + $length1;
    #print STDERR "TR POS: $tr_length\n";
####################################
    
    #
    # DNA level
    #

    my $dna_start = $dna_pos - $seqDiff->offset;
    #my $dna_start = $snps{$snp->id}->{'start'} - $seqDiff->offset;
    $dna_start < 1 && $dna_start--; # no 0 in the coordinate system!

    my $dnamut = Bio::Variation::DNAMutation->new
	(-start => ($dna_start),
	 -end => ($dna_start),
	 );
    
    $dnamut->mut_number(1);
    $dnamut->proof('experimental'); # given coordianate system
    #$dnamut->isMutation(0);
    
    my (@alleles) = split (/\|/, $snp->alleles); 
    my $dnaA1 = Bio::Variation::Allele->new;
    my $a1 = shift @alleles;
    $dnamut->length(CORE::length $a1);
    $dnaA1->seq($a1);
    $dnamut->allele_ori($dnaA1);
    foreach my $alleleseq (@alleles) {
	my $A2 = Bio::Variation::Allele->new;
	$A2->seq($alleleseq);
	$dnamut->add_Allele($A2);
    }

    $seqDiff->add_Variant($dnamut);
    foreach my $link ($snp->each_DBLink) {
	$dnamut->add_DBLink;
    }
    
#Now it would be nice to get the postion of the SNP in transcript coordinates...joy of conversions...

    $dnamut->upStreamSeq
	(lc substr($rna->dna_seq->seq, $tr_length -25, 25));
    $dnamut->dnStreamSeq
	(lc substr($rna->dna_seq->seq, $tr_length +1, 25));


    my $ref_allele;# = lc substr($self->contig->primary_seq->seq, $snp->start, 1 );
    $dnaA1->seq ne $ref_allele &&
	$self->warn("Found DNA ref allele: $ref_allele!");

#what if the SNP is closer to either end than 25 nt?

    #Lets say that we only have exonic SNPs
		
	
    $dnamut->region('exon');

#WHAT SHOULD I DO WITH THAT?????

    #$dnamut->region_value($exon->id. " ($exoncount)");
		

    #
    # RNA level
    #

    my $rnachange = undef;
    #if ($dnamut->region =~ /UTR/ or $dnamut->region eq 'exon' ) {
    if ($dnamut->region eq 'exon' ) { # mRNA affected

#Check exactly what its doing...
	$seqDiff->rna_offset($rna->start_exon->phase -1);

	#print STDERR "OFFSET: ".$seqDiff->rna_offset."\n";
	# new method into Bio::Ensembl::Transcript
	# we shouldn't need this method here
	#my $rna_pos = $rna->rna_pos($snp->start); 
	
	my $rna_pos = $tr_length;
	#print STDERR "RNA $rna_pos\n";
	#print STDERR "RNAPOS: $pos\n";
	#$rnachange = Bio::Variation::RNAChange->new
	#    (-start => $rna_pos - $seqDiff->rna_offset, 
	#     -end =>  $rna_pos - $seqDiff->rna_offset,
	#     );
	#print STDERR "HERE1\n";
	$rnachange = Bio::Variation::RNAChange->new
	    (-start => $rna_pos - $seqDiff->rna_offset, 
	     -end =>  $rna_pos - $seqDiff->rna_offset,
	     );
	#print STDERR "HERE2\n";
	$rnachange->length(1);
	$rnachange->mut_number(1);
	$seqDiff->add_Variant($rnachange);
	$dnamut->RNAChange($rnachange);
	$rnachange->DNAMutation($dnamut);
	$rnachange->proof('computed');
	$rnachange->region('coding');
	

	$rnachange->allele_ori($dnaA1);
	foreach my $alleleseq (@alleles) {
	    my $A2 = Bio::Variation::Allele->new;
	    $A2->seq($alleleseq);
	    $rnachange->add_Allele($A2);
	}
	$rnachange->codon_pos(($rnachange->start -1 )% 3 +1);
	if ($rna_pos < 25 ) {
	    $rnachange->upStreamSeq($dnamut->upStreamSeq);
	} else {
	    $rnachange->upStreamSeq
		(lc substr($seqDiff->rna_ori, $rna_pos -25, 25));
	}    
	if ($rna_pos > $seqDiff->cds_end - 25 ) {
	    $rnachange->dnStreamSeq($dnamut->dnStreamSeq);
	} else {
	    $rnachange->dnStreamSeq
		(lc substr($seqDiff->rna_ori, $rna_pos + 1, 25));
	}
	my $ref_allele = lc substr($seqDiff->rna_ori, $rna_pos,1);
	$dnaA1->seq ne  $ref_allele && 
	    $self->warn("Found RNA ref allele: $ref_allele!");
    }    

    #print STDERR $rnachange->upStreamSeq, "\n";

    #
    # Protein level
    #
    
     if ($dnamut->region eq 'exon' ) { # coding region affected
	
	my $aachange = Bio::Variation::AAChange->new
	    (-start => (int($rnachange->start / 3 + 1))
	     );
	
	$aachange->end($aachange->start);
	$aachange->proof('computed');
	$seqDiff->add_Variant($aachange);
	$rnachange->AAChange($aachange);
	$aachange->RNAChange($rnachange);
	$aachange->mut_number(1);
	
	my $ct = new Bio::Tools::CodonTable;
	my $aa_allele_ori = $ct->translate($rnachange->codon_ori);

	print STDERR "ORI: ".$rnachange->codon_ori, "\n";

	my $aa_o = Bio::Variation::Allele->new;
	$aa_o->seq($aa_allele_ori) if $aa_allele_ori;
	$aachange->allele_ori($aa_o);
	
	my $aa_allele_mut = $ct->translate($rnachange->codon_mut);
	my $aa_m = Bio::Variation::Allele->new;
	$aa_m->seq($aa_allele_mut) if $aa_allele_mut;
	$aachange->add_Allele($aa_m);
	
	my $aa_length_ori = CORE::length($aa_allele_ori);
	$aachange->length($aa_length_ori);
	$aachange->end($aachange->start + $aa_length_ori - 1 );
	
	#terminator codon?
	
	my $ref_allele =  substr($aaseq->seq, $aachange->start -1 ,1);
	$aa_allele_ori ne  $ref_allele && 
	    $self->warn ("Found AA ref allele '$ref_allele' when expected to see $aa_allele_ori");
    }
    return $seqDiff;

}

sub exon_struct {

   my $obj = shift;
   if( @_ ) {
      my $value = shift;
      $obj->{'exon_struct'} = $value;
    }
    return $obj->{'exon_struct'};

}

sub gene {

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

}

sub new {

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

    my ($gene, $transcript, $exonstruct, $snpstruct) =
	    $self->_rearrange([qw(GENE
				  TRANSCRIPT
				  EXON_STRUCT
				  SNP_STRUCT
				)],@args);
    $gene && $self->gene($gene);
    $transcript && $self->transcript($transcript);
    $exonstruct && $self->exon_struct($exonstruct);
    $snpstruct && $self->snp_struct($snpstruct);
    return $self; # success - we hope!

}

sub snp2gene {

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

    $self->gene || $self->throw("Set gene to a Bio::EnsEMBL::Gene object");
    
    my @seqDiffs = ();
    foreach my $trans ($self->gene->each_Transcript) {
	#print STDERR $trans->id,"\n";
	#print STDERR $rna->dna_seq->seq, "\n";
	my $aa = $trans->translation; #Bio::Ensembl::Translation object
	#print $aa->id, ", ", $aa->start, ", ", $aa->end, "\n";
	my $aaseq = $trans->translate;
	#print STDERR $aaseq->seq, "\n";
	
	$seqDiff = $self->_calculate_gene_coordinates($snp, $trans, $aa, $aaseq);
	$seqDiff && push @seqDiffs, $seqDiff;
    }
    return @seqDiffs;

}

sub snp2transcript {

   my ($self,$snp) = @_;
   my $seqDiff = undef;
   $self->transcript || $self->throw("Set transcript to a Bio::EnsEMBL::Transcript object");
   
   my $trans = $self->transcript;

   my $aa = $trans->translation; #Bio::Ensembl::Translation object
   use Data::Dumper;
   print Dumper($aa);
   my $aaseq = $trans->translate;
   $seqDiff = $self->_calculate_gene_coordinates($snp, $trans, $aa, $aaseq);
   return $seqDiff;

}

sub snp_struct {

   my $obj = shift;
   if( @_ ) {
      my $value = shift;
      $obj->{'snp_struct'} = $value;
    }
    return $obj->{'snp_struct'};

}

sub transcript {

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

}

General documentation

Describe contact details here

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