Search in sources :

Example 1 with Atom

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;
}
Also used : Atom(org.biojava.bio.structure.Atom)

Example 2 with Atom

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.");
}
Also used : Chain(org.biojava.bio.structure.Chain) Group(org.biojava.bio.structure.Group) StructureException(org.biojava.bio.structure.StructureException) ResidueNumber(org.biojava.bio.structure.ResidueNumber) Atom(org.biojava.bio.structure.Atom)

Example 3 with Atom

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()));
        }
    }
}
Also used : FileWriter(java.io.FileWriter) PDBCrystallographicInfo(org.biojava.bio.structure.PDBCrystallographicInfo) IOException(java.io.IOException) Atom(org.biojava.bio.structure.Atom) BufferedWriter(java.io.BufferedWriter) SSBond(org.biojava.bio.structure.SSBond) Structure(org.biojava.bio.structure.Structure) File(java.io.File)

Example 4 with Atom

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);
}
Also used : Chain(org.biojava.bio.structure.Chain) Group(org.biojava.bio.structure.Group) ArrayList(java.util.ArrayList) Atom(org.biojava.bio.structure.Atom) AlternativeAlignment(org.biojava.bio.structure.align.pairwise.AlternativeAlignment) Structure(org.biojava.bio.structure.Structure)

Example 5 with Atom

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;
}
Also used : Atom(org.biojava.bio.structure.Atom)

Aggregations

Atom (org.biojava.bio.structure.Atom)6 Structure (org.biojava.bio.structure.Structure)3 Chain (org.biojava.bio.structure.Chain)2 Group (org.biojava.bio.structure.Group)2 BufferedWriter (java.io.BufferedWriter)1 File (java.io.File)1 FileWriter (java.io.FileWriter)1 IOException (java.io.IOException)1 ArrayList (java.util.ArrayList)1 PDBCrystallographicInfo (org.biojava.bio.structure.PDBCrystallographicInfo)1 ResidueNumber (org.biojava.bio.structure.ResidueNumber)1 SSBond (org.biojava.bio.structure.SSBond)1 StructureException (org.biojava.bio.structure.StructureException)1 AlternativeAlignment (org.biojava.bio.structure.align.pairwise.AlternativeAlignment)1