view release on metacpan or  search on metacpan

t/cif.t  view on Meta::CPAN

#!/usr/bin/env perl
# mmCIF.
#
# The claim this file exists to check is a single one: which format a
# structure arrived in does not change the structure.  So most of what is
# below reads the same structure twice, once from a .pdb and once from a .cif,
# and asserts that the two are equal -- not similar, not equal in the fields
# that were thought of, but equal, by is_deeply, over the whole coordinate
# half of the returned hash.
#
# The fixture pairs are written by t/data/generate.pl, which converts the PDB
# records into the mmCIF loop rather than typing the atoms twice, so that any
# difference between the two files is one the generator put there on purpose.
# The one it puts there on purpose is the naming: the .cif files carry
# label_asym_id and label_seq_id that deliberately disagree with the chain ids
# and residue numbers, because auth_* is what a PDB record carries and auth_*
# is what a reader has to use.
require 5.010;
use strict;
use warnings FATAL => 'all';
use Cwd 'abs_path';
use File::Basename 'dirname';
use File::Temp 'tempdir';
use Chem::Structure::Parser;
use Test::Exception;
use Test::More;

my $data = dirname(abs_path(__FILE__)) . '/data';

# coords() -- the half of a parsed structure that the coordinates decide: every
# chain, every residue, every atom, and the counts over them.
#
# What is taken back out is the handful of things _chain_stats() folds into a
# chain from the header, because those are the header's answer and not the
# coordinates', and the header section below tests them against each other one
# at a time.  A fixture pair cannot be equal on all of them anyway -- mini.cif
# names the entity every ligand belongs to and mini.pdb has no record that
# does -- and burying that in this comparison would only make it say
# 'structures differ' about something the header tests say properly.
my @FROM_HEADER = qw(seqres seqres_length n_missing mol_id molecule organism
                     fragment ec dbref);

sub coords {
	my ($i) = @_;
	my %s = %{ $i->{stats} };
	delete $s{n_lines};    # an mmCIF row and a PDB record are not the same line
	my $strip = sub {
		my ($set) = @_;
		return { map {
			my %c = %{ $set->{$_} };
			delete @c{@FROM_HEADER};
			$_ => \%c;
		} keys %$set };
	};
	return {
		chains      => $strip->($i->{chains}),
		chain_order => $i->{chain_order},
		stats       => \%s,
		model       => $i->{model},
		n_models    => $i->{n_models},
		(exists $i->{models}
			? (models => { map { $_ => { chains      => $strip->($i->{models}{$_}{chains}),
			                             chain_order => $i->{models}{$_}{chain_order} } }
			               keys %{ $i->{models} } })
			: ()),
	};
}

#--------------------------------------------------------------------
# the same structure, both ways
#--------------------------------------------------------------------
for my $pair ([ 'mini', 'one of everything' ],
              [ 'nmr',  'a three model ensemble' ],
              [ 'bare', 'coordinates and nothing else' ]) {
	my ($stem, $what) = @$pair;
	my $p = structure_info("$data/$stem.pdb");
	my $c = structure_info("$data/$stem.cif");

	is($p->{format}, 'pdb',   "$stem.pdb is read as PDB");
	is($c->{format}, 'mmcif', "$stem.cif is read as mmCIF");
	is_deeply(coords($c), coords($p),
		"$stem: $what -- every chain, residue, atom and count is the same from either file");

	# and the views over it, which are what most callers actually touch
	is_deeply(structure_atoms($c),    structure_atoms($p),    "$stem: structure_atoms agrees");
	is_deeply(structure_residues($c), structure_residues($p), "$stem: structure_residues agrees");
	is_deeply(structure_ligands($c),  structure_ligands($p),  "$stem: structure_ligands agrees");
	is_deeply(structure_sequences($c), structure_sequences($p), "$stem: structure_sequences agrees");
	is_deeply(structure_sequences("$data/$stem.cif"), structure_sequences("$data/$stem.pdb"),
		"$stem: structure_sequences agrees when handed the file name");
}

# bare.cif has no _atom_site.type_symbol, as bare.pdb has no element columns,
# so both readers have to get the element out of the atom name -- and get the
# same answer, which is what stops a CA becoming calcium in one and carbon in
# the other
{
	my $c = structure_info("$data/bare.cif");
	is($c->{chains}{A}{residues}{1}{atoms}{CA}{element}, 'C',
		'with no type_symbol the element is worked out from the atom name');
	is_deeply($c->{stats}{elements}, structure_info("$data/bare.pdb")->{stats}{elements},