001/* This class is based on the original FATCAT implementation by
002 * <pre>
003 * Yuzhen Ye & Adam Godzik (2003)
004 * Flexible structure alignment by chaining aligned fragment pairs allowing twists.
005 * Bioinformatics vol.19 suppl. 2. ii246-ii255.
006 * https://www.ncbi.nlm.nih.gov/pubmed/14534198
007 * </pre>
008 *
009 * Thanks to Yuzhen Ye and A. Godzik for granting permission to freely use and redistribute this code.
010 *
011 * This code may be freely distributed and modified under the
012 * terms of the GNU Lesser General Public Licence.  This should
013 * be distributed with the code.  If you do not have a copy,
014 * see:
015 *
016 *      http://www.gnu.org/copyleft/lesser.html
017 *
018 * Copyright for this code is held jointly by the individual
019 * authors.  These should be listed in @author doc comments.
020 *
021 *
022 * Created on Jun 17, 2009
023 * Created by Andreas Prlic - RCSB PDB
024 *
025 */
026
027package org.biojava.nbio.structure.align.fatcat.calc;
028
029import org.biojava.nbio.structure.*;
030import org.biojava.nbio.structure.align.model.AFP;
031import org.biojava.nbio.structure.align.model.AFPChain;
032import org.biojava.nbio.structure.jama.Matrix;
033
034import java.util.ArrayList;
035import java.util.List;
036
037/** a class that performs calculations on AFPCHains
038 *
039 * @author Andreas Prlic
040 *
041 */
042public class AFPCalculator
043{
044        public static final boolean debug = FatCatAligner.debug;
045
046
047        public static final  void extractAFPChains(FatCatParameters params, AFPChain afpChain,Atom[] ca1,Atom[] ca2) throws StructureException {
048
049
050
051                List<AFP> afpSet = new ArrayList<AFP>();
052                afpChain.setAfpSet(afpSet);
053
054                if ( debug )
055                        System.err.println("nr of atoms ca1: " + ca1.length + " ca2: " +  ca2.length);
056
057
058
059                int     p1, p2;
060                @SuppressWarnings("unused")
061                int n0, n, n1, n2;
062                double  filter1;
063                double rmsd = 0;
064
065                Matrix r = new Matrix(3,3);
066                Atom   t = new AtomImpl();
067
068
069                int sparse = params.getSparse();
070                int maxTra = params.getMaxTra();
071                int fragLen = params.getFragLen();
072                double disFilter = params.getDisFilter();
073                double rmsdCut = params.getRmsdCut();
074                double badRmsd = params.getBadRmsd();
075                double fragScore = params.getFragScore();
076
077                int     add = sparse + 1; //if add > 1, use sparse sampling
078                n0 = n = n1 = n2 = 0;
079
080                int minLen = 0;
081
082                int prot1Length = ca1.length;
083                int prot2Length = ca2.length;
084
085                if(prot1Length < prot2Length)
086                        minLen = prot1Length;
087                else
088                        minLen = prot2Length;
089                afpChain.setMinLen(minLen);
090
091                afpChain.setBlockResList(new int[maxTra+1][2][minLen]);
092                afpChain.setFocusRes1(new int[minLen]);
093                afpChain.setFocusRes2(new int[minLen]);
094
095                for(p1 = 0; p1 < prot1Length - fragLen; p1 += add )    {
096                        for(p2 = 0; p2 < prot2Length - fragLen; p2 += add)     {
097                                n0 ++;
098                                filter1 = getEnd2EndDistance(ca1, ca2, p1, p1 + fragLen - 1, p2, p2 + fragLen - 1);
099                                //difference bewteen end-to-end distances
100                                if(filter1 > disFilter) { n1 ++; continue; }
101                                boolean filter2 = filterTerminal(ca1,ca2, p1, p1 + fragLen - 1, p2, p2 + fragLen - 1, fragLen, minLen);
102                                if(filter2)     {
103                                        n2 ++;
104                                        continue;
105
106                                } //be cautious to use this filter !!
107
108                                // here FATCAT does a a jacobi transformation
109                                //rmsd = kearsay(fragLen, ca1[p1], ca2[p2], r, t);
110                                // we use the BioJava SVD instead...
111
112                                //
113                                rmsd = getRmsd(ca1,ca2,fragLen, p1,p2,r,t);
114
115                                if(rmsd < rmsdCut)      {
116                                        AFP     afptmp = new AFP();
117                                        afptmp.setP1(p1);
118                                        afptmp.setP2(p2);
119                                        afptmp.setFragLen(fragLen);
120                                        afptmp.setRmsd(rmsd);
121                                        afptmp.setM(r);
122                                        afptmp.setT(t.getCoords());
123                                        afptmp.setScore(scoreAfp(afptmp,badRmsd,fragScore));
124                                        afpSet.add(afptmp);
125                                        n ++;
126                                }
127                        }
128                }
129
130                int afpNum = afpSet.size();
131
132                if(debug) {
133                        String msg = String.format("possible AFP-pairs %d, remain %d after filter 1 remove %d; filter 2 remove %d\n",
134                                        n0, afpNum, n1, n2);
135                        System.err.println(msg);
136                }
137
138
139        }
140
141        /**
142         * filter 1 for AFP extration: the distance of end-to-end
143         * @param p1b
144         * @param p1e
145         * @param p2b
146         * @param p2e
147         * @return
148         */
149        private static final double getEnd2EndDistance(Atom[] ca1, Atom[] ca2, int p1b, int p1e, int p2b, int p2e)
150        {
151
152                double min = 99;
153                        double dist1 = Calc.getDistance(ca1[p1b], ca1[p1e]);
154                        double dist2 = Calc.getDistance(ca2[p2b], ca2[p2e]);
155                        min = dist1 - dist2;
156
157                return Math.abs(min);
158        }
159
160        /**
161         * filter 2 for AFP extration: the context
162         * @param p1b
163         * @param p1e
164         * @param p2b
165         * @param p2e
166         * @return
167         */
168
169        private static final  boolean filterTerminal(Atom[] ca1, Atom[] ca2, int p1b, int p1e, int p2b, int p2e, int fragLen, int minLen)
170        {
171                int     d1 = (p1b < p2b)?p1b:p2b;
172                int     d2 = (ca1.length - p1e) < (ca2.length - p2e)?(ca1.length - p1e):(ca2.length - p2e);
173                int     d3 = d1 + d2 + fragLen; //maximum alignment length from current AFP
174
175
176                /// DO NOT DO Math.round() this will give different results to FATCAT....
177                int     d4 = (int)(0.3 * minLen);
178
179                return d3 < d4;
180
181        }
182
183        private static final double getRmsd(Atom[] ca1, Atom[] ca2, int fragLen, int p1, int p2, Matrix m, Atom t) throws StructureException {
184
185
186                double rmsd = 99.9;
187                        Atom[] catmp1 = getFragment(ca1, p1, fragLen,false);
188                        Atom[] catmp2 = getFragment(ca2, p2, fragLen,true); // clone the atoms for fragment 2 so we can manipulate them...
189
190                        if ( catmp1 == null) {
191                                System.err.println("could not get fragment for ca1 " + p1 + " " + fragLen );
192                                return rmsd;
193
194                        }
195
196                        if ( catmp2 == null) {
197                                System.err.println("could not get fragment for ca2 " + p2 + " " + fragLen );
198                                return rmsd;
199
200                        }
201
202                        SVDSuperimposer svd = new SVDSuperimposer(catmp1, catmp2);
203
204                        m = svd.getRotation();
205                        t = svd.getTranslation();
206
207                        for (Atom a : catmp2){
208                                Calc.rotate(a,m);
209                                Calc.shift(a,t);
210
211                        }
212
213                        rmsd = SVDSuperimposer.getRMS(catmp1,catmp2);
214
215                return rmsd;
216        }
217
218        /** get a continue subset of Atoms based by the starting position and the length
219         *
220         * @param caall
221         * @param pos ... the start position
222         * @param fragmentLength .. the length of the subset to extract.
223         * @param clone: returns a copy of the atom (in case the coordinate get manipulated...)
224         * @return an Atom[] array
225         */
226        private static final Atom[] getFragment(Atom[] caall, int pos, int fragmentLength , boolean clone){
227
228                if ( pos+fragmentLength > caall.length)
229                        return null;
230
231                Atom[] tmp = new Atom[fragmentLength];
232
233                for (int i=0;i< fragmentLength;i++){
234                        if (clone){
235                                tmp[i] = (Atom)caall[i+pos].clone();
236                        } else {
237                                tmp[i] = caall[i+pos];
238                        }
239                }
240                return tmp;
241
242        }
243
244
245        /**
246         * Assign score to each AFP
247         */
248
249        private static final double scoreAfp(AFP afp, double badRmsd, double fragScore)
250        {
251                //longer AFP with low rmsd is better
252                double  s, w;
253                //s = (rmsdCut - afptmp.rmsd) * afptmp.len; //the same scroing strategy as that in the post-processing
254                w = afp.getRmsd() / badRmsd;
255                w = w * w;
256                s = fragScore * (1.0 - w);
257                return s;
258        }
259
260        //------------------------------------------------------------------
261        //Sort the AFPs in increase of their diagonals(i,j)
262        //------------------------------------------------------------------
263        public static final  void sortAfps(AFPChain afpChain, Atom[] ca1, Atom[] ca2)
264        {
265
266
267                List<AFP> afpSet = afpChain.getAfpSet();
268
269                if ( debug)
270                        System.err.println("entering sortAfps");
271
272                int pro1Len = ca1.length;
273                int pro2Len = ca2.length;
274
275                afpChain.setAfpIndex(      new int[pro1Len][pro2Len]); //the index of (i,j) pair in AFP list, otherwise -1
276                afpChain.setAfpAftIndex(   new int[pro1Len][pro2Len]);  //the index of AFP (i,j*) nearest to (i,j), j*<j. if a AFP exits for (i,j), it equals to afpIndex
277                afpChain.setAfpBefIndex(   new int[pro1Len][pro2Len]); //the index of AFP (i,j*) nearest to (i,j), j*>j. if a AFP exits for (i,j), it equals to afpIndex
278
279                int[][] afpIndex       = afpChain.getAfpIndex();
280                int[][] afpAftIndex    = afpChain.getAfpAftIndex();
281                int[][] afpBefIndex    = afpChain.getAfpBefIndex();
282
283                for(int i = 0; i < pro1Len; i ++)   {
284                        for(int j = 0; j < pro2Len; j ++)   {
285
286                                afpIndex[i][j] = afpAftIndex[i][j] = afpBefIndex[i][j] = -1;
287                        }
288                }
289
290                //index the AFP for easy extraction of compatible AFPs
291                int afpNum = afpSet.size();
292
293                int     b0 = 0;
294                for(int a = 0; a < afpNum; a ++)    {
295                        if(a == afpNum - 1 || afpSet.get(a).getP1() != afpSet.get(a+1).getP1())   {
296                                int i = afpSet.get(a).getP1();
297                                for(int b = b0; b <= a; b ++)       {
298                                        int j = afpSet.get(b).getP2();
299                                        afpIndex[i][j]=b ;
300                                        afpBefIndex[i][j]=b;
301                                        afpAftIndex[i][j]=b;
302                                        if(afpSet.get(b).getP1() != i)    {
303                                                System.err.println(String.format("Warning: wrong afp index %d %d\n", i, afpSet.get(b).getP1()));
304                                                return;
305                                        }
306                                }
307                                for(int k = 1; k < pro2Len; k ++)   {
308                                        if( afpBefIndex[i][k] == -1){
309                                                afpBefIndex[i][k] = afpBefIndex[i][k-1];
310                                        }
311                                }
312                                for(int k = pro2Len - 2; k >= 0; k --)      {
313                                        if(afpAftIndex[i][k] == -1) {
314                                                afpAftIndex[i][k] =  afpAftIndex[i][k+1];
315                                        }
316                                }
317                                b0 = a + 1;
318                        }
319                }
320
321                if ( debug)
322                        System.err.println("done sortAfps");
323
324
325        }
326}