File Coverage

blib/lib/Chemistry/OpenSMILES/Writer.pm
Criterion Covered Total %
statement 158 166 95.1
branch 74 84 88.1
condition 25 27 92.5
subroutine 19 20 95.0
pod 0 2 0.0
total 276 299 92.3


line stmt bran cond sub pod time code
1             package Chemistry::OpenSMILES::Writer;
2              
3 13     13   6482 use strict;
  13         40  
  13         392  
4 13     13   90 use warnings;
  13         38  
  13         394  
5              
6 13         735 use Chemistry::OpenSMILES qw(
7             is_aromatic
8             is_chiral
9             toggle_cistrans
10 13     13   974 );
  13         30  
11 13     13   1046 use Chemistry::OpenSMILES::Parser;
  13         34  
  13         374  
12 13     13   5360 use Graph::Traversal::DFS;
  13         3836  
  13         407  
13 13     13   82 use List::Util qw( any uniq );
  13         28  
  13         29033  
14              
15             # ABSTRACT: OpenSMILES format writer
16             our $VERSION = '0.8.5'; # VERSION
17              
18             require Exporter;
19             our @ISA = qw( Exporter );
20             our @EXPORT_OK = qw(
21             write_SMILES
22             );
23              
24             sub write_SMILES
25             {
26 155     155 0 31210 my( $what, $order_sub ) = @_;
27              
28 155 100       489 if( ref $what eq 'HASH' ) {
29             # subroutine will also accept and properly represent a single
30             # atom:
31 76         195 return _pre_vertex( $what );
32             }
33              
34 79 100       249 my @moieties = ref $what eq 'ARRAY' ? @$what : ( $what );
35 79         132 my @components;
36              
37 79 100       224 $order_sub = \&_order unless $order_sub;
38              
39 79         157 for my $graph (@moieties) {
40 81         140 my @symbols;
41             my %vertex_symbols;
42 81         141 my $nrings = 0;
43 81         205 my %seen_rings;
44             my @chiral;
45 81         0 my %discovered_from;
46              
47 81         140 my $rings = {};
48              
49             my $operations = {
50 423     423   68140 tree_edge => sub { my( $seen, $unseen, $self ) = @_;
51 423 50       1252 if( $vertex_symbols{$unseen} ) {
52 0         0 ( $seen, $unseen ) = ( $unseen, $seen );
53             }
54 423         1003 push @symbols, _tree_edge( $seen, $unseen, $self, $order_sub );
55 423         1573 $discovered_from{$unseen} = $seen },
56              
57 56     56   7215 non_tree_edge => sub { my @sorted = sort { $vertex_symbols{$a} <=>
58 56         231 $vertex_symbols{$b} }
59             @_[0..1];
60             $rings->{$vertex_symbols{$_[0]}}
61             {$vertex_symbols{$_[1]}} =
62             $rings->{$vertex_symbols{$_[1]}}
63 56         134 {$vertex_symbols{$_[0]}} =
64             _depict_bond( @sorted, $graph ) },
65              
66 504     504   18832 pre => sub { my( $vertex, $dfs ) = @_;
67 504 100       1213 push @chiral, $vertex if is_chiral $vertex;
68             push @symbols,
69 1586         4051 _pre_vertex( { map { $_ => $vertex->{$_} }
70 504         1763 grep { $_ ne 'chirality' }
  1628         3253  
71             keys %$vertex } );
72 504         2345 $vertex_symbols{$vertex} = $#symbols },
73              
74 504     504   108116 post => sub { push @symbols, ')' },
75 81         924 next_root => undef,
76             };
77              
78 81 50       240 if( $order_sub ) {
79             $operations->{first_root} =
80 81     81   267 sub { return $order_sub->( $_[1], $_[0]->graph ) };
  81         3450  
81             $operations->{next_successor} =
82 81     423   251 sub { return $order_sub->( $_[1], $_[0]->graph ) };
  423         76278  
83             }
84              
85 81         520 my $traversal = Graph::Traversal::DFS->new( $graph, %$operations );
86 81         41852 $traversal->dfs;
87              
88 81 100       7349 if( scalar keys %vertex_symbols != scalar $graph->vertices ) {
89 1         25 warn scalar( $graph->vertices ) - scalar( keys %vertex_symbols ) .
90             ' unreachable atom(s) detected in moiety' . "\n";
91             }
92              
93 81 50       1913 next unless @symbols;
94 81         153 pop @symbols;
95              
96             # Dealing with chirality
97 81         199 for my $atom (@chiral) {
98 42 100       239 next unless $atom->{chirality} =~ /^@@?$/;
99              
100 40         146 my @neighbours = $graph->neighbours($atom);
101 40 100 66     4468 if( scalar @neighbours < 3 || scalar @neighbours > 4 ) {
102             # TODO: process also configurations other than tetrahedral
103 1         19 warn "chirality '$atom->{chirality}' observed for atom " .
104             'with ' . scalar @neighbours . ' neighbours, can only ' .
105             'process tetrahedral chiral centers with possible ' .
106             'lone pairs' . "\n";
107 1         6 next;
108             }
109              
110 39         117 my $chirality_now = $atom->{chirality};
111 39 50       119 if( $atom->{chirality_neighbours} ) {
112 39 50       74 if( scalar @neighbours !=
113 39         110 scalar @{$atom->{chirality_neighbours}} ) {
114 0         0 warn 'number of neighbours does not match the length ' .
115             "of 'chirality_neighbours' array, cannot process " .
116             'such chiral centers' . "\n";
117 0         0 next;
118             }
119              
120 39         78 my %indices;
121 39         61 for (0..$#{$atom->{chirality_neighbours}}) {
  39         133  
122 148         224 my $pos = $_;
123 148 100 100     182 if( scalar @{$atom->{chirality_neighbours}} == 3 && $_ != 0 ) {
  148         402  
124             # Lone pair is always second in the chiral neighbours array
125 16         26 $pos++;
126             }
127 148         469 $indices{$vertex_symbols{$atom->{chirality_neighbours}[$_]}} = $pos;
128             }
129              
130 39         97 my @order_new;
131             # In newly established order, the atom from which this one
132             # is discovered (left hand side) will be the first, if any
133 39 100       116 if( $discovered_from{$atom} ) {
134             push @order_new,
135 35         113 $indices{$vertex_symbols{$discovered_from{$atom}}};
136             }
137             # Second, there will be ring bonds as they are added the
138             # first of all the neighbours
139 39 100       115 if( $rings->{$vertex_symbols{$atom}} ) {
140 2         6 push @order_new, map { $indices{$_} }
141 0         0 sort { $a <=> $b }
142 2         4 keys %{$rings->{$vertex_symbols{$atom}}};
  2         9  
143             }
144             # Finally, all neighbours are added, uniq will remove duplicates
145 148         289 push @order_new, map { $indices{$_} }
146 164         282 sort { $a <=> $b }
147 39         91 map { $vertex_symbols{$_} }
  148         399  
148             @neighbours;
149 39         294 @order_new = uniq @order_new;
150              
151 39 100       135 if( scalar @order_new == 3 ) {
152             # Accommodate the lone pair
153 8 100       24 if( $discovered_from{$atom} ) {
154 6         19 @order_new = ( $order_new[0], 1, @order_new[1..2] );
155             } else {
156 2         5 unshift @order_new, 1;
157             }
158             }
159              
160 39 100       114 if( join( '', _permutation_order( @order_new ) ) ne '0123' ) {
161 7 100       29 $chirality_now = $chirality_now eq '@' ? '@@' : '@';
162             }
163             }
164              
165 39         178 my $parser = Chemistry::OpenSMILES::Parser->new;
166 39         212 my( $graph_reparsed ) = $parser->parse( $symbols[$vertex_symbols{$atom}],
167             { raw => 1 } );
168 39         130 my( $atom_reparsed ) = $graph_reparsed->vertices;
169 39         705 $atom_reparsed->{chirality} = $chirality_now;
170 39         251 $symbols[$vertex_symbols{$atom}] =
171             write_SMILES( $atom_reparsed );
172             }
173              
174             # Adding ring numbers
175 81         502 my @ring_ids = ( 1..99, 0 );
176 81         141 my @ring_ends;
177 81         242 for my $i (0..$#symbols) {
178 1350 100       2452 if( $rings->{$i} ) {
179 54         115 for my $j (sort { $a <=> $b } keys %{$rings->{$i}}) {
  2         9  
  54         193  
180 56 100       160 next if $i > $j;
181 28 50       83 if( !@ring_ids ) {
182             # All 100 rings are open now. There is no other
183             # solution but to terminate the program.
184 0         0 die 'cannot represent more than 100 open ring' .
185             ' bonds';
186             }
187 28 50       182 $symbols[$i] .= $rings->{$i}{$j} .
188             ($ring_ids[0] < 10 ? '' : '%') .
189             $ring_ids[0];
190             $symbols[$j] .= ($rings->{$i}{$j} eq '/' ? '\\' :
191             $rings->{$i}{$j} eq '\\' ? '/' :
192 28 100       131 $rings->{$i}{$j}) .
    100          
    50          
193             ($ring_ids[0] < 10 ? '' : '%') .
194             $ring_ids[0];
195 28         49 push @{$ring_ends[$j]}, shift @ring_ids;
  28         93  
196             }
197             }
198 1350 100       2488 if( $ring_ends[$i] ) {
199             # Ring bond '0' must stay in the end
200 2770 50       4600 @ring_ids = sort { ($a == 0) - ($b == 0) || $a <=> $b }
201 28         64 (@{$ring_ends[$i]}, @ring_ids);
  28         117  
202             }
203             }
204              
205 81         2719 push @components, join '', @symbols;
206             }
207              
208 79         402 return join '.', @components;
209             }
210              
211             # DEPRECATED
212             sub write
213             {
214 0     0 0 0 &write_SMILES;
215             }
216              
217             sub _tree_edge
218             {
219 423     423   785 my( $u, $v, $self ) = @_;
220              
221 423         1061 return '(' . _depict_bond( $u, $v, $self->graph );
222             }
223              
224             sub _pre_vertex
225             {
226 580     580   1129 my( $vertex ) = @_;
227              
228 580         1149 my $atom = $vertex->{symbol};
229 580   100     2976 my $is_simple = $atom =~ /^[bcnosp]$/i ||
230             $atom =~ /^(F|Cl|Br|I|\*)$/;
231              
232 580 100       1326 if( exists $vertex->{isotope} ) {
233 3         9 $atom = $vertex->{isotope} . $atom;
234 3         8 $is_simple = 0;
235             }
236              
237 580 100       1339 if( is_chiral $vertex ) {
238 41         126 $atom .= $vertex->{chirality};
239 41         75 $is_simple = 0;
240             }
241              
242 580 100       1427 if( $vertex->{hcount} ) { # if non-zero
243 6 100       25 $atom .= 'H' . ($vertex->{hcount} == 1 ? '' : $vertex->{hcount});
244 6         12 $is_simple = 0;
245             }
246              
247 580 100       1196 if( $vertex->{charge} ) { # if non-zero
248 8 100       39 $atom .= ($vertex->{charge} > 0 ? '+' : '') . $vertex->{charge};
249 8         44 $atom =~ s/([-+])1$/$1/;
250 8         16 $is_simple = 0;
251             }
252              
253 580 50       1167 if( $vertex->{class} ) { # if non-zero
254 0         0 $atom .= ':' . $vertex->{class};
255 0         0 $is_simple = 0;
256             }
257              
258 580 100       3438 return $is_simple ? $atom : "[$atom]";
259             }
260              
261             sub _depict_bond
262             {
263 479     479   2271 my( $u, $v, $graph ) = @_;
264              
265 479 100       1279 if( !$graph->has_edge_attribute( $u, $v, 'bond' ) ) {
266 345 100 100     71063 return is_aromatic $u && is_aromatic $v ? '-' : '';
267             }
268              
269 134         27752 my $bond = $graph->get_edge_attribute( $u, $v, 'bond' );
270 134 100 100     26875 return $bond if $bond ne '/' && $bond ne '\\';
271 43 100       197 return $bond if $u->{number} < $v->{number};
272 13         46 return toggle_cistrans $bond;
273             }
274              
275             # Reorder a permutation of elements 0, 1, 2 and 3 by taking an element
276             # and moving it two places either forward or backward in the line. This
277             # subroutine is used to check whether a sign change of tetragonal
278             # chirality is required or not.
279             sub _permutation_order
280             {
281             # Safeguard against endless cycles due to undefined values
282 64 100 100 64   13541 if( (scalar @_ != 4) ||
      66        
283 247   100 247   1293 (any { !defined || !/^[0-3]$/ } @_) ||
284             (join( ',', sort @_ ) ne '0,1,2,3') ) {
285 5         53 warn '_permutation_order() accepts only permutations of numbers ' .
286             "'0', '1', '2' and '3', unexpected input received";
287 5         53 return 0..3; # Return original order
288             }
289              
290 59   100     398 while( $_[2] == 0 || $_[3] == 0 ) {
291 50         227 @_ = ( $_[0], @_[2..3], $_[1] );
292             }
293 59 100       161 if( $_[0] != 0 ) {
294 32         105 @_ = ( @_[1..2], $_[0], $_[3] );
295             }
296 59         138 while( $_[1] != 1 ) {
297 47         152 @_[1..3] = ( @_[2..3], $_[1] );
298             }
299 59         299 return @_;
300             }
301              
302             sub _order
303             {
304 364     364   1988 my( $vertices ) = @_;
305 364         1440 my @sorted = sort { $vertices->{$a}{number} <=>
306 905         1922 $vertices->{$b}{number} } keys %$vertices;
307 364         1257 return $vertices->{shift @sorted};
308             }
309              
310             1;