File Coverage

blib/lib/HackaMol/X/Vina.pm
Criterion Covered Total %
statement 66 90 73.3
branch 15 22 68.1
condition n/a
subroutine 16 20 80.0
pod 3 5 60.0
total 100 137 72.9


line stmt bran cond sub pod time code
1             package HackaMol::X::Vina;
2             $HackaMol::X::Vina::VERSION = '0.012';
3             #ABSTRACT: HackaMol extension for running Autodock Vina
4 2     2   210399 use Moose;
  2         915441  
  2         19  
5 2     2   17396 use MooseX::StrictConstructor;
  2         53843  
  2         14  
6 2     2   19296 use Moose::Util::TypeConstraints;
  2         12  
  2         20  
7 2     2   6427 use Math::Vector::Real;
  2         32183  
  2         244  
8 2     2   1452 use MooseX::Types::Path::Tiny qw(AbsPath) ;
  2         529938  
  2         27  
9 2     2   6872 use HackaMol; # for building molecules
  2         8164153  
  2         85  
10 2     2   1346 use File::chdir;
  2         5310  
  2         244  
11 2     2   17 use namespace::autoclean;
  2         3  
  2         21  
12 2     2   147 use Carp;
  2         2  
  2         1760  
13              
14             with qw(HackaMol::X::Roles::ExtensionRole);
15              
16             has $_ => (
17             is => 'rw',
18             isa => AbsPath,
19             predicate => "has_$_",
20             required => 1,
21             coerce => 1,
22             ) foreach ( qw( receptor ligand ) );
23              
24             has 'save_mol' => (
25             is => 'rw',
26             isa => 'Bool',
27             default => 0,
28             );
29            
30              
31             has $_ => (
32             is => 'rw',
33             isa => 'Num',
34             predicate => "has_$_",
35             ) foreach qw(center_x center_y center_z size_x size_y size_z);
36              
37             has 'num_modes' => (
38             is => 'rw',
39             isa => 'Int',
40             predicate => "has_num_modes",
41             default => 1,
42             lazy => 1,
43             );
44              
45             has $_ => (
46             is => 'rw',
47             isa => 'Int',
48             predicate => "has_$_",
49             ) foreach qw(energy_range exhaustiveness seed cpu);
50              
51              
52             has 'center' => (
53             is => 'rw',
54             isa => 'Math::Vector::Real',
55             predicate => "has_center",
56             trigger => \&_set_center,
57             );
58              
59             has 'size' => (
60             is => 'rw',
61             isa => 'Math::Vector::Real',
62             predicate => "has_size",
63             trigger => \&_set_size,
64             );
65              
66             sub BUILD {
67 7     7 0 11 my $self = shift;
68              
69 7 100       236 if ( $self->has_scratch ) {
70 5 100       158 $self->scratch->mkpath unless ( $self->scratch->exists );
71             }
72              
73             # build in some defaults
74 7 100       1785 $self->in_fn("conf.txt") unless ($self->has_in_fn);
75 7 100       867 $self->exe($ENV{"HOME"}."/bin/vina") unless $self->has_exe;
76              
77 7 100       317 unless ( $self->has_out_fn ) {
78 5         187 my $outlig = $self->ligand->basename;
79 5         80 $outlig =~ s/\.pdbqt/\_out\.pdbqt/;
80 5         152 $self->out_fn($outlig);
81             }
82              
83 7 100       694 unless ( $self->has_command ) {
84 5         41 my $cmd = $self->build_command;
85 5         153 $self->command($cmd);
86             }
87              
88 7         277 return;
89             }
90              
91             sub _set_center {
92 3     3   7 my ( $self, $center, $old_center ) = @_;
93 3         137 $self->center_x( $center->[0] );
94 3         82 $self->center_y( $center->[1] );
95 3         91 $self->center_z( $center->[2] );
96             }
97              
98             sub _set_size {
99 3     3   6 my ( $self, $size, $old_size ) = @_;
100 3         79 $self->size_x( $size->[0] );
101 3         101 $self->size_y( $size->[1] );
102 3         87 $self->size_z( $size->[2] );
103             }
104              
105             #required methods
106             sub build_command {
107 7     7 0 1388 my $self = shift;
108 7         12 my $cmd;
109 7         238 $cmd = $self->exe;
110 7         230 $cmd .= " --config " . $self->in_fn->stringify;
111              
112             # we always capture output
113 7         143 return $cmd;
114             }
115              
116             sub _build_map_in {
117             # this builds the default behavior, can be set anew via new
118 1     1   18 return sub { return ( shift->write_input ) };
  1     1   288  
119             }
120              
121             sub _build_map_out {
122             # this builds the default behavior, can be set anew via new
123             my $sub_cr = sub {
124 0     0   0 my $self = shift;
125 0         0 my $qr = qr/^\s+\d+\s+(-*\d+\.\d)/;
126 0         0 my ( $stdout, $sterr ) = $self->capture_sys_command;
127 0         0 my @be = map { m/$qr/; $1 }
  0         0  
  0         0  
128 0         0 grep { m/$qr/ }
129             split( "\n", $stdout );
130 0         0 return (@be);
131 0     0   0 };
132 0         0 return $sub_cr;
133             }
134              
135             sub dock {
136 0     0 1 0 my $self = shift;
137 0         0 my $num_modes = shift;
138 0 0       0 $self->num_modes($num_modes) if defined($num_modes);
139 0         0 $self->map_input;
140 0         0 return $self->map_output;
141             }
142              
143             sub dock_mol {
144             # want this to return configurations of the molecule
145 0     0 1 0 my $self = shift;
146 0         0 my $num_modes = shift;
147 0 0       0 $self->num_modes($num_modes) if defined($num_modes);
148 0         0 $self->map_input;
149 0 0       0 local $CWD = $self->scratch if ( $self->has_scratch );
150 0         0 my @bes = $self->map_output; # this is fragile... broken if map_out changed...
151 0         0 my $mol = HackaMol -> new(hush_read => 1)
152             -> read_file_mol($self->out_fn->stringify);
153 0         0 $mol->push_score(@bes);
154 0         0 return ($mol);
155             }
156              
157             sub write_input {
158 1     1 1 2 my $self = shift;
159 1         2 my $input;
160 1         40 $input .= sprintf( "%-15s = %-55s\n", 'out', $self->out_fn->stringify );
161 1 50       90 $input .= sprintf( "%-15s = %-55s\n", 'log', $self->log_fn->stringify )
162             if $self->has_log_fn;
163 1         12 foreach my $cond (
164             qw(receptor ligand cpu num_modes energy_range exhaustiveness seed))
165             {
166 7         22 my $condition = "has_$cond";
167 7 100       356 $input .= sprintf( "%-15s = %-55s\n", $cond, $self->$cond )
168             if $self->$condition;
169             }
170 1         4 foreach my $metric (qw(center_x center_y center_z size_x size_y size_z)) {
171 6         368 $input .= sprintf( "%-15s = %-55s\n", $metric, $self->$metric );
172             }
173 1         43 $self->in_fn->spew($input);
174 1         842 return ($input);
175             }
176              
177             __PACKAGE__->meta->make_immutable;
178              
179             1;
180              
181             __END__
182              
183             =pod
184              
185             =head1 NAME
186              
187             HackaMol::X::Vina - HackaMol extension for running Autodock Vina
188              
189             =head1 VERSION
190              
191             version 0.012
192              
193             =head1 SYNOPSIS
194              
195             use Modern::Perl;
196             use HackaMol;
197             use HackaMol::X::Vina;
198             use Math::Vector::Real;
199              
200             my $receptor = "receptor.pdbqt";
201             my $ligand = "lig.pdbqt",
202             my $rmol = HackaMol -> new( hush_read=>1 ) -> read_file_mol( $receptor );
203             my $lmol = HackaMol -> new( hush_read=>1 ) -> read_file_mol( $ligand );
204             my $fh = $lmol->print_pdb("lig_out.pdb");
205              
206             my @centers = map {$_ -> xyz}
207             grep {$_ -> name eq "OH" }
208             grep {$_ -> resname eq "TYR"} $rmol -> all_atoms;
209              
210             foreach my $center ( @centers ){
211              
212             my $vina = HackaMol::X::Vina -> new(
213             receptor => $receptor,
214             ligand => $ligand,
215             center => $center,
216             size => V( 20, 20, 20 ),
217             cpu => 4,
218             exhaustiveness => 12,
219             exe => '~/bin/vina',
220             scratch => 'tmp',
221             );
222              
223             my $mol = $vina->dock_mol(3); # fill mol with 3 binding configurations
224              
225             printf ("Score: %6.1f\n", $mol->get_score($_) ) foreach (0 .. $mol->tmax);
226              
227             $mol->print_pdb_ts([0 .. $mol->tmax], $fh);
228              
229             }
230              
231             $_->segid("hgca") foreach $rmol->all_atoms; #for vmd rendering cartoons.. etc
232             $rmol->print_pdb("receptor.pdb");
233              
234             =head1 DESCRIPTION
235              
236             HackaMol::X::Vina provides an interface to AutoDock Vina, which is a widely used program for docking small molecules
237             (ligands) into biological molecules (receptors). This class provides methods for writing configuration files and for
238             processing output. The input/output associated with running Vina is pretty simple, but there is still a fair amount of
239             scripting required to apply the program to virtual drug-screens that often involve sets of around 100,000 ligands, several
240             sites (centers) within a given receptor, which may also have multiple configurations. The goal of this interface is to reduce
241             the amount of scripting needed to set up massive drug screens, provide flexibility in analysis/application, and improve
242             control of what is written into files that can quickly accumulate. For example, the synopsis docks a ligand into a
243             receptor for a collection of centers located at the hydroxy group of tyrosine residues; there are a multitude of binding
244             site prediction software that can be used to provide a collection of centers. Loops over ligands, receptors, centers are
245             straightforward to implement, but large screens on a computer cluster will require splitting the loops into chunks that
246             can be spread across the queueing system. See the examples.
247              
248             This class does not include the AutoDock Vina program, which is
249             L<released under a very permissive Apache license|http://vina.scripps.edu/manual.html#license>, with few
250             restrictions on commercial or non-commercial use, or on the derivative works, such is this. Follow these
251             L<instructions | http://vina.scripps.edu/manual.html#installation> to acquire the program. Most importantly, if
252             you use this interface effectively, please be sure to cite AutoDock Vina in your work:
253              
254             O. Trott, A. J. Olson, AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading, Journal of Computational Chemistry 31 (2010) 455-461
255              
256             Since HackaMol has no pdbqt writing capabilities (yet, HackaMol can read pdbqt files; hush_read=>1 recommended, see
257             synopsis), the user is required to provide those files. L<OpenBabel| http://openbabel.org/wiki/Main_Page> and L<MGLTools| http://mgltools.scripps.edu> are popular
258             and effective. This is still a work in progress and the API may change. Documentation will improve as API
259             gets more stable... comments/contributions welcome! The automated testing reported on metacpan will likely give a bunch
260             of fails until I have time to figure out how to skip tests calling on the vina program to run.
261              
262             =head1 METHODS
263              
264             =head2 write_input
265              
266             This method takes no arguments; it returns, as a scalar, the input constructed from attributes. This method is called by map_input method via the map_in attribute to write the configuration file for running Vina.
267              
268             =head2 map_input
269              
270             provided by L<HackaMol::X::Roles::ExtensionRole>. Writes the configuration file for Vina. See dock and dock_mol methods.
271              
272             =head2 map_output
273              
274             provided by L<HackaMol::X::Roles::ExtensionRole>. By default, this method returns the docking scores as an array.
275              
276             =head2 dock_mol
277              
278             this method takes the number of binding modes (Integer) as an argument (Int). The argument is optional, and the num_modes attribute is rewritten if passed. This method calls the map_input and map_output methods for preparing and running Vina. It loads the resulting pdbqt and scores into a L<HackaMol::Molecule> object. The scores are stored into the score attribute provided by the L<HackaMol::QmMolRole>. See the synopsis for an example.
279              
280             =head2 dock
281              
282             this method is similar to dock_mol, but returns only the scores.
283              
284             =head1 ATTRIBUTES
285              
286             =head2 mol
287              
288             isa L<HackaMol::Molecule> object that is 'ro' and provided by L<HackaMol::X::Roles::ExtensionRole>.
289              
290             =head2 map_in map_out
291              
292             these attributes are 'ro' CodeRefs that can be adjusted in a given instance of a class. These are provided by L<HackaMol::X::Roles::ExtensionRole>. Setting the map_in and map_out attributes are for advanced use. Defaults are provided that are used in the map_input and map_output methods.
293              
294             =head2 receptor ligand
295              
296             these attributes are 'rw' and coerced into L<Path::Tiny> objects using the AbsPath type provided by L<MooseX::Types::Path::Tiny>. Thus, setting the receptor or ligand attributes with a string will store the entire path to the file, which
297             is provided to Vina via the input configuration file. The receptor and ligand attributes typically point to pdbqt
298             files used for running the docking calculations.
299              
300             =head2 save_mol
301              
302             this attribute isa 'Bool' that is 'rw'.
303              
304             =head2 center
305              
306             this attribute isa Math::Vector::Real object that is 'rw'. This attribute comes with a trigger that writes the
307             center_x, center_y, and center_z attributes that are used in Vina configuration files.
308              
309             =head2 center_x center_y center_z
310              
311             this attribute isa Num that is 'rw'. These attributes provide the center for the box that (with size_x, size_y, size_z) define the docking space searched by Vina. Using the center attribute may be more convenient since it has the same
312             type as the coordinates in atoms. See the synopsis.
313              
314             =head2 size_x size_y size_z
315              
316             this attribute isa Num that is 'rw'. These attributes provide the edgelengths of the the box that (with center_x,
317             center_y, center_z) define the docking space searched by Vina.
318              
319             =head2 num_modes
320              
321             this attribute isa Int that is 'rw'. It provides the requested number of binding modes (ligand configurations) for
322             Vina via the configuration file. Vina may return a smaller number of configurations depending on energy_range
323             or other factors (that need documentation).
324              
325             =head2 energy_range
326              
327             this attribute isa Int that is 'rw'. In kcal/mol, provides a window for the number of configurations to return.
328              
329             =head2 exhaustiveness
330              
331             this attribute isa Int that is 'rw'. The higher the number the more time Vina will take looking for optimal
332             docking configurations.
333              
334             =head2 cpu
335              
336             this attribute isa Int that is 'rw'. By default Vina will try to use all the cores available. Setting this
337             attribute will limit the number of cores used by Vina.
338              
339             =head2 scratch
340              
341             this attribute isa L<Path::Tiny> that is 'ro' and provided by L<HackaMol::PathRole>. Setting this attribute return a
342             Path::Tiny object with absolute path that will be created if needed and then used for
343             all Vina calculations to be run.
344              
345             =head2 in_fn
346              
347             this attribute isa L<Path::Tiny> that is 'rw' and provided by L<HackaMol::PathRole>. The default is set to conf.txt
348             when the object is built using the new method. If many instances of Vina will be running at the same time in the
349             same directory, this conf.txt will need to be unique for each one!!! The same applies to out_fn which is
350             described next.
351              
352             =head2 out_fn
353              
354             this attribute isa L<Path::Tiny> that is 'rw' and provided by L<HackaMol::PathRole>. The default is set to a value
355             derived from the the basename of the ligand attribute. i.e. out_fn is set to lig_out.pdbqt from
356             /some/big/path/lig.pdbqt. The Vina default behavior is to write to /some/big/path/lig_out.pdbqt, is usually not
357             wanted (by me anyway); thus, the default is always set and written to the configuration file. If many instances
358             of Vina will be running at the same time in the same directory, the output will need to be unique for each one as
359             described above.
360              
361             =head1 SEE ALSO
362              
363             =over 4
364              
365             =item *
366              
367             L<HackaMol>
368              
369             =item *
370              
371             L<HackaMol::X::Roles::ExtensionRole>
372              
373             =item *
374              
375             L<HackaMol::X::Calculator>
376              
377             =item *
378              
379             L<PBS::Client>
380              
381             =item *
382              
383             L<Vina | http://vina.scripps.edu>
384              
385             =item *
386              
387             L<MGLTools | http://mgltools.scripps.edu>
388              
389             =item *
390              
391             L<Open Babel | http://openbabel.org>
392              
393             =back
394              
395             =head1 EXTENDS
396              
397             =over 4
398              
399             =item * L<Moose::Object>
400              
401             =back
402              
403             =head1 CONSUMES
404              
405             =over 4
406              
407             =item * L<HackaMol::Roles::ExeRole>
408              
409             =item * L<HackaMol::Roles::ExeRole|HackaMol::Roles::PathRole>
410              
411             =item * L<HackaMol::Roles::PathRole>
412              
413             =item * L<HackaMol::X::Roles::ExtensionRole>
414              
415             =back
416              
417             =head1 AUTHOR
418              
419             Demian Riccardi <demianriccardi@gmail.com>
420              
421             =head1 COPYRIGHT AND LICENSE
422              
423             This software is copyright (c) 2015 by Demian Riccardi.
424              
425             This is free software; you can redistribute it and/or modify it under
426             the same terms as the Perl 5 programming language system itself.
427              
428             =cut