← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to scripts/utilities/bp_pairwise_kaks.pl

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/utilities/bp_pairwise_kaks.pl
 
1
#!perl
 
2
use strict;
 
3
use warnings;
 
4
# Author Jason Stajich <jason-at-bioperl-dot-org>
 
5
 
 
6
=head1 NAME
 
7
 
 
8
bp_pairwise_kaks - script to calculate pairwise Ka,Ks for a set of sequences
 
9
 
 
10
=head1 SYNOPSIS
 
11
 
 
12
bp_pairwise_kaks.PLS -i t/data/worm_fam_2785.cdna [-f fasta/genbank/embl...] [-msa tcoffee/clustal] [-kaks yn00/codeml]
 
13
 
 
14
=head1 DESCRIPTION 
 
15
 
 
16
  This script will take as input a dataset of cDNA sequences verify
 
17
 that they contain no stop codons, align them in protein space,
 
18
 project the alignment back into cDNA and estimate the Ka
 
19
 (non-synonymous) and Ks (synonymous) substitutions based on the ML
 
20
 method of Yang with the PAML package.
 
21
 
 
22
 Requires:
 
23
 * bioperl-run package
 
24
 * PAML program codeml or yn00
 
25
 * Multiple sequence alignment programs Clustalw OR T-Coffee
 
26
 
 
27
 Often there are specific specific parameters you want to run when you
 
28
 a computing Ka/Ks ratios so consider this script a starting point and
 
29
 do not rely it on for every situation.
 
30
 
 
31
=head1 FEEDBACK
 
32
 
 
33
=head2 Mailing Lists
 
34
 
 
35
User feedback is an integral part of the evolution of this and other
 
36
Bioperl modules. Send your comments and suggestions preferably to
 
37
the Bioperl mailing list.  Your participation is much appreciated.
 
38
 
 
39
  bioperl-l@bioperl.org                  - General discussion
 
40
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 
41
 
 
42
=head2 Reporting Bugs
 
43
 
 
44
Report bugs to the Bioperl bug tracking system to help us keep track
 
45
of the bugs and their resolution. Bug reports can be submitted via the
 
46
web:
 
47
 
 
48
  https://redmine.open-bio.org/projects/bioperl/
 
49
 
 
50
=head1 AUTHOR
 
51
 
 
52
 Jason Stajich jason-at-bioperl-dot-org
 
53
 
 
54
=cut
 
55
 
 
56
eval {
 
57
    # Ka/Ks estimators
 
58
    require Bio::Tools::Run::Phylo::PAML::Codeml;
 
59
    require Bio::Tools::Run::Phylo::PAML::Yn00;
 
60
    
 
61
    # Multiple Sequence Alignment programs
 
62
    require Bio::Tools::Run::Alignment::Clustalw;
 
63
    require Bio::Tools::Run::Alignment::TCoffee;
 
64
};
 
65
if( $@ ) {
 
66
    die("Must have bioperl-run pkg installed to run this script");
 
67
}
 
68
# for projecting alignments from protein to R/DNA space
 
69
use Bio::Align::Utilities qw(aa_to_dna_aln);
 
70
 
 
71
# for input of the sequence data
 
72
use Bio::SeqIO;
 
73
use Bio::AlignIO;
 
74
 
 
75
# for the command line argument parsing
 
76
use Getopt::Long;
 
77
 
 
78
my ($aln_prog, $kaks_prog,$format, $output,
 
79
    $cdna,$verbose,$help) = qw(clustalw codeml fasta);
 
80
 
 
81
GetOptions(
 
82
           'i|input:s'      => \$cdna,
 
83
           'o|output:s'     => \$output,
 
84
           'f|format:s'     => \$format,
 
85
           'msa:s'          => \$aln_prog,
 
86
           'kaks:s'         => \$kaks_prog,
 
87
           'v|verbose'      => \$verbose,
 
88
           'h|help'         => \$help,
 
89
           );
 
90
 
 
91
if( $help ) {
 
92
    exec('perldoc',$0);
 
93
    exit(0);
 
94
}
 
95
$verbose = -1 unless $verbose;
 
96
my ($aln_factory,$kaks_factory);
 
97
if( $aln_prog =~ /clus/i ) {
 
98
    $aln_factory = Bio::Tools::Run::Alignment::Clustalw->new(-verbose => $verbose);
 
99
} elsif( $aln_prog =~ /t\_?cof/i ) {
 
100
    $aln_factory = Bio::Tools::Run::Alignment::TCoffee->new(-verbose => $verbose);
 
101
} else { 
 
102
    warn("Did not provide either 'clustalw' or 'tcoffee' as alignment program names");
 
103
    exit(0);
 
104
}
 
105
unless( $aln_factory->executable ) {
 
106
    warn("Could not find the executable for $aln_prog, make sure you have installed it and have either set ".uc($aln_prog)."DIR or it is in your PATH");
 
107
    exit(0);
 
108
}
 
