scripts/utilities/bp_mrtrans.pl

   1 #!perl
   2 use strict;
   3
   4 # Author:      Jason Stajich <jason-at-bioperl-dot-org>
   5 # Description: Perl implementation of Bill Pearson's mrtrans
   6 #              to project protein alignment back into cDNA coordinates
   7 #
   8
   9 =head1 NAME
  10
  11 bp_mrtrans - implement a transformer of alignments from protein to mrna coordinates
  12
  13 =head1 SYNOPSIS
  14
  15 Usage:
  16   bp_mrtrans -i inputfile -o outputfile [-if input format] [-of output format] [-s cDNA sequence database]  [-sf cDNA sequence format] [-h]
  17
  18 =head1 DESCRIPTION
  19
  20 This script will convert a protein alignment back into a cDNA.  Loosely
  21 based on Bill Pearson's mrtrans.
  22
  23 The options are:
  24
  25    -o filename          - the output filename [default STDOUT]
  26    -of format           - output sequence format
  27                           (multiple sequence alignment)
  28                           [default phylip]
  29    -i filename          - the input filename [required]
  30    -if format           - input sequence format
  31                           (multiple sequence alignment)
  32                           [ default clustalw]
  33    -s --seqdb filename  - the cDNA sequence database file
  34    -sf --seqformat      - the cDNA seq db format (flatfile sequence format)
  35    -h                   - this help menu
  36
  37 =head1 AUTHOR
  38
  39 Jason Stajich, jason-at-bioperl-dot-org
  40
  41 =cut
  42
  43 use strict;
  44 use warnings;
  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 fasta\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__