use of org.biojava.bio.structure.Atom in project ffx by mjschnie.
the class PDBFileMatcher method simpleRMSD.
private double simpleRMSD(Atom[] atomsA, Atom[] atomsB) throws IllegalArgumentException {
double rmsd = 0.0;
// double comp = 0.0; PDB file precision more important than addition error.
int numMatches = 0;
int numAtomsA = atomsA.length;
for (int i = 0; i < numAtomsA; i++) {
Atom atomFromA = atomsA[i];
try {
Atom atomFromB = getMatchingAtom(atomFromA, atomsB, i);
// Done here in case of exceptions thrown.
++numMatches;
rmsd += getSqDistance(atomFromA, atomFromB);
/*double dist2 = getSqDistance(atomFromA, atomFromB) - comp; // Kahan summation algorithm.
double temp = rmsd + dist2;
comp = (temp - rmsd) - dist2;
rmsd = temp;*/
} catch (IllegalArgumentException ex) {
if (!atomFromA.getElement().isHydrogen()) {
// Will eventually move logger statement here.
}
logger.info(String.format(" Error in finding mate for atom %s", atomFromA.toString()));
}
}
if (numMatches == 0) {
throw new IllegalArgumentException(" No atomic matches found when calculating RMSD.");
}
rmsd = Math.sqrt(rmsd / numMatches);
return rmsd;
}
use of org.biojava.bio.structure.Atom in project ffx by mjschnie.
the class PDBFileMatcher method getMatchingAtom.
/**
* Finds atom1's match in structure2; uses atom1's residue number, atom
* type, and PDB serial number as a guess; if robustMatch is set true, will
* search all Atoms in structure2.
*
* @param atom1 An Atom
* @param structure2 A Structure to search
* @param searchAll Whether to search all Atoms in structure2.
* @return atom1's match in structure2
* @throws IllegalArgumentException If no match can be found.
*/
private Atom getMatchingAtom(Atom atom1, Structure structure2, boolean searchAll) throws IllegalArgumentException {
ResidueNumber res1 = atom1.getGroup().getResidueNumber();
String chainID = res1.getChainId();
Atom atom2 = null;
try {
Chain chain2 = structure2.getChainByPDB(chainID);
Group group2 = chain2.getGroupByPDB(res1);
try {
atom2 = group2.getAtom(atom1.getName());
} catch (StructureException ex) {
if (atom1.getName().equalsIgnoreCase("H")) {
atom2 = group2.getAtom("H1");
} else if (atom1.getName().equalsIgnoreCase("H1")) {
atom2 = group2.getAtom("H");
}
atom2 = group2.getAtom(atom1.getPDBserial());
}
} catch (StructureException ex) {
if (!searchAll) {
throw new IllegalArgumentException("Matching atom not found.");
}
for (Chain chain : structure2.getChains()) {
for (Group group : chain.getAtomGroups()) {
for (Atom atom : group.getAtoms()) {
if (compareAtoms(atom1, atom)) {
return atom2;
}
}
}
}
}
if (atom2 != null && compareAtoms(atom1, atom2)) {
return atom2;
}
throw new IllegalArgumentException("Matching atom not found.");
}
use of org.biojava.bio.structure.Atom in project ffx by mjschnie.
the class PDBFileMatcher method fixFile.
private void fixFile(FileFilePair currentPair, PDBFileReader filereader) throws IOException {
File matchFile = currentPair.getMatchedFile();
Structure matchStructure = currentPair.getStructure();
if (matchStructure == null) {
matchStructure = filereader.getStructure(matchFile);
}
File sourceFile = currentPair.getSourceFile();
if (sourceFile == null) {
throw new IOException(String.format("No source file was matched to file %s", matchFile.toString()));
}
Structure sourceStructure = null;
if (fixAtoms) {
sourceStructure = filereader.getStructure(sourceFile);
Atom[] matchAtoms = StructureTools.getAllAtomArray(matchStructure);
for (Atom matchAtom : matchAtoms) {
Atom sourceAtom = getMatchingAtom(matchAtom, sourceStructure, robustMatch);
if (fixBFactors) {
matchAtom.setTempFactor(sourceAtom.getTempFactor());
}
}
// Other methods can go here.
}
if (fixSSBonds) {
if (sourceStructure == null) {
sourceStructure = filereader.getStructure(sourceFile);
}
List<SSBond> sourceBonds = sourceStructure.getSSBonds();
List<SSBond> matchBonds = matchStructure.getSSBonds();
for (SSBond sourceBond : sourceBonds) {
boolean isContained = false;
for (SSBond matchBond : matchBonds) {
if (compareSSBonds(matchBond, sourceBond)) {
isContained = true;
break;
}
}
if (!isContained) {
matchStructure.addSSBond(sourceBond.clone());
}
}
}
if (fixCryst) {
if (sourceStructure == null) {
sourceStructure = filereader.getStructure(sourceFile);
}
PDBCrystallographicInfo crystalInfo = sourceStructure.getCrystallographicInfo();
try {
PDBCrystallographicInfo duplicateInfo = cloneCrystalInfo(crystalInfo);
matchStructure.setCrystallographicInfo(duplicateInfo);
} catch (IllegalArgumentException ex) {
logger.warning(String.format(" No crystal information for source structure " + "%s: nothing attached to file %s", sourceFile.toString(), matchFile.toString()));
}
}
String pdb = matchStructure.toPDB();
if (headerLink) {
StringBuilder pdbBuilder = new StringBuilder(pdb);
int position = pdbBuilder.lastIndexOf("REMARK ");
int remarkNumber = 4;
if (position >= 0) {
String nextLine = pdbBuilder.substring(position, position + 1000);
int offset = nextLine.indexOf("%n");
if (offset < 0) {
nextLine = pdbBuilder.substring(position);
offset = nextLine.indexOf("%n");
}
position += offset;
String[] tok = nextLine.split(" +", 3);
try {
remarkNumber = Integer.parseInt(tok[1]) + 1;
} catch (NumberFormatException ex) {
// Silent.
}
}
String toInsert = String.format("REMARK%4d SOURCE FILE: %s", remarkNumber, sourceFile.getName());
toInsert = toInsert.concat(String.format("REMARK%4d RMSD:%11.6f ANGSTROMS", remarkNumber, currentPair.getRMSD()));
pdbBuilder.insert(position, toInsert);
pdb = pdbBuilder.toString();
}
File newFile = createVersionedCopy(matchFile);
try (BufferedWriter bw = new BufferedWriter(new FileWriter(newFile))) {
try {
bw.write(pdb);
} catch (IOException ex) {
logger.warning(String.format(" Error writing to file %s", newFile.getName()));
}
}
}
use of org.biojava.bio.structure.Atom in project ffx by mjschnie.
the class PDBFileMatcher method calculateRMSD.
private double calculateRMSD(FileFilePair matchFile, Structure sourceStructure, StructurePairAligner aligner) throws StructureException {
Structure matchStructure = matchFile.getStructure();
if (superpose) {
double rmsd = Double.MAX_VALUE;
aligner.align(matchStructure, matchStructure);
AlternativeAlignment[] alignments = aligner.getAlignments();
for (AlternativeAlignment alignment : alignments) {
double alignRMSD = alignment.getRmsd();
rmsd = alignRMSD < rmsd ? alignRMSD : rmsd;
}
return rmsd;
}
Atom[] matchArray;
Atom[] sourceArray;
switch(atomsUsed) {
case 1:
String[] atomNames = { "CA", "N1", "N9" };
matchArray = StructureTools.getAtomArray(matchStructure, atomNames);
sourceArray = StructureTools.getAtomArray(sourceStructure, atomNames);
break;
case 2:
List<Atom> matchAtoms = new ArrayList<>();
List<Chain> matchChains = matchStructure.getChains();
for (Chain chain : matchChains) {
List<Group> matchGroups = chain.getAtomGroups(GroupType.AMINOACID);
matchGroups.addAll(chain.getAtomGroups(GroupType.NUCLEOTIDE));
for (Group group : matchGroups) {
matchAtoms.addAll(group.getAtoms());
}
}
matchArray = matchAtoms.toArray(new Atom[matchAtoms.size()]);
List<Atom> sourceAtoms = new ArrayList<>();
List<Chain> sourceChains = sourceStructure.getChains();
for (Chain chain : sourceChains) {
List<Group> sourceGroups = chain.getAtomGroups(GroupType.AMINOACID);
sourceGroups.addAll(chain.getAtomGroups(GroupType.NUCLEOTIDE));
for (Group group : sourceGroups) {
sourceAtoms.addAll(group.getAtoms());
}
}
sourceArray = sourceAtoms.toArray(new Atom[sourceAtoms.size()]);
break;
case 3:
matchAtoms = new ArrayList<>();
matchChains = matchStructure.getChains();
for (Chain chain : matchChains) {
List<Group> matchGroups = chain.getAtomGroups();
for (Group group : matchGroups) {
if (!group.isWater()) {
matchAtoms.addAll(group.getAtoms());
}
}
}
matchArray = matchAtoms.toArray(new Atom[matchAtoms.size()]);
sourceAtoms = new ArrayList<>();
sourceChains = sourceStructure.getChains();
for (Chain chain : sourceChains) {
List<Group> matchGroups = chain.getAtomGroups();
for (Group group : matchGroups) {
if (!group.isWater()) {
sourceAtoms.addAll(group.getAtoms());
}
}
}
sourceArray = sourceAtoms.toArray(new Atom[sourceAtoms.size()]);
break;
case 4:
matchArray = StructureTools.getAllAtomArray(matchStructure);
sourceArray = StructureTools.getAllAtomArray(sourceStructure);
break;
default:
throw new IllegalArgumentException("atomsUsed is not 1-4: this has not been properly checked at an earlier stage.");
}
/*List<Atom> matchAtoms = new ArrayList<>();
Structure matchStructure = matchFile.getStructure();
List<Chain> matchChains = matchStructure.getChains();
if (atomsUsed == 4) {
matchAtoms = StructureTools.getAllAtomArray(matchStructure);
}
for (Chain chain : matchChains) {
List<Group> matchGroups = chain.getSeqResGroups();
for (Group group : matchGroups) {
String groupType = group.getType();
if (groupType.equalsIgnoreCase("HETATOM")) {
if (atomsUsed < 3) {
continue;
} else if (atomsUsed < 4 && group.isWater()) {
continue;
}
}
if (atomsUsed != 1) {
matchAtoms.addAll(group.getAtoms());
} else {
try {
matchAtoms.add(getReferenceAtom(group));
} catch (StructureException ex) {
String refAtomType = (groupType.equalsIgnoreCase("AminoAcid") ? "CA" : "N1/9");
// refAtomType should be binary between amino acid and nucleic acid, as HETATOM is eliminated.
logger.info(String.format(" Reference atom %s could not be found for group %s",
refAtomType, group.toString()));
}
}
}
}
Atom[] match = new Atom[matchAtoms.size()];
matchAtoms.toArray(match);
matchAtoms.clear();
List<Atom> sourceAtoms = new ArrayList<>();
List<Chain> sourceChains = sourceStructure.getChains();
for (Chain chain : sourceChains) {
List<Group> sourceGroups = chain.getSeqResGroups();
for (Group group : sourceGroups) {
String groupType = group.getType();
if (groupType.equalsIgnoreCase("HETATOM")) {
if (atomsUsed < 3) {
continue;
} else if (atomsUsed < 4 && group.isWater()) {
continue;
}
}
if (atomsUsed != 1) {
sourceAtoms.addAll(group.getAtoms());
} else {
try {
sourceAtoms.add(getReferenceAtom(group));
} catch (StructureException ex) {
logger.info(String.format(" Reference atom could not be found for group %s", group.toString()));
}
}
}
}
Atom[] sourceArray = new Atom[sourceAtoms.size()];
sourceAtoms.toArray(sourceArray);
sourceAtoms.clear();*/
return simpleRMSD(sourceArray, matchArray);
}
use of org.biojava.bio.structure.Atom in project ffx by mjschnie.
the class PDBFileMatcher method getMatchingAtom.
/**
* Finds atom1's match in atoms2[], given an array index to search first.
*
* @param atom1 An Atom
* @param atoms2 An Atom[] to search
* @param i Index in atoms2 to check first.
* @return atom1's match in atoms2.
* @throws IllegalArgumentException If no match can be found.
*/
private Atom getMatchingAtom(Atom atom1, Atom[] atoms2, int i) throws IllegalArgumentException {
Atom atom2;
try {
atom2 = atoms2[i];
} catch (ArrayIndexOutOfBoundsException ex) {
// May be important if one structure lacks hydrogens.
atom2 = atoms2[0];
}
if (compareAtoms(atom1, atom2)) {
return atom2;
}
atom2 = getMatchingAtom(atom1, atoms2);
return atom2;
}
Aggregations