Chem-Structure-Parser
view release on metacpan or search on metacpan
lib/Chem/Structure/Parser.pm view on Meta::CPAN
}
return $c unless defined $lo;
for my $rk (@{ $c->{residue_order} }) {
my $r = $c->{residues}{$rk};
next unless $r->{hetero};
next unless $r->{type} eq 'amino_acid' || $r->{type} eq 'nucleotide';
next unless defined $r->{number};
# one either side, so that a modified residue capping a terminus is
# still part of the chain
next if $r->{number} >= $lo - 1 && $r->{number} <= $hi + 1;
$r->{type} = 'ligand';
$r->{one} = '';
$r->{modified} = 0;
$r->{free} = 1; # a free amino acid, not part of the polymer
}
return $c;
}
sub _chain_type {
my ($count, $poly) = @_;
my $aa = $count->{amino_acid} || 0;
my $nuc = $count->{nucleotide} || 0;
if ($aa || $nuc) {
return 'protein' if $aa >= $nuc;
my $deoxy = grep { $_->{resname} =~ /\AD[ACGTUI]\z/ } @$poly;
return $deoxy * 2 >= $nuc ? 'dna' : 'rna';
}
return 'water' if ($count->{water} || 0) && !($count->{ligand} || 0) && !($count->{ion} || 0);
return 'hetero' if ($count->{ligand} || 0) || ($count->{ion} || 0) || ($count->{water} || 0);
return 'unknown';
}
# SEQRES, COMPND and SOURCE all describe chains; fold them in once the chains
# exist, so that everything about a chain is in one place
sub _chain_stats {
my ($info) = @_;
# The free-text COMPND of an old file names no chains, so there was no
# chain to file it under when the header was read and there is one now: the
# entry is the one molecule and every chain in it is that molecule. This is
# the only place a chain is added to entity_of_chain, because it is the only
# entity that could not say for itself which chains it means.
if (($info->{compound}{1} || {})->{free_text} && !%{ $info->{entity_of_chain} }) {
my $s = $info->{source}{1} || {};
$info->{entity_of_chain}{$_} = {
mol_id => 1,
molecule => $info->{compound}{1}{molecule},
organism => $s->{organism_scientific},
} for @{ $info->{chain_order} };
}
for my $cid (@{ $info->{chain_order} }) {
my $c = $info->{chains}{$cid};
if (my $s = $info->{seqres}{$cid}) {
$c->{seqres} = $s->{sequence};
$c->{seqres_length} = $s->{length};
$c->{n_missing} = $s->{length} - $c->{n_polymer} if defined $s->{length};
}
if (my $e = $info->{entity_of_chain}{$cid}) {
$c->{mol_id} = $e->{mol_id};
$c->{molecule} = $e->{molecule} if defined $e->{molecule};
$c->{organism} = $e->{organism} if defined $e->{organism};
$c->{fragment} = $e->{fragment} if defined $e->{fragment};
$c->{ec} = $e->{ec} if defined $e->{ec};
}
$c->{dbref} = $info->{dbref}{$cid} if $info->{dbref}{$cid};
}
return $info;
}
sub _id_from {
my ($info, $file) = @_;
return $info->{header}{id_code} if length($info->{header}{id_code} || '');
return undef unless defined $file;
my ($base) = $file =~ m{([^/\\]+)\z};
$base =~ s/\.(gz|bz2|z)\z//i;
$base =~ s/\.(pdb|ent|cif|mmcif)\z//i;
$base =~ s/\.ent\z//i;
$base =~ s/\Apdb//i;
return uc $base;
}
#
# Header records
#
# Every one of these is a fixed-column record too, but there are only a few
# dozen lines of them in a file, they are irregular, and they are where a new
# quirk turns up every few hundred structures. That is Perl's job, not C's.
#
# The keys a parsed structure always has, whatever was in the file and
# whichever format it was in. Set before either reader runs, so that a caller
# can read $info->{resolution} without first asking whether the file was an
# mmCIF, and get undef for "the file does not say" in both.
sub _meta_defaults {
my ($info) = @_;
$info->{$_} = undef for qw(title resolution r_work r_free);
$info->{$_} = [] for qw(keywords experiment authors);
$info->{$_} = {} for qw(header compound source seqres het hetnam formul
remarks dbref entity_of_chain cryst1 journal
modres);
$info->{$_} = [] for qw(helix sheet ssbond link cispep revdat site conect);
return $info;
}
sub _parse_meta {
my ($info, $meta) = @_;
_meta_defaults($info);
if (my $h = $meta->{HEADER}) {
my $l = $h->[0];
$info->{header} = {
classification => _c($l, 10, 40),
deposit_date => _c($l, 50, 9),
id_code => _c($l, 62, 4),
};
}
# The entry id, for _untail(): the text records of an old file end in it, and
# it is the only thing that tells the stationery from the text.
my $eid = $info->{header}{id_code};
# a record that is not in the file reads as undef, not as an empty string:
# "there was no TITLE" and "the TITLE was blank" are different answers
$info->{title} = $meta->{TITLE} ? _joined($meta->{TITLE}, 10, $eid) : undef;
lib/Chem/Structure/Parser.pm view on Meta::CPAN
}
if (my $n = $meta->{NUMMDL}) {
my $v = _c($n->[0], 10, 4);
$info->{n_models_declared} = $v + 0 if $v =~ /\A\d+\z/;
}
$info->{records} = { map { $_ => scalar @{ $meta->{$_} } } keys %$meta };
return $info;
}
# COMPND and SOURCE are "TOKEN: value;" lists broken into MOL_ID groups
sub _mol_records {
my ($lines, $from, $free_key, $entry_id) = @_;
return {} unless $lines;
my $text = _joined($lines, $from, $entry_id);
my %mol;
my $id = 1;
for my $piece (split /;/, $text) {
next unless $piece =~ /\S/;
my ($k, $v) = $piece =~ /\A\s*([A-Z0-9_ ]+?)\s*:\s*(.*)\z/;
next unless defined $k;
$k = lc $k;
$k =~ s/\s+/_/g;
$v = _t($v);
if ($k eq 'mol_id') {
$id = $v;
$mol{$id}{mol_id} = $v;
next;
}
$mol{$id}{mol_id} = $id unless exists $mol{$id};
if ($k eq 'chain') {
$mol{$id}{chain} = [ grep { length } map { _t($_) } split /,/, $v ];
} else {
$mol{$id}{$k} = exists $mol{$id}{$k} ? "$mol{$id}{$k} $v" : $v;
}
}
# A file older than the MOL_ID convention writes the record as free text --
# 'COMPND GAMMA DELTA RESOLVASE', 'SOURCE (ESCHERICHIA COLI)' -- and a
# reader that knows only about 'MOLECULE:' throws away the one thing the
# record says. There is no chain list in that form because there was
# nothing to distinguish: the whole entry is the one molecule, which is what
# free_text says and what _chain_stats() does with it.
if (!%mol && $free_key && $text =~ /\S/) {
my $v = _t($text);
$v =~ s/\A\((.*)\)\z/$1/; # SOURCE used to parenthesise the organism
%mol = (1 => { mol_id => 1, $free_key => $v, free_text => 1 });
}
return \%mol;
}
# one flat record per chain, so a chain hash can say what molecule it is
sub _entities {
my ($info) = @_;
my %by_chain;
for my $id (keys %{ $info->{compound} }) {
my $c = $info->{compound}{$id};
my $s = $info->{source}{$id} || {};
for my $cid (@{ $c->{chain} || [] }) {
$by_chain{$cid} = {
mol_id => $id,
molecule => $c->{molecule},
fragment => $c->{fragment},
ec => $c->{ec_number} || $c->{ec},
organism => $s->{organism_scientific},
taxid => $s->{organism_taxid},
expressed_in => $s->{expression_system},
};
}
}
$info->{entity_of_chain} = \%by_chain;
return $info;
}
#
# mmCIF header categories
#
# The same facts, filed differently. A PDB file says the resolution on a
# REMARK 2 line and an mmCIF file says it in _refine.ls_d_res_high, and a
# caller who wants to know the resolution should not have to care which. So
# this fills in the same $info keys _parse_meta() fills in, from the
# categories that carry the same information.
#
# Where a fact exists in one format and not the other it is left alone rather
# than invented: an mmCIF file has no REMARK records, so $info->{remarks} stays
# empty, and reading it gets the same "nothing there" a PDB file with no
# remarks would give.
#
# Identifiers are the auth_* ones throughout -- pdbx_strand_id, auth_asym_id --
# because those are the chain ids the coordinates were read under and the ones
# the PDB record carried. Using label_asym_id here would file the annotations
# under chains that the chains hash does not have.
#
sub _parse_cif_meta {
my ($info, $p) = @_;
_meta_defaults($info);
my $cif = $p->{cif} || {};
my $loops = $p->{cif_loops} || {};
my $id = _cif1($p, '_entry', 'id');
$info->{header} = {
classification => _cif1($p, '_struct_keywords', 'pdbx_keywords'),
deposit_date => _cif1($p, '_pdbx_database_status', 'recvd_initial_deposition_date'),
id_code => defined $id ? uc $id : '',
};
# The data_ block name is deliberately not a key of its own. It is usually
# the entry id, which $info->{id} already has, and where it is not -- a file
# written by a simulation program calls its block 'cell' -- it is worse than
# the file name _id_from() falls back to. A key only one of the two formats
# could ever fill in is a key a caller has to test the format for.
$info->{title} = _cif1($p, '_struct', 'title');
$info->{keywords} = [ grep { length } map { _t($_) }
split /,/, (_cif1($p, '_struct_keywords', 'text') || '') ];
$info->{experiment} = [ grep { defined && length }
map { $_->{method} } @{ _cif_rows($p, '_exptl') } ];
$info->{authors} = [ grep { defined && length }
map { $_->{name} } @{ _cif_rows($p, '_audit_author') } ];
# resolution: refined structures say so in _refine; the others say it
# wherever their method says it
for my $where ([ '_refine', 'ls_d_res_high' ],
lib/Chem/Structure/Parser.pm view on Meta::CPAN
h('structure_info'); # by name
h(*res_type); # by name, unquoted
h(\&aa3to1); # by reference
h(); # the list of documented functions
perl -MChem::Structure::Parser -e 'h(*structure_info)' # straight from the shell
Note that C<h(res_type)>, with no quotes and no sigil, cannot be made to work:
every function here is exported, so Perl parses the bareword as a call to
C<res_type()> before C<h> is ever reached. Use one of the three forms above.
=head1 Functions/Subroutines
=head2 structure_info
my $info = structure_info($file, %options);
Reads C<$file> and returns a hash reference. The format is worked out from the
file name — C<.pdb>, C<.ent>, C<.cif>, C<.mmcif>, C<.pdbx> — and from the first
records in the file when the name gives nothing away. C<.gz> files are read as
they are, without unpacking to a temporary file.
=head3 What comes back
Laid out the way C<tree> lays out a directory, this is C<1a22.ent.pdb> — a real
file, real values, the long lists cut short:
$info
├── file '1a22.ent.pdb' the path it was read from
├── format 'pdb' or 'mmcif'
├── id '1A22' from HEADER, or from the file name
├── title 'HUMAN GROWTH HORMONE BOUND TO SINGLE RECEPTOR'
├── header
│ ├── classification 'COMPLEX (HORMONE/RECEPTOR)'
│ ├── deposit_date '15-JAN-98'
│ └── id_code '1A22'
├── experiment [ 'X-RAY DIFFRACTION' ]
├── resolution 2.6 REMARK 2
├── r_work 0.187 REMARK 3
├── r_free undef this entry does not report one
├── temperature 287 REMARK 200
├── ph 6.5
├── keywords [ 'COMPLEX (HORMONE-RECEPTOR)', 'PITUITARY HORMONE', ... ]
├── authors [ 'A.M.DE VOS', 'M.ULTSCH' ]
├── journal
│ ├── auth [ 'T.CLACKSON', 'M.H.ULTSCH', 'J.A.WELLS', 'A.M.DE VOS' ]
│ ├── titl 'STRUCTURAL AND FUNCTIONAL ANALYSIS OF THE 1:1 GROWTH...'
│ ├── ref 'J.MOL.BIOL. V. 277 1111 1998'
│ ├── refn 'ISSN 0022-2836'
│ ├── pmid '9571026'
│ └── doi '10.1006/JMBI.1998.1669'
├── compound COMPND, by MOL_ID
│ ├── 1
│ │ ├── mol_id '1'
│ │ ├── molecule 'GROWTH HORMONE'
│ │ ├── chain [ 'A' ]
│ │ ├── engineered 'YES'
│ │ └── mutation 'YES'
│ └── 2 { molecule 'GROWTH HORMONE RECEPTOR', chain [ 'B' ],
│ fragment 'EXTRACELLULAR DOMAIN', engineered 'YES' }
├── source SOURCE, by MOL_ID
│ └── 1 { organism_scientific 'HOMO SAPIENS', organism_common
│ 'HUMAN', organism_taxid '9606', mol_id '1',
│ expression_system 'ESCHERICHIA COLI',
│ expression_system_taxid '562' }
├── entity_of_chain COMPND and SOURCE, by chain
│ ├── A { mol_id '1', molecule 'GROWTH HORMONE', fragment undef,
│ │ ec undef, organism 'HOMO SAPIENS', taxid '9606',
│ │ expressed_in 'ESCHERICHIA COLI' }
│ └── B { ..., fragment 'EXTRACELLULAR DOMAIN' }
├── seqres what SEQRES says was in the crystal
│ ├── A
│ │ ├── sequence 'FPTIPLSRLFDNAMLRAHRLHQLAFDTYQEFEEAYIPKEQKYSFLQ...'
│ │ ├── residues [ 'PHE', 'PRO', 'THR', 'ILE', ... ] 191 of them
│ │ └── length 191
│ └── B { sequence, residues, length 238 }
├── dbref
│ └── A [ { database 'UNP', accession 'P01241',
│ db_id 'SOMA_HUMAN', seq_begin '1', seq_end '191',
│ db_begin '27', db_end '217', chain 'A' } ]
├── seqadv [ { chain 'A', resseq '120', resname 'ARG',
│ db_res 'GLY', db_seq '146', comment 'ENGINEERED' } ]
├── modres { } no MSE-style residues here
├── het
│ └── HOH { het_id 'HOH', formula '69(H2 O)', water 1 }
├── hetnam { }
├── formul { }
├── helix [ { id '1', class '1', length '29',
│ init_chain 'A', init_resname 'SER', init_resseq '7',
│ end_chain 'A', end_resname 'TYR', end_resseq '35' },
│ ... ] 12 of them
├── sheet [ ... ] 12
├── ssbond [ { chain1 'A', resseq1 '53',
│ chain2 'A', resseq2 '165', length '2.02' }, ... ] 5
├── link [ ]
├── cispep [ ]
├── site [ ]
├── cryst1 { a '67.7', b '67.7', c '228',
│ alpha '90', beta '90', gamma '90',
│ sgroup 'P 43 21 2', z '8' }
├── biological_assembly [ 32 lines of REMARK 350, verbatim ]
├── revdat [ { num '3', date '18-APR-18', id '1A22',
│ type '1', what 'REMARK' }, ... ]
├── remarks every REMARK, by number
│ ├── 2 [ '', 'RESOLUTION. 2.60 ANGSTROMS.' ]
│ ├── 350 [ ... ] 32 lines
│ └── ... 1, 3, 4, 100, 200, 280, 290, 300, 465, 470, 500
├── conect [ [ 448, 1255 ], ... ] 10
├── records every record type, counted
│ ├── REMARK 365
│ ├── SEQRES 34
│ ├── HELIX 12
│ └── ... AUTHOR, COMPND, CONECT, CRYST1, DBREF, SOURCE, SSBOND, ...
├── n_models 1 how many MODEL records the file has
├── model 1 which one the chains below are
├── models there only with model => 'all'
├── stats
│ ├── n_atoms 3113 atoms kept: this model, less what was filtered
│ ├── total_atoms 3113 atoms the file has, every model, unfiltered
│ ├── n_hetatm 69 of n_atoms, the ones written as HETATM
│ ├── n_hydrogens 0
│ ├── n_water_atoms 69
│ ├── n_lines 3605
│ ├── n_atom_records 3044 ATOM lines seen, whether kept or not
│ ├── n_hetatm_records 69 HETATM lines, likewise
│ ├── n_anisou 0
│ ├── n_skipped 0 coordinate lines the options threw away
│ ├── elements { C 1946, O 643, N 507, S 17 }
│ │ every element in the file, keyed by its IUPAC
│ │ symbol; the counts add up to n_atoms
│ ├── bfactor { min '2.7', max '85.39', mean 30.83, n 3113 }
│ ├── bbox { xmin '12.142', xmax '80.34', ymin '2.011', ... }
│ └── center [ '46.241', '29.135', '134.559' ]
├── chain_order [ 'A', 'B' ] the order the file has them in
└── chains
├── A
│ ├── id 'A'
│ ├── type 'protein' protein dna rna water hetero unknown
│ ├── sequence 'FPTIPLSRLFDNAMLRAHRLHQLAFDTYQEFEEAYIPKEQ...'
│ │ single-letter, what has coordinates
│ ├── seqres 'FPTIPLSRLFDNAMLRAHRLHQLAFDTYQEFEEAYIPKEQ...'
│ ├── seqres_length 191
│ ├── n_residues 206
│ ├── n_polymer 180
│ ├── n_water 26
│ ├── n_ligand 0
│ ├── n_atoms 1492
│ ├── n_hetatm 26
│ ├── elements { C 938, O 301, N 246, S 7 }
│ │ the same tally for this chain alone;
│ │ adds up to the chain's n_atoms
│ ├── n_missing 11 SEQRES less what was modelled
│ ├── gaps [ { after 129, before 136, missing 6 },
│ │ { after 148, before 154, missing 5 } ]
│ ├── n_gaps 2
│ ├── missing_residues
│ │ [ 130, 131, 132, 133, 134, 135,
│ │ 149, 150, 151, 152, 153 ]
│ ├── first 1 the first and last polymer residue keys
│ ├── last 191
│ ├── residue_types { amino_acid 180, water 26 }
│ ├── molecule 'GROWTH HORMONE' from COMPND
│ ├── organism 'HOMO SAPIENS' from SOURCE
│ ├── mol_id '1'
│ ├── dbref [ { ... } ] as in the top-level dbref
│ │ ec and fragment are here too, in a
│ │ chain whose file gives them
│ ├── residue_order [ '1', '2', '3', ... '574' ] file order, 206
│ └── residues keyed number + insertion code
│ ├── 54
│ │ ├── resname 'PHE'
│ │ ├── number 54
│ │ ├── icode ''
│ │ ├── key '54'
│ │ ├── chain 'A'
│ │ ├── one 'F' '' when there is no letter for it
│ │ ├── type 'amino_acid'
│ │ │ nucleotide water ligand ion
│ │ ├── standard 1 one of the twenty, or a standard base
│ │ ├── modified 0 1 for MSE, still an M in the sequence
│ │ ├── hetero 0 1 when it was written as HETATM
│ │ ├── free not here; 1 for a free amino acid
│ │ │ bound in a site (see below)
│ │ ├── n_atoms 11
│ │ ├── b_mean 22.55
│ │ ├── center [ 65.311, 17.127, 140.515 ]
│ │ ├── atom_order [ 'N', 'CA', 'C', 'O', 'CB', ... ]
│ │ └── atoms
│ │ ├── CA
│ │ │ ├── name 'CA'
│ │ │ ├── serial 450
│ │ │ ├── element 'C'
│ │ │ ├── charge ''
│ │ │ ├── x '66.446'
│ │ │ ├── y '18.25'
│ │ │ ├── z '141.982'
│ │ │ ├── occupancy '1'
│ │ │ ├── bfactor '24.53'
│ │ │ ├── altloc ''
│ │ │ ├── hetero 0
│ │ │ └── altlocs [ { altloc, x, y, z, occupancy,
│ │ │ bfactor }, ... ]
│ │ │ present only when the atom has
│ │ │ alternate conformers; every conformer
│ │ │ is listed, the chosen one included,
│ │ │ and one of them having no letter at
│ │ │ all does not take it off the list
│ │ └── ... N, C, O, CB, CG, CD1, CD2, CE1, CE2, CZ
│ └── ... 1 .. 191, then the waters at 512 .. 574
└── B the same again: 235 residues, 1621 atoms
A record that is not in the file reads as C<undef>, and a list that is not in
the file reads as an empty arrayref — C<title> being C<undef> means there was no
TITLE, which is a different thing from a TITLE that was blank.
Everything the module does not take apart is still in C<remarks> and in the
raw record counts, so nothing in the file is lost.
=head3 The two sequences
C<sequence> and C<seqres> are the two different questions people mean by "the
sequence": what was modelled, and what was in the crystal. They differ
wherever a terminus or a loop went unmodelled, which is what C<gaps> counts and
C<n_missing> totals — eleven residues of chain A above, in two stretches.
C<missing_residues> is the same eleven one number at a time, in ascending
order, for asking whether a particular residue was modelled without walking
( run in 0.705 second using v1.01-cache-2.11-cpan-364913b4093 )