view release on metacpan or  search on metacpan

t/real_cif.t  view on Meta::CPAN

# Windows has no HOME, and a smoker runs with warnings fatal, so reading it
# to build the maintainer's default paths is what failed the file there rather
# than skipping it.  USERPROFILE is the Windows spelling; with neither set
# there is no default path to try, and the entry drops out.
my $home = defined $ENV{HOME} ? $ENV{HOME} : $ENV{USERPROFILE};

# canon_charges() -- both readers' charges in one spelling, in place.
#
# The two formats write a formal charge differently and each reader reports
# what its own format wrote: a PDB file saying "1-" reads back as '1-', and the
# mmCIF rendering of that same charge, -1, reads back as '1-' too.  That is the
# conversion working.  But the archive holds PDB files that write the mmCIF
# spelling in the PDB column -- 4byf has "-1" in columns 79-80 -- and there the
# PDB reader passes through what it found while the mmCIF reader normalises, so
# the same charge comes back spelled two ways.
#
# Both are right.  So the comparison asks whether it is the same charge rather
# than whether it is the same string, and a wrong charge -- 1- against 2-, or
# 1- against 1+ -- still fails.  In place because $info is dropped straight
# after, and copying every atom to rewrite one field of it would cost more than
# the comparison.
sub canon_charges {
	my ($i) = @_;
	for my $cid (@{ $i->{chain_order} }) {
		my $c = $i->{chains}{$cid};
		for my $rk (@{ $c->{residue_order} }) {
			my $r = $c->{residues}{$rk};
			for my $an (@{ $r->{atom_order} }) {
				my $a = $r->{atoms}{$an};
				my $q = $a->{charge};
				next unless defined $q && length $q;
				# sign and magnitude, whichever order the field put them in.
				# Zero has no sign: 4iu3 writes "0-" in columns 79-80, and a
				# charge of minus nothing is a charge of nothing.
				my ($sign, $mag) =
					  $q =~ /\A(\d)([-+])\z/ ? ($2, $1)
					: $q =~ /\A([-+])(\d)\z/ ? ($1, $2)
					: $q =~ /\A(\d)\z/       ? ('+', $1)
					:                          (undef, undef);
				next unless defined $mag;              # not a charge; leave it be
				$a->{charge} = $mag == 0 ? '0' : "$sign$mag";
			}
		}
	}
	return $i;
}

# the coordinate half, which is what has to survive the change of format
sub coords {
	my ($i) = @_;
	canon_charges($i);
	my %s = %{ $i->{stats} };
	delete $s{n_lines};    # an mmCIF row and a PDB record are not the same line
	my %c;
	for my $cid (keys %{ $i->{chains} }) {
		my %x = %{ $i->{chains}{$cid} };
		# what _chain_stats() folded in from the header, which the conversion
		# below does not carry across and which t/cif.t tests on its own
		delete @x{qw(seqres seqres_length n_missing mol_id molecule organism
		             fragment ec dbref)};
		$c{$cid} = \%x;
	}
	return { chains => \%c, chain_order => $i->{chain_order}, stats => \%s };
}

# --- consistency, for a structure with nothing to compare it against -------
sub adds_up {
	my ($info, $name) = @_;
	my $atoms = 0;
	$atoms += $info->{chains}{$_}{n_atoms} for @{ $info->{chain_order} };
	return "the chains hold $atoms atoms but the file had $info->{stats}{n_atoms}"
		if $atoms != $info->{stats}{n_atoms};
	for my $cid (@{ $info->{chain_order} }) {
		my $c = $info->{chains}{$cid};
		my ($n, $poly) = (0, 0);
		for my $rk (@{ $c->{residue_order} }) {
			my $r = $c->{residues}{$rk};
			$n += $r->{n_atoms};
			$poly++ if $r->{type} eq 'amino_acid' || $r->{type} eq 'nucleotide';
			return "chain $cid residue $rk has no atoms" unless $r->{n_atoms};
			return "chain $cid residue $rk is keyed wrong" unless $r->{key} eq $rk;
		}
		return "chain $cid: residues hold $n atoms, the chain says $c->{n_atoms}"
			if $n != $c->{n_atoms};
		return "chain $cid: $poly polymer residues but a sequence of "
		       . length($c->{sequence})
			if length($c->{sequence}) != $poly;
		return "chain $cid: $c->{n_residues} residues but "
		       . scalar(@{ $c->{residue_order} }) . ' in the order'
			if $c->{n_residues} != @{ $c->{residue_order} };
	}

	# and the header, for the files that have one.  The .cif files a
	# simulation writes carry no header at all, so this says nothing about
	# them; point STRUCTURE_INFO_TEST_CIF_DIR at a directory of archive
	# entries and it is the half of the reader that gets exercised.
	for my $k (qw(resolution r_work r_free temperature ph)) {
		next unless defined $info->{$k};
		return "$k is '$info->{$k}', which is not a number"
			unless $info->{$k} =~ /\A-?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?\z/;
	}
	for my $k (qw(keywords experiment authors helix sheet ssbond link cispep)) {
		return "$k is not an array reference" unless ref $info->{$k} eq 'ARRAY';
	}
	for my $k (qw(header seqres het remarks modres dbref cryst1 journal)) {
		return "$k is not a hash reference" unless ref $info->{$k} eq 'HASH';
	}
	for my $cid (keys %{ $info->{seqres} }) {
		my $s = $info->{seqres}{$cid};
		return "seqres $cid has no sequence" unless defined $s->{sequence};
		return "seqres $cid: length $s->{length} but a sequence of "
		       . length($s->{sequence})
			if defined $s->{length} && @{ $s->{residues} }
			&& $s->{length} != length $s->{sequence};
	}
	return undef;
}

#---------------------------------------------------------------------------
# 1. real .cif files