Hi David,

On Tue, Dec 23, 2008 at 10:31:28AM -0500, Borhani, David wrote:
> Hi Clemens,
> 
> Thanks for all your tests; the scripts/keywords you used to run LSQKAB
> with these test systems would help to clarify what may be going right
> vs. going wrong.

That was just a simple run with

  lsqkab workcd work.pdb refrcd ref.pdb <<EOF
  FIT RESI ALL 1 TO 5
  MATCH 1 TO 5
  EOF

> A few points that I hope may be helpful:
> 
> 1. Atom names are 4 characters. If a true match is desired, all 4
> characters (even the first one that is often a " ") must be compared.

I agree: the 3-character test is not something I'm involved with at
all. This is what LSQKAB is doing - and I _think_ it has to do partly
with MMDB (which does some shifting as far as I can remember).
 
> 2. I think the new code is not correct, as your examples show, and as
> others have found in using the program. The logic, in those cases where
> alt coded atoms are present, seems to be either wrong, unexpected, or
> ill-defined (i.e., it matters which coord set is work vs. reference), or
> perhaps even all of the above.

Maybe my email wasn't quite clear: e.g. the problem with "which coord
set is work vs. reference" happens in the OLD code (CCP3 6.0.2) and is
fixed in the NEW code (6.1). All the unexpected behaviour is present
in the old 6.0.2 version of LSQKAB - the 6.1 version is fixed and
behaves as expected.

The original problem reported was that the output seems to suggest
that LSQKAB was using AltConf atoms in the superposition. And so it
was: but it did it wrongly (i.e. it mattered which PDB file was
defined as reference and which as work). The new code fixes that - but
if you want a feature like

  match up with AltConf "A" in case one PDB has no AltConf and the
  other has

then this needs to be newly introduced into LSQKAB: it never was in
there and the impression that it might have worked like that in the
old version was due to some wrong logic introducing this random
behaviour.

> 3. I agree that chain, residue number, and atom name do not by
> themselves specify a unique atom; insertion code must also be used. Alt
> code *may* be used, and that's where (IMHO) it gets tricky (and
> apparently the older versions of LSQKAB just used an implicit logic
> (i.e., likely matched the first found alt coded atom, without any
> checks):
> 	A. Option one, no defined logic: just ignore the alt code; use
> any (random) atom you get (first).

Yes, that is 6.0.2 behaviour.

> 	B. Option two, defined logic (what that logic should be is the
> key point to discuss, I think).

Yes, that is 6.1 behaviour: exact matching. But yes: there should be
some better message (one can deduce this from the number of atoms in
the working set compared to the number of atoms used - printed just
above the RMS message).

> Rigid potential logic:
> 	1. User must explicitly specify what will constitute a match.
> Absent such a specification, program stops with
> 		an error if alt coded atoms are found (or if they don't
> meet the specification).
> (I don't recommend this!)
> 
> Flexible ("intelligent"?) potential logic:
> 	1. Match the atom without the alt code (i.e., " "), if it
> exists; else match an atom with an altcode.
> 	2. Now matching alt codes:
> 		A. Is there an atom with alt code "A"? Use it, else look
> for "B", "C", etc., in sort order.
> 			(FYI: documentation (6.0.2-03) says: "If there
> are two or more conformations, the first (labelled A) 
> 				is chosen for comparison.")

Ah: that might have been true for pre-mmdb coordinate library use (but
even then it might have depended on the way the PDB file was
written). But with mmdb (as far as I understand) the atoms might be
stored in a different (random?) order.

> 		B. ALTERNATIVELY, use the alt coded atom with the
> highest occupancy; use sort order to 
> 			resolve ties (A > B > C... [usually, one tries,
> at least, to put the most significant 
> 			atom as the " " alt code or the "A" atom]).
> (I prefer using occupancy instead of B factor, because once one is
> modeling alt conformations, occupancy receives some conscious attention;
> the B will then just refine to where it needs to be given the
> user-assigned occupancy [unless one is refining occupancies in SHELX].
> So, in most cases, I suspect, occupancy trumps B factor. Others may
> disagree.)

That would be the nicest way.

> There also may need to be a new keyword/keyvalue to allow the user to
> specify which of several potential alternative logics to use.
> 
> 4. It appears to me that the new version (6.1.0) doesn't have any
> changes to the FIT/MATCH keywords to handle insertion codes. If the user
> specifies, for example:
> 	FIT RESIDU SIDE 155 TO 156 CHAIN A 
> 	MATCH RESIDU 155 TO 156 CHAIN A
> then IF there exists a residue 155A in the working coords, there must
> also be a residue 155A in the ref coords, else error.

Possible: I'm not sure how the sequence "155 TO 156" is interpreted
when using mmdb-based coordinate library. Is 155A stored in "155 155A
156" or are residues with insertion code stored at the end or
... something completely different?

> There are good reasons to allow users to alter this behavior, e.g.
> fitting immunoglobulin hypervariable regions, which often have (a
> variable number of) insertions. Current LSQKAB logic would appear to
> make this task difficult. To be more explicit, if I want to fit residues
> 25-40, knowing that there is a variable loop, with insertion codes after
> residue 30, i.e. I want to fit 25-30 and 31-40, it would be nice to be
> able to specify "25-40" and "SKIP INSERTIONS" or something similar.

I completely agree: insertion- and alternate conformation code
handling in LSQKAB is far from ideal and difficult. But at least
version 6.1 gives the correct answer for the limited scenarios one can
do with LSQKAB.
 
> 5. Finally, ensuring that whatever logic is chosen works no matter which
> coordinate set is specified as work or reference would be highly
> desireable, as your examples clearly point out!

Absolutely - I can't find an official author for LSQKAB, so I guess
those improvements are up for grabs to anyone ;-)

Cheers

Clemens

-- 

***************************************************************
* Clemens Vonrhein, Ph.D.     vonrhein AT GlobalPhasing DOT com
*
*  Global Phasing Ltd.
*  Sheraton House, Castle Park 
*  Cambridge CB3 0AX, UK
*--------------------------------------------------------------
* BUSTER Development Group      (http://www.globalphasing.com)
***************************************************************