File Coverage

blib/lib/HackaMol/Roles/ReadPdbxRole.pm
Criterion Covered Total %
statement 48 205 23.4
branch 13 108 12.0
condition 0 9 0.0
subroutine 7 10 70.0
pod 0 2 0.0
total 68 334 20.3


line stmt bran cond sub pod time code
1             $HackaMol::Roles::ReadPdbxRole::VERSION = '0.052';
2             # ABSTRACT: parse PDBx/mmCIF files
3             use Moose::Role;
4 11     11   6545 use HackaMol::PeriodicTable qw(_element_name _trim _qstring_num);
  11         29  
  11         83  
5 11     11   50997 use Math::Vector::Real;
  11         31  
  11         701  
6 11     11   66 use Carp;
  11         30  
  11         484  
7 11     11   80  
  11         22  
  11         29575  
8             requires 'readline_func';
9              
10             # EXPERIMENTAL: use with caution, getting to work, then need to refactor to min DRY
11             my $self = shift;
12             my $fh = shift;
13 0     0 0 0 my $info = shift;
14 0         0 my $in_loop = 0;
15 0         0 my $atom_site_position = 0;
16 0         0  
17 0         0 while (<$fh>) {
18              
19 0         0 # set loop flag if loop is open
20             if (/^loop_/) {
21             $in_loop = 1;
22 0 0       0 }
23 0         0 if ( /^#/ && $in_loop ) {
24             $in_loop = 0;
25 0 0 0     0 }
26 0         0  
27             if (/^_pdbx_database_status.recvd_initial_deposition_date\s+(\S+)/) {
28             $info->{deposition_date} = $1;
29 0 0       0 }
30 0         0 if (/^_pdbx_audit_revision_history.revision_date\s*(\S+)?/) {
31             my $revision_date = $1;
32 0 0       0 if ($in_loop) {
33 0         0 while ( my $line = <$fh> ) {
34 0 0       0  
35 0         0 if ( $line =~ /^\#/ ) {
36             $in_loop = 0;
37 0 0       0 last;
38 0         0 }
39 0         0  
40             chomp($line);
41              
42 0         0 ($revision_date) = $line =~ /.+(\d{4}\-\d{2}\-\d{2})/;
43             }
44 0         0 }
45             $info->{last_revision_date} = $revision_date;
46             }
47 0         0 if (/^_pdbx_struct_assembly_gen.asym_id_list\s*(\S+)?/){
48             my @asym_list;
49 0 0       0 if ($1){
50 0         0 @asym_list = split /\s*,\s*/, $1;
51 0 0       0 }
52 0         0 else {
53             while (my $line = <$fh>){
54             chomp($line);
55 0         0 last if $line =~ /^(;\s*$|#)/;
56 0         0 if( $line =~ /^;/ ){
57 0 0       0 @asym_list = split /\s*,\s*/, $line;
58 0 0       0 }
59 0         0 @asym_list = split /\s*,\s*/, $line;
60             }
61 0         0  
62             }
63             $info->{auth_asym_id_list}{$_} = 0 foreach @asym_list; # use hash to connect with entities?
64             }
65 0         0 # this works for 3j7p and others but not using for now because not needed by me
66             # if (/^_struct_conn.id\s*$/) {
67             # die "assumed loop exception" unless $in_loop;
68             # my @labels = ("_struct_conn.id");
69              
70             # while ( my $line = <$fh> ) {
71             # if ( $line =~ /^(_struct\S+)/ ) {
72             # push @labels, $1;
73             # next;
74             # }
75             # if ( $line =~ /^#/ ) {
76             # $in_loop = 0;
77             # last;
78             # }
79              
80             # chomp($line);
81              
82             # # watch out for things like 'C-U MISPAIR' in 3j7p
83             # $line = _quote_hack($line);
84             # my @tokens = split /\s+/, $line;
85              
86             # # suck up lines until @tokens == @lables
87             # while ( scalar(@tokens) < scalar(@labels) ) {
88             # my $next_line = <$fh>;
89             # chomp($next_line);
90             # $next_line = _quote_hack($next_line);
91             # push @tokens, split /\s+/, $next_line;
92             # }
93             # die "too many tokens! [@labels] [@tokens] " if scalar(@tokens) != scalar(@labels);
94             # next unless ($tokens[0] =~ /disulf/);
95             # my $connect;
96             # foreach my $i ( 1 .. $#labels ) {
97             # next if ( $tokens[$i] eq '.' || $tokens[$i] eq '?' );
98             # $connect->{ $labels[$i] } = _unhack_quote($tokens[$i]);
99             # }
100             # $info->{connect}{ $tokens[0] } = $connect;
101              
102             # }
103             # }
104             if (/^_entity_poly.entity_id\s*(\d+)?/) {
105             chomp;
106 0 0       0 # suck up the entity sequences, do it until # if in a loop
107 0         0 my $entity_id = $1;
108              
109 0         0 if ($entity_id) {
110             die "should not be in loop..." if $in_loop;
111 0 0       0 my $seq;
112 0 0       0 my $seq_can;
113 0         0 while ( my $line = <$fh> ) {
114             chomp($line);
115 0         0 if ( $line =~ /^_entity_poly\.pdbx_seq_one_letter_code\s+(\S+)?\s*$/ )
116 0         0 {
117 0 0       0 if ($1){
118             $seq = $1;
119 0 0       0 }
120 0         0 else {
121             while (my $next_line = <$fh>){
122             chomp $next_line;
123 0         0 last if $next_line =~ /^;\s*$/;
124 0         0 $seq .= $next_line;
125 0 0       0 }
126 0         0 }
127             }
128              
129             if ( $line =~ /^_entity_poly\.pdbx_seq_one_letter_code_can\s+(\S+)?\s*$/ )
130             {
131 0 0       0 if ($1){
132             $seq_can = $1;
133 0 0       0 }
134 0         0 else {
135             while (my $next_line = <$fh>){
136             chomp $next_line;
137 0         0 last if $next_line =~ /^;\s*$/;
138 0         0 $seq_can .= $next_line;
139 0 0       0 }
140 0         0 }
141             }
142            
143             last if $line =~ /^#/;
144            
145 0 0       0 }
146             $seq =~ s/[;'"]//g;
147             $seq_can =~ s/[;'"]//g;
148 0         0 $info->{entity}{$entity_id}{'_entity_poly.pdbx_seq_one_letter_code'} = $seq;
149 0         0 $info->{entity}{$entity_id}{'_entity_poly.pdbx_seq_one_letter_code_can'} = $seq_can;
150 0         0 }
151 0         0 else {
152             die "expected loop exception" unless ($in_loop);
153             # we are in a loop
154 0 0       0 my @labels = ($_);
155             #my @labels = ("_entity_poly.entity_id");
156 0         0 # TODO FACTOR OUT DRY
157             while ( my $line = <$fh> ) {
158             if ( $line =~ /^(_entity_poly.\S+)/ ) {
159 0         0 push @labels, $1;
160 0 0       0 next;
161 0         0 }
162 0         0 if ( $line =~ /^\#/ ) {
163             $in_loop = 0;
164 0 0       0 last;
165 0         0 }
166 0         0 chomp($line);
167             $line = _quote_hack($line);
168 0         0 my @tokens = split /\s+/, $line;
169 0         0 # suck up lines until @tokens == @lables
170 0         0 while ( scalar(@tokens) < scalar(@labels) ) {
171             my $next_line = <$fh>;
172 0         0 $next_line = _quote_hack ($next_line);
173 0         0 chomp $next_line;
174 0         0 my $seq;
175 0         0 if( $next_line =~ /^;/ ){
176 0         0 $seq .= $next_line;
177 0 0       0 while (my $seq_line = <$fh>){
178 0         0 chomp($seq_line);
179 0         0 #$seq_line = _quote_hack($seq_line);
180 0         0 last if $seq_line =~ /;\s*$/;
181             $seq .= $seq_line;
182 0 0       0 }
183 0         0 $tokens[$#tokens+1] = $seq;
184             }
185 0         0 else {
186              
187             push @tokens, split /\s+/, $next_line;
188             }
189 0         0 }
190            
191             die "too many tokens! [@labels] [@tokens]" if scalar(@tokens) != scalar(@labels);
192              
193 0 0       0 my $entity;
194             foreach my $i ( 1 .. $#labels ) {
195 0         0 next if ( $tokens[$i] eq '.' || $tokens[$i] eq '?' );
196 0         0 ($entity->{ $labels[$i] } = _unhack_quote($tokens[$i])) =~ s/[;'"]//g;
197 0 0 0     0 }
198 0         0 $info->{entity}{ $tokens[0] } = $entity;
199             }
200 0         0  
201             # $pdbx_seq_one_letter_code = 0;
202             # while ( my $line = <$fh> ) {
203             # chomp($line);
204             # if ( $line =~ /^(\d+)\s/ ) {
205             # $pdbx_seq_one_letter_code = 1;
206             # $entity_id = $1;
207             # next;
208             # }
209             # if ($pdbx_seq_one_letter_code) {
210             # if ( $line =~ /^;$/ )
211             # { # taking the first sequence that ends with ^;\n
212             # $pdbx_seq_one_letter_code = 0;
213             # next;
214             # }
215             # $seq = $line;
216             # $seq =~ s/(\s|;)//g;
217             # $info->{entity}{$entity_id} .= $seq;
218             # }
219             # if ( $line =~ /^\#/ ) {
220             # $in_loop = 0;
221             # last;
222             # }
223             # }
224             }
225             }
226             if (/^_struct_keywords.text\s*(?:\'?(.+)\'?)?\s*$/) {
227             if ($1){
228 0 0       0 $info->{keywords} = $1;
229 0 0       0 }
230 0         0 else {
231             while (my $line = <$fh>){
232             chomp($line);
233 0         0 last if $line =~ /^#/;
234 0         0 $info->{keywords} .= $line
235 0 0       0 }
236 0         0 }
237             $info->{keywords} =~ s/^;|;$|['"]|\s+$//g;
238             }
239 0         0 if (/^_exptl\.(\S+)\s*(?:\'(.+)\')?/){ # 2ah8 changes order
240             if ($1 eq 'method'){
241 0 0       0 $info->{exp_method} = $2;
242 0 0       0 }
    0          
243 0         0 elsif ($in_loop){
244             my @labels = ($_);
245             my $exptl = {};
246 0         0 while (my $line = <$fh>){
247 0         0 if ( $line =~ /^(_exptl\S+)/ ) {
248 0         0 push @labels, $1;
249 0 0       0 next;
250 0         0 }
251 0         0  
252             if ( $line =~ /^\#/ ) {
253             $in_loop = 0;
254 0 0       0 last;
255 0         0 }
256 0         0 chomp($line);
257             $line = _quote_hack($line);
258 0         0 my @tokens = split /\s+/, $line;
259 0         0 die "too many tokens! [@labels] [@tokens]" if scalar(@tokens) != scalar(@labels);
260 0         0  
261 0 0       0 foreach my $i (0 .. $#labels){
262             push @{$exptl->{$labels[$i]}},$tokens[$i]
263 0         0 }
264 0         0 }
  0         0  
265            
266             $info->{exp_method} = join ';', map {
267             my $str = _unhack_quote($_);
268             $str =~ s/['"]//g; $str } @{$exptl->{'_exptl.method'}};
269 0         0 }
270 0         0 else {
  0         0  
  0         0  
  0         0  
271             while (my $line = <$fh>){
272             last if $line =~ /^#/;
273 0         0 if ($line =~ /^_exptl\.method\s+(?:\'(.+)\')?/) {
274 0 0       0 die "unable to parse exp_method" unless($1);
275 0 0       0 $info->{exp_method} = $1;
276 0 0       0 last;
277 0         0 }
278 0         0
279             }
280             }
281             }
282             if (/^_refine_hist\.d_res_high\s+(\d+\.\d+)/) {
283             $info->{resolution} = sprintf("%s",$1);
284 0 0       0 }
285 0         0 if (/^_em_3d_reconstruction.resolution\s+(\d+\.\d+)/) {
286             $info->{resolution} = $1;
287 0 0       0 }
288 0         0 if (/^_citation.pdbx_database_id_DOI\s+(\S+)/) {
289             $info->{doi} = $1;
290 0 0       0 }
291 0         0 if (/^_atom_site.group_PDB/){
292             $atom_site_position = $. - 2 ;
293 0 0       0 $info->{fh_position_atom_site} = $atom_site_position ;
294 0         0 }
295 0         0 }
296            
297             # move $fh to atom_site_position
298             seek ($fh, 0, 0 );
299             $. = 0;
300 0         0 my $LINE;
301 0         0 do { $LINE = <$fh> } until $. == $atom_site_position || eof;
302 0         0  
303 0   0     0 return $info;
  0         0  
304             }
305 0         0  
306             # watch out for things like 'C-U MISPAIR' in 3j7p
307             my $line = shift;
308             my @q_strs = $line =~ /['"](?:[^"'\\]++|\\.)*+["']/g; # via man perlre
309             foreach my $q_str (@q_strs){
310 0     0   0 (my $hack_str = $q_str) =~ s/(\s+)/__HACK__/g;
311 0         0 $line =~ s/$q_str/$hack_str/;
312 0         0 }
313 0         0 return $line;
314 0         0 }
315              
316 0         0 # watch out for things like 'C-U MISPAIR' in 3j7p
317             my $line = shift;
318             $line =~ s/__HACK__/ /g ;
319             return $line;
320             }
321 0     0   0  
322 0         0  
323 0         0 # merge atoms here???
324             my $self = shift;
325             my ( $atoms, $models, $count_check ) = $self->_read_cif_atoms(@_);
326             my @sets;
327             foreach my $model_num ( @{$models} ) {
328             push @sets, [ grep { $_->model_num == $model_num } @$atoms ];
329 1     1 0 2 }
330 1         6 wantarray ? return @sets : return $sets[0];
331 1         3 }
332 1         2  
  1         5  
333 2         4 my %line_parser_dp = (
  124         2367  
334             'def' => sub {
335 1 50       30 my (
336             $record_name, $serial, $element, $name,
337             $altloc, $res_name, $chain_id, $entity_id,
338             $res_id, $icod, $x, $y,
339             $z, $occ, $B, $charge,
340             $auth_seq_id, $auth_comp_id, $auth_asym_id, $auth_atom_id,
341             $model_num
342             ) = split;
343              
344             $element = ucfirst( lc($element) );
345             my $xyz = V( $x, $y, $z );
346              
347             my $atom = HackaMol::Atom->new(
348             name => $name,
349             record_name => $record_name,
350             serial => $serial,
351             chain => $chain_id,
352             symbol => $element,
353             coords => [$xyz],
354             occ => $occ * 1,
355             bfact => $B * 1,
356             resname => $res_name,
357             resid => $res_id,
358             auth_seq_id => $auth_seq_id,
359             auth_comp_id => $auth_comp_id,
360             auth_asym_id => $auth_asym_id,
361             auth_atom_id => $auth_atom_id,
362             entity_id => $entity_id,
363             model_num => $model_num,
364             );
365             $atom->altloc($altloc) if ( $altloc ne '.' );
366             $atom->icode($icod) if ( $icod ne '?' );
367             return ($atom, $model_num);
368             },
369             'esd' => sub {
370             my (
371             $record_name, $serial, $element, $name,
372             $altloc, $res_name, $chain_id, $entity_id,
373             $res_id, $icod, $x, $y,
374             $z, $occ, $B,
375             undef, undef, undef, undef, undef, # ignor esd for now
376             $charge,
377             $auth_seq_id, $auth_comp_id, $auth_asym_id, $auth_atom_id,
378             $model_num
379             ) = split;
380              
381             $element = ucfirst( lc($element) );
382             my $xyz = V( $x, $y, $z );
383              
384             my $atom = HackaMol::Atom->new(
385             name => $name,
386             record_name => $record_name,
387             serial => $serial,
388             chain => $chain_id,
389             symbol => $element,
390             coords => [$xyz],
391             occ => $occ * 1,
392             bfact => $B * 1,
393             resname => $res_name,
394             resid => $res_id,
395             auth_seq_id => $auth_seq_id,
396             auth_comp_id => $auth_comp_id,
397             auth_asym_id => $auth_asym_id,
398             auth_atom_id => $auth_atom_id,
399             entity_id => $entity_id,
400             model_num => $model_num,
401             );
402             $atom->altloc($altloc) if ( $altloc ne '.' );
403             $atom->icode($icod) if ( $icod ne '?' );
404             return ($atom, $model_num);
405             },
406              
407             );
408              
409              
410             # read pdb file and generate list of Atom objects
411             # no model branching
412             my $self = shift;
413             my $fh = shift;
414             my @atoms;
415             my %track_models;
416              
417             my @expected_labels = (
418 1     1   2 "_atom_site.group_PDB",# (record_name) _atom_site.group_PDB
419 1         2 "_atom_site.id",# (serial) _atom_site.id
420 1         2 "_atom_site.type_symbol",# (element) _atom_site.type_symbol
421             "_atom_site.label_atom_id",# (name) _atom_site.label_atom_id
422             "_atom_site.label_alt_id",# (altloc) _atom_site.label_alt_id
423 1         8 "_atom_site.label_comp_id",# (res_name) _atom_site.label_comp_id
424             "_atom_site.label_asym_id",# (chain_id) _atom_site.label_asym_id
425             "_atom_site.label_entity_id",# (entity_id) _atom_site.label_entity_id
426             "_atom_site.label_seq_id",# (res_id) _atom_site.label_seq_id
427             "_atom_site.pdbx_PDB_ins_code",# (icod) _atom_site.pdbx_PDB_ins_code
428             "_atom_site.Cartn_x",# (x) _atom_site.Cartn_x
429             "_atom_site.Cartn_y",# (y) _atom_site.Cartn_y
430             "_atom_site.Cartn_z",# (z) _atom_site.Cartn_z
431             "_atom_site.occupancy",# (occ) _atom_site.occupancy
432             "_atom_site.B_iso_or_equiv",# (B) _atom_site.B_iso_or_equiv
433             "_atom_site.pdbx_formal_charge",# (charge) #_atom_site.Cartn_x_esd #
434             "_atom_site.auth_seq_id",# (auth_seq_id) #_atom_site.Cartn_y_esd #
435             "_atom_site.auth_comp_id",# (auth_comp_id) #_atom_site.Cartn_z_esd #
436             "_atom_site.auth_asym_id",# (auth_asym_id) #_atom_site.occupancy_esd #
437             "_atom_site.auth_atom_id",# (auth_atom_id) #_atom_site.B_iso_or_equiv_esd #
438             "_atom_site.pdbx_PDB_model_num",# (model_num) #_atom_site.pdbx_formal_charge
439             #_atom_site.auth_seq_id
440             ); #_atom_site.auth_comp_id
441             #_atom_site.auth_asym_id
442             #_atom_site.auth_atom_id
443             #_atom_site.pdbx_PDB_model_num
444            
445             # look for the labels to determine order of attrs
446             my $site_type = 'def'; # 'esd'
447             LABEL: while (<$fh>){
448             if (/^_atom_site.Cartn_x_esd/){
449             $site_type = 'esd';
450             last;
451             }
452 1         3 last if (/^_atom_site.pdbx_PDB_model_num/);
453 1         5  
454 9 50       19 if (/^_atom_site.group_PDB/){
455 0         0 my $i = 1;
456 0         0 while (my $line = <$fh>){
457             die "labels out of order $. $line $expected_labels[$i]" unless $line =~ /$expected_labels[$i]\s?/;
458 9 100       16 last if ($line =~ /^_atom_site.B_iso_or_equiv/);
459             $i++;
460 8 100       22 }
461 1         2 }
462 1         4 }
463 14 50       156  
464 14 100       33 while (<$fh>) {
465 13         29 if ( $self->has_readline_func ) {
466             next if $self->readline_func->($_) eq 'PDB_SKIP';
467             }
468              
469             last if (/^#/); # this works in combination with order search above
470 1         8  
471 62 50       1620 if (/^(?:HETATM|ATOM)\s/) {
472 0 0       0 my ($atom, $model_num) = $line_parser_dp{$site_type}->($_);
473             $track_models{$model_num}++;
474             push @atoms, $atom;
475 62 50       124 }
476             }
477 62 50       168  
478 62         130 my %model_atom_count = reverse %track_models;
479 62         93  
480 62         299 if ( keys(%model_atom_count) > 1 ) {
481             warn "atom number inconsistencies in models. use third return value\n";
482             }
483              
484 1         7 return ( \@atoms, [ sort { $a <=> $b } keys %track_models ],
485             \%model_atom_count );
486 1 50       6 }
487 0         0  
488             no Moose::Role;
489              
490 1         9 1;
  1         15  
491              
492              
493             =pod
494 11     11   137  
  11         29  
  11         81  
495             =head1 NAME
496              
497             HackaMol::Roles::ReadPdbxRole - parse PDBx/mmCIF files
498              
499             =head1 VERSION
500              
501             version 0.052
502              
503             =head1 SYNOPSIS
504              
505             my @atoms = HackaMol->new
506             ->read_cif_atoms("some.cif");
507              
508             =head1 DESCRIPTION
509              
510             The HackaMol::Roles::ReadPdbxRole provides methods (such as read_cif_atoms for building HM atoms) for pulling information from PDBx/mmcif.
511             More info available in cif and more attributes are populated than via ReadPdbRole.
512              
513             =head1 METHODS
514              
515             =head2 read_pdb_atoms
516              
517             One argument: the filehandle
518             Returns a list of HackaMol::Atom objects.
519              
520             The implementation differs from the PDB parser in how models are handled. Each atom tracks the model number through the model_num attribute;
521             the read_pdb_atoms adds configurations to each subsequent atom into t (with identical metadata) as model number increases. Still considering the
522             best way to handle this.
523              
524             =head1 SEE ALSO
525              
526             =over 4
527              
528             =item *
529              
530             L<HackaMol>
531              
532             =item *
533              
534             L<HackaMol::Atom>
535              
536             =item *
537              
538             L<HackaMol::Roles::MolReadRole>
539              
540             =item *
541              
542             L<Protein Data Bank|http://pdb.org>
543              
544             =back
545              
546             =head1 AUTHOR
547              
548             Demian Riccardi <demianriccardi@gmail.com>
549              
550             =head1 COPYRIGHT AND LICENSE
551              
552             This software is copyright (c) 2017 by Demian Riccardi.
553              
554             This is free software; you can redistribute it and/or modify it under
555             the same terms as the Perl 5 programming language system itself.
556              
557             =cut