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.geometry.SuperPositions; 025 026import javax.vecmath.Matrix4d; 027import javax.vecmath.Point3d; 028 029import java.util.List; 030 031/** 032 * 033 * @author Peter 034 */ 035public class QuatSuperpositionScorer { 036 037 /** 038 * Returns minimum, mean, and maximum RMSD and TM-Score for two superimposed sets of subunits 039 * 040 * TM score: Yang Zhang and Jeffrey Skolnick, PROTEINS: Structure, Function, and Bioinformatics 57:702–710 (2004) 041 * @param subunits subunits to be scored 042 * @param transformation transformation matrix 043 * @param permutations permutation that determines which subunits are superposed 044 * @return 045 */ 046 public static QuatSymmetryScores calcScores(QuatSymmetrySubunits subunits, Matrix4d transformation, List<Integer> permutation) { 047 QuatSymmetryScores scores = new QuatSymmetryScores(); 048 049 double minTm = Double.MAX_VALUE; 050 double maxTm = Double.MIN_VALUE; 051 double minRmsd = Double.MAX_VALUE; 052 double maxRmsd = Double.MIN_VALUE; 053 054 double totalSumTm = 0; 055 double totalSumDsq = 0; 056 double totalLength = 0; 057 058 Point3d t = new Point3d(); 059 List<Point3d[]> traces = subunits.getTraces(); 060 061 // loop over the Calpha atoms of all subunits 062 for (int i = 0; i < traces.size(); i++) { 063 // in helical systems not all permutations involve all subunit. -1 indicates subunits that should not be permuted. 064 if (permutation.get(i) == -1) { 065 continue; 066 } 067 // get original subunit 068 Point3d[] orig = traces.get(i); 069 totalLength += orig.length; 070 071 // get permuted subunit 072 Point3d[] perm = traces.get(permutation.get(i)); 073 074 // calculate TM specific parameters 075 int tmLen = Math.max(orig.length, 17); // don't let d0 get negative with short sequences 076 double d0 = 1.24 * Math.cbrt(tmLen - 15.0) - 1.8; 077 double d0Sq = d0 * d0; 078 079 double sumTm = 0; 080 double sumDsq = 0; 081 for (int j = 0; j < orig.length; j++) { 082 // transform coordinates of the permuted subunit 083 t.set(perm[j]); 084 transformation.transform(t); 085 086 double dSq = orig[j].distanceSquared(t); 087 sumTm += 1.0/(1.0 + dSq/d0Sq); 088 sumDsq += dSq; 089 } 090 091 // scores for individual subunits 092 double sTm = sumTm/tmLen; 093 minTm = Math.min(minTm, sTm); 094 maxTm = Math.max(maxTm, sTm); 095 096 double sRmsd = Math.sqrt(sumDsq/orig.length); 097 minRmsd = Math.min(minRmsd, sRmsd); 098 maxRmsd = Math.max(maxRmsd, sRmsd); 099 100 totalSumTm += sumTm; 101 totalSumDsq += sumDsq; 102 } 103 104 // save scores for individual subunits 105 scores.setMinRmsd(minRmsd); 106 scores.setMaxRmsd(maxRmsd); 107 scores.setMinTm(minTm); 108 scores.setMaxTm(maxTm); 109 110 // save mean scores over all subunits 111 scores.setTm(totalSumTm/totalLength); 112 scores.setRmsd(Math.sqrt(totalSumDsq/totalLength)); 113 114 // add intra subunit scores 115 calcIntrasubunitScores(subunits, transformation, permutation, scores); 116 117 return scores; 118 } 119 120 private static void calcIntrasubunitScores(QuatSymmetrySubunits subunits, Matrix4d transformation, List<Integer> permutation, QuatSymmetryScores scores) { 121 double totalSumTm = 0; 122 double totalSumDsq = 0; 123 double totalLength = 0; 124 125 List<Point3d[]> traces = subunits.getTraces(); 126 127 // loop over the Calpha atoms of all subunits 128 for (int i = 0; i < traces.size(); i++) { 129 // in helical systems not all permutations involve all subunit. -1 indicates subunits that should not be permuted. 130 if (permutation.get(i) == -1) { 131 continue; 132 } 133 // get original subunit 134 Point3d[] orig = traces.get(i); 135 totalLength += orig.length; 136 137 // get permuted subunit 138 Point3d[] perm = traces.get(permutation.get(i)); 139 140 // calculate TM specific parameters 141 int tmLen = Math.max(orig.length, 17); // don't let d0 get negative with short sequences 142 double d0 = 1.24 * Math.cbrt(tmLen - 15.0) - 1.8; 143 double d0Sq = d0 * d0; 144 145 double sumTm = 0; 146 double sumDsq = 0; 147 Point3d[] trans = new Point3d[orig.length]; 148 for (int j = 0; j < orig.length; j++) { 149 trans[j] = new Point3d(perm[j]); 150 } 151 152 // superpose individual subunits 153 SuperPositions.superposeAndTransform(orig, trans); 154 for (int j = 0; j < orig.length; j++) { 155 double dSq = orig[j].distanceSquared(trans[j]); 156 sumTm += 1.0/(1.0 + dSq/d0Sq); 157 sumDsq += dSq; 158 } 159 160 totalSumTm += sumTm; 161 totalSumDsq += sumDsq; 162 } 163 // save mean scores over all subunits 164 scores.setRmsdIntra(Math.sqrt(totalSumDsq/totalLength)); 165 scores.setTmIntra(totalSumTm/totalLength); 166 } 167 168}