branch-1-6/scripts/utilities/bp_mrtrans.PLS

   1 #!perl
   2 # $Id$
   3 use strict;
   4
   5 # Author:      Jason Stajich <jason-at-bioperl-dot-org>
   6 # Description: Perl implementation of Bill Pearson's mrtrans
   7 #              to project protein alignment back into cDNA coordinates
   8 #
   9
  10 =head1 NAME
  11
  12 bp_mrtrans - implement a transformer of alignments from protein to mrna coordinates
  13
  14 =head1 SYNOPSIS
  15
  16 Usage:
  17   bp_mrtrans -i inputfile -o outputfile [-if input format] [-of output format] [-s cDNA sequence database]  [-sf cDNA sequence format] [-h]
  18
  19 =head1 DESCRIPTION
  20
  21 This script will convert a protein alignment back into a cDNA.  Loosely
  22 based on Bill Pearson's mrtrans.
  23
  24 The options are:
  25
  26    -o filename          - the output filename [default STDOUT]
  27    -of format           - output sequence format
  28                           (multiple sequence alignment)
  29                           [default phylip]
  30    -i filename          - the input filename [required]
  31    -if format           - input sequence format
  32                           (multiple sequence alignment)
  33                           [ default clustalw]
  34    -s --seqdb filename  - the cDNA sequence database file
  35    -sf --seqformat      - the cDNA seq db format (flatfile sequence format)
  36    -h                   - this help menu
  37
  38 =head1 AUTHOR
  39
  40 Jason Stajich, jason-at-bioperl-dot-org
  41
  42 =cut
  43
  44 use strict;
  45 use Bio::Align::Utilities qw(aa_to_dna_aln);
  46 use Bio::AlignIO;
  47 use Bio::SeqIO;
  48 use Getopt::Long;
  49
  50 # TODO - finish documentation,
  51 #      - add support for extra options in output alignment formats
  52 #        such as idnewline in phylip out to support Molphy input files
  53
  54 my ($iformat,$seqformat,$oformat,$seqdb,$input,$output) = ('clustalw','fasta',
  55                                                            'phylip');
  56 my ($help,$usage);
  57
  58 $usage = "usage: bp_mrtrans.pl -i prot_alignment -if align_format -o out_dna_align -of output_format -s cDNA_seqdb -sf cDNA_seqdb\n".
  59 "defaults: -if clustalw
  60           -of phylip
  61           -sf fasta\n";
  62
  63 GetOptions(
  64            'if|iformat:s'  => \$iformat,
  65            'i|input:s'     => \$input,
  66            'o|output:s'    => \$output,
  67            'of|outformat:s'=> \$oformat,
  68            's|seqdb|db:s'  => \$seqdb,
  69            'sf|seqformat:s'=> \$seqformat,
  70            'h|help'        => sub{ exec('perldoc',$0);
  71                                    exit(0)
  72                                    },
  73            );
  74
  75 $input ||= shift;
  76 $seqdb ||= shift;
  77 $output ||= shift;
  78 if( ! defined $seqdb ) {
  79     die("cannot proceed without a valid seqdb\n$usage");
  80 }
  81 if( ! defined $input ) {
  82     die("cannot proceed without an input file\n$usage");
  83 }
  84 my $indb = new Bio::SeqIO(-file => $seqdb,
  85                           -format=>$seqformat);
  86 my %seqs;
  87 while( my $seq = $indb->next_seq ) {
  88     $seqs{$seq->id} = $seq;
  89 }
  90
  91 my $in = new Bio::AlignIO(-format => $iformat,
  92                           -file   => $input);
  93 my $out = new Bio::AlignIO(-format => $oformat,
  94                            -idlength => 22,
  95                            -interleaved => 0,
  96                            defined $output ? (-file   => ">$output") : () );
  97
  98 while( my $aln = $in->next_aln ) {
  99     my $dnaaln = aa_to_dna_aln($aln,\%seqs);
 100     $dnaaln->set_displayname_flat(1);
 101     $out->write_aln($dnaaln);
 102 }
 103
 104 __END__