109
 
 
110
 
 
111
if( $kaks_prog =~ /yn00/i ) {
 
112
    $kaks_factory = Bio::Tools::Run::Phylo::PAML::Yn00->new(-verbose => $verbose);
 
113
} elsif( $kaks_prog =~ /codeml/i ) {
 
114
    # change the parameters here if you want to tweak your Codeml running!
 
115
    $kaks_factory = Bio::Tools::Run::Phylo::PAML::Codeml->new
 
116
        (-verbose => $verbose,
 
117
         -params => { 'runmode' => -2,
 
118
                      'seqtype' => 1,
 
119
                  }
 
120
         );
 
121
}
 
122
unless ( $kaks_factory->executable ) {
 
123
    warn("Could not find the executable for $kaks_prog, make sure you have installed it and you have defined PAMLDIR or it is in your PATH");
 
124
    exit(0);
 
125
}
 
126
 
 
127
unless ( $cdna && -f $cdna && -r $cdna &&  ! -z $cdna ) {
 
128
    warn("Did not specify a valid cDNA sequence file as input"); 
 
129
    exit(0);
 
130
}
 
131
 
 
132
my $seqin = new Bio::SeqIO(-file   => $cdna, 
 
133
                           -format => $format);
 
134
 
 
135
my %seqs;
 
136
my @prots;
 
137
while( my $seq = $seqin->next_seq ) {
 
138
    $seqs{$seq->display_id} = $seq;
 
139
    my $protein = $seq->translate();
 
140
    my $pseq = $protein->seq();
 
141
    
 
142
    $pseq =~ s/\*$//;
 
143
    if( $pseq =~ /\*/ ) {
 
144
        warn("provided a cDNA (".$seq->display_id.") sequence with a stop codon, PAML will choke!");
 
145
        exit(0);
 
146
    }
 
147
    # Tcoffee can't handle '*'
 
148
    $pseq =~ s/\*//g;
 
149
    $protein->seq($pseq);
 
150
    push @prots, $protein;
 
151
}
 
152
if( @prots < 2 ) {
 
153
    warn("Need at least 2 cDNA sequences to proceed");
 
154
    exit(0);
 
155
}
 
156
 
 
157
local * OUT;
 
158
if( $output ) {
 
159
    open(OUT, ">$output") || die("cannot open output $output for writing");
 
160
} else { 
 
161
    *OUT = *STDOUT;
 
162
}
 
163
 
 
164
my $aa_aln = $aln_factory->align(\@prots);
 
165
my $dna_aln = &aa_to_dna_aln($aa_aln, \%seqs);
 
166
 
 
167
my @each = $dna_aln->each_seq();
 
168
 
 
169
 
 
170
$kaks_factory->alignment($dna_aln);
 
171
 
 
172
my ($rc,$parser) = $kaks_factory->run();
 
173
if( $rc <= 0 ) { 
 
174
    warn($kaks_factory->error_string,"\n");
 
175
    exit;
 
176
}
 
177
my $result = $parser->next_result;
 
178
 
 
179
if ($result->version =~ m/3\.15/) {
 
180
        warn("This script does not work with v3.15 of PAML!  Please use 3.14 instead.");
 
181
        exit(0);
 
182
}
 
183
 
 
184
my $MLmatrix = $result->get_MLmatrix();
 
185
 
 
186
my @otus = $result->get_seqs();
 
187
 
 
188
my @pos = map { 
 
189
    my $c= 1;
 
190
    foreach my $s ( @each ) {
 
191
        last if( $s->display_id eq $_->display_id );
 
192
        $c++;
 
193
    }
 
194
    $c; 
 
195
} @otus; 
 
196
 
 
197
print OUT join("\t", qw(SEQ1 SEQ2 Ka Ks Ka/Ks PROT_PERCENTID CDNA_PERCENTID)), "\n";
 
198
for( my $i = 0; $i < (scalar @otus -1) ; $i++) {
 
199
    for( my $j = $i+1; $j < (scalar @otus); $j++ ) {
 
200
        my $sub_aa_aln = $aa_aln->select_noncont($pos[$i],$pos[$j]);
 
201
        my $sub_dna_aln = $dna_aln->select_noncont($pos[$i],$pos[$j]);
 
202
        print OUT join("\t",  
 
203
                       $otus[$i]->display_id,
 
204
                       $otus[$j]->display_id,$MLmatrix->[$i]->[$j]->{'dN'},
 
205
                       $MLmatrix->[$i]->[$j]->{'dS'},
 
206
                       $MLmatrix->[$i]->[$j]->{'omega'},
 
207
                       sprintf("%.2f",$sub_aa_aln->percentage_identity),
 
208
                       sprintf("%.2f",$sub_dna_aln->percentage_identity),
 
209
                       ), "\n";
 
210
    }
 
211
}