view release on metacpan or  search on metacpan

bin/BS_CodonJuggler.pl  view on Meta::CPAN

{
  die "BSERROR: The start and stop coordinates do not parse\n";
}

$p{ITERATE} = $p{ITERATE} || 'chromosome';
if ($p{ITERATE} ne 'genome' && $p{ITERATE} ne 'chromosome')
{
  die "BSERROR: Argument to iterate must be 'genome' or 'chromosome'.\n";
}

################################################################################
################################# CONFIGURING ##################################
################################################################################
my $newchr = $chr->iterate(-version => $p{ITERATE});
my $GD = $newchr->GD;

$p{SWAPTYPE} = $GD->codon_change_type(-from => $p{FROM}, -to => $p{TO});
$p{MORF} = $GD->rcomplement($p{FROM});
$p{OT} = $GD->rcomplement($p{TO});
$p{NEWAA} = $GD->codontable->{$p{TO}};
$p{OLDAA} = $GD->codontable->{$p{FROM}};

my $chrseq = $newchr->sequence();
my $state = $BS->SGE
  ? farm_juggling($newchr, \%p)
  : serial_juggling($newchr, \%p);

#Do error checking
$chrseq = $newchr->sequence();
my @genes = $newchr->db->features(
  -seq_id     => $newchr->seq_id,
  -types      => 'gene',
);
foreach my $gene (@genes)
{
  my $gname = $gene->display_name;
  my $gstart = $gene->start();
  my $gend = $gene->end();
  my $glen = $gend - $gstart + 1;

  my $newseq = substr $chrseq, $gstart - 1, $glen;
  my $oldseq = substr $oldchrseq, $gstart - 1, $glen;
  if ($newseq eq $oldseq && ($state->{$gname}->[0] || $state->{$gname}->[1]))
  {
    $state->{$gname}->[2] .= ' No change in sequence;';
  }
  my $cdna = $chr->make_cDNA($gene);
  my $oldpep = $GD->translate(-sequence => $cdna);
	my $newcdna = $newchr->make_cDNA($gene);
  my $newpep = $GD->translate(-sequence => $newcdna);
  if ($newpep ne $oldpep)
  {
    $state->{$gname}->[2] .= ' Change in amino acid sequence;';
  }
}

#Do reporting
print "\nReport:\n";
foreach my $gname (sort keys %{$state})
{
  my $gid = $state->{$gname}->[3];
  my @results = @{$state->{$gname}};
  next unless($results[0] || $results[1] || $results[2]);
  print $gname, q{ : };
  if ($results[0])
  {
    my $plural = $results[0] > 1  ? 's' : q{};
    print "$results[0] $p{FROM} codon$plural changed; ";
  }
  if ($results[1])
  {
    my $plural = $results[1] > 1  ? 's' : q{};
    print "$results[1] codon$plural changed; ";
  }
  print "$results[2]" if ($results[2]);
  print "\n";
}

#Tell chromosome to write itself
$newchr->add_reason($p{EDITOR}, $p{MEMO});
$newchr->write_chromosome();

exit;

__END__

=head1 NAME

  BS_CodonJuggler.pl

=head1 VERSION

    Version 2.10

=head1 DESCRIPTION

  This utility switches any one codon to any other. By default, it will not make
    any change to a gene that will cause a translation change in an overlapping
    gene; this behavior can be overridden with the -D and -V flags, which will
    allow the utility to make nonsynonymous changes to ORFs marked dubious and
    verified, respectively.

	If a stop codon is changed to a different stop codon, the change will be
    marked "stop_retained_variant". Otherwise synonymous changes are marked
    "synonymous_codon". If a stop is changed to a non stop, it is a "stop_lost".
    If a non stop is changed to a stop it is "stop_gained"; any other change is
    a "non_synonymous_codon".

=head1 USAGE

Required arguments:

  -C,  --CHROMOSOME : The chromosome to be modified
  -E,  --EDITOR : The person responsible for the edits
  -M,  --MEMO   : Justification for the edits
  -F,  --FROM   : The codon to be replaced
  -T,  --TO     : The codon to be introduced

Optional arguments:

  --ITERATE : [genome, chromosome (def)] Which version number to increment?