001/* 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 */ 021 022package org.biojava.nbio.structure.symmetry.core; 023 024import org.biojava.nbio.structure.symmetry.geometry.SuperPosition; 025 026import javax.vecmath.Matrix4d; 027import javax.vecmath.Point3d; 028import java.util.List; 029 030/** 031 * 032 * @author Peter 033 */ 034public class QuatSuperpositionScorer { 035 036 /** 037 * Returns minimum, mean, and maximum RMSD and TM-Score for two superimposed sets of subunits 038 * 039 * TM score: Yang Zhang and Jeffrey Skolnick, PROTEINS: Structure, Function, and Bioinformatics 57:702–710 (2004) 040 * @param subunits subunits to be scored 041 * @param transformation transformation matrix 042 * @param permutations permutation that determines which subunits are superposed 043 * @return 044 */ 045 public static QuatSymmetryScores calcScores(Subunits subunits, Matrix4d transformation, List<Integer> permutation) { 046 QuatSymmetryScores scores = new QuatSymmetryScores(); 047 048 double minTm = Double.MAX_VALUE; 049 double maxTm = Double.MIN_VALUE; 050 double minRmsd = Double.MAX_VALUE; 051 double maxRmsd = Double.MIN_VALUE; 052 053 double totalSumTm = 0; 054 double totalSumDsq = 0; 055 double totalLength = 0; 056 057 Point3d t = new Point3d(); 058 List<Point3d[]> traces = subunits.getTraces(); 059 060 // loop over the Calpha atoms of all subunits 061 for (int i = 0; i < traces.size(); i++) { 062 // in helical systems not all permutations involve all subunit. -1 indicates subunits that should not be permuted. 063 if (permutation.get(i) == -1) { 064 continue; 065 } 066 // get original subunit 067 Point3d[] orig = traces.get(i); 068 totalLength += orig.length; 069 070 // get permuted subunit 071 Point3d[] perm = traces.get(permutation.get(i)); 072 073 // calculate TM specific parameters 074 int tmLen = Math.max(orig.length, 17); // don't let d0 get negative with short sequences 075 double d0 = 1.24 * Math.cbrt(tmLen - 15.0) - 1.8; 076 double d0Sq = d0 * d0; 077 078 double sumTm = 0; 079 double sumDsq = 0; 080 for (int j = 0; j < orig.length; j++) { 081 // transform coordinates of the permuted subunit 082 t.set(perm[j]); 083 transformation.transform(t); 084 085 double dSq = orig[j].distanceSquared(t); 086 sumTm += 1.0/(1.0 + dSq/d0Sq); 087 sumDsq += dSq; 088 } 089 090 // scores for individual subunits 091 double sTm = sumTm/tmLen; 092 minTm = Math.min(minTm, sTm); 093 maxTm = Math.max(maxTm, sTm); 094 095 double sRmsd = Math.sqrt(sumDsq/orig.length); 096 minRmsd = Math.min(minRmsd, sRmsd); 097 maxRmsd = Math.max(maxRmsd, sRmsd); 098 099 totalSumTm += sumTm; 100 totalSumDsq += sumDsq; 101 } 102 103 // save scores for individual subunits 104 scores.setMinRmsd(minRmsd); 105 scores.setMaxRmsd(maxRmsd); 106 scores.setMinTm(minTm); 107 scores.setMaxTm(maxTm); 108 109 // save mean scores over all subunits 110 scores.setTm(totalSumTm/totalLength); 111 scores.setRmsd(Math.sqrt(totalSumDsq/totalLength)); 112 113 // add intra subunit scores 114 calcIntrasubunitScores(subunits, transformation, permutation, scores); 115 116 return scores; 117 } 118 119 private static void calcIntrasubunitScores(Subunits subunits, Matrix4d transformation, List<Integer> permutation, QuatSymmetryScores scores) { 120 double totalSumTm = 0; 121 double totalSumDsq = 0; 122 double totalLength = 0; 123 124 List<Point3d[]> traces = subunits.getTraces(); 125 126 // loop over the Calpha atoms of all subunits 127 for (int i = 0; i < traces.size(); i++) { 128 // in helical systems not all permutations involve all subunit. -1 indicates subunits that should not be permuted. 129 if (permutation.get(i) == -1) { 130 continue; 131 } 132 // get original subunit 133 Point3d[] orig = traces.get(i); 134 totalLength += orig.length; 135 136 // get permuted subunit 137 Point3d[] perm = traces.get(permutation.get(i)); 138 139 // calculate TM specific parameters 140 int tmLen = Math.max(orig.length, 17); // don't let d0 get negative with short sequences 141 double d0 = 1.24 * Math.cbrt(tmLen - 15.0) - 1.8; 142 double d0Sq = d0 * d0; 143 144 double sumTm = 0; 145 double sumDsq = 0; 146 Point3d[] trans = new Point3d[orig.length]; 147 for (int j = 0; j < orig.length; j++) { 148 trans[j] = new Point3d(perm[j]); 149 } 150 151 // superpose individual subunits 152 SuperPosition.superposeWithTranslation(trans, orig); 153 for (int j = 0; j < orig.length; j++) { 154 double dSq = orig[j].distanceSquared(trans[j]); 155 sumTm += 1.0/(1.0 + dSq/d0Sq); 156 sumDsq += dSq; 157 } 158 159 totalSumTm += sumTm; 160 totalSumDsq += sumDsq; 161 } 162 // save mean scores over all subunits 163 scores.setRmsdIntra(Math.sqrt(totalSumDsq/totalLength)); 164 scores.setTmIntra(totalSumTm/totalLength); 165 } 166 167}