Fixing issue in getAlignedResIndex, now the mapping will work for

josemduarte · josemduarte · commit 7f1a7f3cb8cf · 2015-08-03T16:09:54.000+02:00
residues that are observed in a chain but unobserved in others
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/Compound.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/Compound.java
@@ -490,14 +490,39 @@ public int getAlignedResIndex(Group g, Chain c) {
 		Map<ResidueNumber,Integer> map = chains2pdbResNums2ResSerials.get(c.getChainID());
 		int serial;
 		if (map!=null) {
-			Integer alignedSerial = map.get(g.getResidueNumber());
-
-			if (alignedSerial==null) {
-				// the map doesn't contain this group, something's wrong: return -1
+			
+			ResidueNumber resNum = g.getResidueNumber();
+			// the resNum will be null for groups that are SEQRES only and not in ATOM, 
+			// still it can happen that a group is in ATOM in one chain but not in other of the same compound.
+			// This is what we try to find out here (analogously to what we do in initResSerialsMap() ):			
+			if (resNum==null && c.getSeqResGroups()!=null && !c.getSeqResGroups().isEmpty()) {
+				int index = -1;
+				for (int i=0;i<c.getSeqResGroups().size();i++) {
+					if (g==c.getSeqResGroup(i)) {
+						index = i; break; 
+					}
+				}
+				
+				resNum = findResNumInOtherChains(index, c);
+								
+			}  
+
+			if (resNum == null) { 
+				// still null, we really can't map
 				serial = -1;
-			} else {
-				serial = alignedSerial;
 			}
+			else {
+
+				Integer alignedSerial = map.get(resNum);
+
+				if (alignedSerial==null) {
+					// the map doesn't contain this group, something's wrong: return -1
+					serial = -1;
+				} else {
+					serial = alignedSerial;
+				}
+			}
+			
 		} else {
 			// no seqres groups available we resort to using the pdb residue numbers are given
 			serial = g.getResidueNumber().getSeqNum();
@@ -519,9 +544,44 @@ private void initResSerialsMap(Chain c) {
 		chains2pdbResNums2ResSerials.put(c.getChainID(), resNums2ResSerials);
 		
 		for (int i=0;i<c.getSeqResGroups().size();i++) {
-			resNums2ResSerials.put(c.getSeqResGroup(i).getResidueNumber(), i+1);
+			
+			// The seqres group will have a null residue number whenever its corresponding atom group doesn't exist
+			// because it is missing in the electron density.
+			// However, it can be observed in the density in other chains of the same compound,
+			// to be complete we go and look for the residue number in other chains, so that we have a 
+			// seqres to atom mapping as complete as possible (with all known atom groups of any chain of this compound)
+			
+			ResidueNumber resNum = c.getSeqResGroup(i).getResidueNumber();
+			
+			if (resNum==null) {
+				resNum = findResNumInOtherChains(i,c);
+			}
+			
+			// NOTE that resNum will still be null here for cases where the residue 
+			// is missing in atom groups (not observed in density) in all chains
+			// Thus the mapping will not be possible for residues that are only in SEQRES groups
+			resNums2ResSerials.put(resNum, i+1);
 		}
 	}
+	
+	private ResidueNumber findResNumInOtherChains(int i, Chain chain) {
+		for (Chain c: getChains()) {
+			if (c == chain) continue;
+			
+			Group seqResGroup = c.getSeqResGroup(i);
+			
+			if (seqResGroup==null) {
+				logger.warn("The SEQRES group is null for index {} in chain {}, whilst it wasn't null in chain {}",
+						 i, c.getChainID(), chain.getChainID());
+				continue;
+			}
+			
+			if (seqResGroup.getResidueNumber()!=null) return seqResGroup.getResidueNumber();
+			
+		}
+		
+		return null;
+	}
 
 	/**
 	 * Return the ref chain id value.
