Source code

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 * Created on Jan 28, 2006
021 *
022 */
023package org.biojava.nbio.structure.align.pairwise;
024
025import org.biojava.nbio.structure.Atom;
026import org.biojava.nbio.structure.Calc;
027import org.biojava.nbio.structure.SVDSuperimposer;
028import org.biojava.nbio.structure.StructureException;
029import org.biojava.nbio.structure.align.StrucAligParameters;
030import org.biojava.nbio.structure.align.helper.AlignTools;
031import org.biojava.nbio.structure.align.helper.JointFragments;
032import org.biojava.nbio.structure.jama.Matrix;
033
034import java.util.ArrayList;
035import java.util.Collections;
036import java.util.Comparator;
037import java.util.List;
038import java.util.logging.Logger;
039
040
041/** Joins the initial Fragments together to larger Fragments
042 *
043 * @author Andreas Prlic
044 * @author Peter Lackner
045 * @since 1.5
046 * @version %I% %G%
047 */
048public class FragmentJoiner {
049
050        public static Logger logger =  Logger.getLogger("org.biojava.nbio.structure.align");
051
052        public FragmentJoiner() {
053                super();
054
055        }
056
057        /**
058         * Reallocates an array with a new size, and copies the contents
059         * of the old array to the new array.
060         * @param oldArray  the old array, to be reallocated.
061         * @param newSize   the new array size.
062         * @return          A new array with the same contents.
063         */
064        @SuppressWarnings("rawtypes")
065        public static Object resizeArray (Object oldArray, int newSize) {
066                int oldSize = java.lang.reflect.Array.getLength(oldArray);
067                Class elementType = oldArray.getClass().getComponentType();
068                Object newArray = java.lang.reflect.Array.newInstance(
069                                  elementType,newSize);
070                int preserveLength = Math.min(oldSize,newSize);
071                if (preserveLength > 0)
072                        System.arraycopy (oldArray,0,newArray,0,preserveLength);
073                return newArray;
074        }
075
076
077        /**
078         *  In helices often many similar fragments can be found. To reduce these to a few
079         *  representative ones this check can be used. It does a distance check between
080         *  all known Fragments and a new one. If this one is on a similar diagonal and it
081         *  has a lower rms, this one is a better representation. Note: shifts of one are
082         *  not allowed.
083         *
084         * @param fragments
085         * @param f
086         * @param rmsmat
087         * @return true - if this is a better representant for a group of locala fragments.
088         */
089        public static boolean reduceFragments(List<FragmentPair> fragments, FragmentPair f, Matrix rmsmat){
090                boolean doNotAdd =false;
091                int i = f.getPos1();
092                int j = f.getPos2();
093
094                for ( int p =0; p < fragments.size(); p++){
095                        FragmentPair tmp = fragments.get(p);
096
097                        int di1 = Math.abs(f.getPos1() - tmp.getPos1());
098                        int di2 = Math.abs(f.getPos2() - tmp.getPos2());
099                        if (( Math.abs(di1-di2) == 2)) {
100                                double rms1 = rmsmat.get(tmp.getPos1(),tmp.getPos2());
101                                double rms2 = rmsmat.get(i,j);
102                                doNotAdd = true;
103                                if ( rms2 < rms1){
104                                        fragments.remove(p);
105                                        fragments.add(f);
106                                        break;
107                                }
108                                p = fragments.size();
109                        }
110                }
111                return doNotAdd;
112        }
113
114
115        public JointFragments[] approach_ap3(Atom[] ca1, Atom[] ca2, FragmentPair[] fraglst,
116                                                                                                         StrucAligParameters params) throws StructureException {
117
118                //the final list of joined fragments stores as apairs
119                List<JointFragments> fll = new ArrayList<JointFragments>();
120
121                FragmentPair[] tmpfidx = new FragmentPair[fraglst.length];
122                for ( int i=0 ; i < fraglst.length; i++){
123                        tmpfidx[i] = (FragmentPair)fraglst[i].clone();
124                }
125
126                int n  = tmpfidx.length;
127
128                // if two fragments after having been joint have rms < X
129                // keep the new fragment.
130
131                for (int i=0;i< fraglst.length;i++){
132
133                        boolean[] used = new boolean[n];
134
135                        int p1i  = tmpfidx[i].getPos1();
136                        int p1j  = tmpfidx[i].getPos2();
137                        int l1   = tmpfidx[i].getLength();
138
139                        JointFragments f = new JointFragments();
140
141                        int maxi=p1i+l1-1;
142
143                        f.add(p1i,p1j,0,l1);
144                        used[i] = true;
145
146                        //n = tmpfidx.length;
147
148                        for (int j=(i+1);j<n;j++){
149
150                                if ( used[j])
151                                        continue;
152
153                                int p2i = tmpfidx[j].getPos1();
154                                int p2j = tmpfidx[j].getPos2();
155                                int l2  = tmpfidx[j].getLength();
156
157                                if ( p2i > maxi) {
158
159
160                                        // TODO: replace this with plo angle calculation
161                                        if ( params.isDoAngleCheck()){
162
163                                                // was 0.174
164                                                if (! angleCheckOk(tmpfidx[i], tmpfidx[j], 0.4f))
165                                                        continue;
166                                        }
167
168                                        if ( params.isDoDistanceCheck()) {
169
170                                                if (! distanceCheckOk(tmpfidx[i],tmpfidx[j], params.getFragmentMiniDistance()))
171                                                        continue;
172                                        }
173
174                                        if ( params.isDoDensityCheck()) {
175
176                                                if ( ! densityCheckOk(ca1,ca2, f.getIdxlist(), p2i,p2j, l2, params.getDensityCutoff()))
177                                                        continue;
178                                        }
179
180
181                                        if ( params.isDoRMSCheck()) {
182
183                                                double rms = rmsCheck(ca1,ca2, f.getIdxlist(), p2i, p2j, l2);
184                                                if ( rms > params.getJoinRMSCutoff())
185                                                        continue;
186                                                f.setRms(rms);
187                                        }
188
189                                        f.add(p2i,p2j,0,l2);
190                                        used[j] = true;
191                                        maxi = p2i+l2-1;
192
193
194                                }
195                        }
196                        fll.add(f);
197                }
198
199
200                Comparator<JointFragments> comp = new JointFragmentsComparator();
201                Collections.sort(fll,comp);
202                Collections.reverse(fll);
203                int m = Math.min(params.getMaxrefine(),fll.size());
204                List<JointFragments> retlst = new ArrayList<JointFragments>();
205                for ( int i = 0 ; i < m ; i++){
206                        JointFragments jf = fll.get(i);
207                        retlst.add(jf);
208                }
209
210                return retlst.toArray(new JointFragments[retlst.size()]);
211
212        }
213
214        private boolean densityCheckOk(Atom[] aa1, Atom[] aa2, List<int[]> idxlist,
215                                                                                         int p2i, int p2j, int l2,
216                                                                                         double densityCutoff) throws StructureException {
217                JointFragments ftmp = new JointFragments();
218                ftmp.setIdxlist(idxlist);
219                ftmp.add(p2i,p2j,0,l2);
220                AlternativeAlignment ali = new AlternativeAlignment();
221                ali.apairs_from_idxlst(ftmp);
222
223                Atom[] aa3 = aa2.clone();
224
225                int[] idx1 = ali.getIdx1();
226                int[] idx2 = ali.getIdx2();
227
228                Atom[] ca1subset = AlignTools.getFragmentFromIdxList(aa1, idx1);
229
230                Atom[] ca2subset = AlignTools.getFragmentFromIdxList(aa3,idx2);
231
232                double density = getDensity(ca1subset, ca2subset);
233
234                return density <= densityCutoff;
235
236
237        }
238
239
240        /** this is probably useless
241         *
242         * @param ca1subset
243         * @param ca2subset
244         * @return a double
245         * @throws StructureException
246         */
247        private double getDensity(Atom[] ca1subset, Atom[] ca2subset ) throws StructureException{
248
249                Atom centroid1 =  Calc.getCentroid(ca1subset);
250                Atom centroid2 = Calc.getCentroid(ca2subset);
251
252                // get Average distance to centroid ...
253
254                double d1 = 0;
255                double d2 = 0;
256
257                for ( int i = 0 ; i < ca1subset.length;i++){
258                        double dd1 = Calc.getDistance(centroid1, ca1subset[i]);
259                        double dd2 = Calc.getDistance(centroid2, ca2subset[i]);
260
261                        d1 += dd1;
262                        d2 += dd2;
263
264                }
265
266                double avd1 = d1 / ca1subset.length;
267                double avd2 = d2 / ca2subset.length;
268
269                return Math.min(avd1,avd2);
270        }
271
272        private double rmsCheck(Atom[] a1, Atom[] a2,List<int[]> idxlist, int p2i, int p2j, int l2) throws StructureException {
273
274                //System.out.println(dd);
275                // check if a joint fragment has low rms ...
276                JointFragments ftmp = new JointFragments();
277                ftmp.setIdxlist(idxlist);
278                ftmp.add(p2i,p2j,0,l2);
279                Atom[] a3 = new Atom[a2.length];
280                for (int i=0;i < a2.length;i++){
281                        a3[i] = (Atom)a2[i].clone();
282                }
283                return getRMS(a1,a3,ftmp);
284        }
285
286        /** get the RMS of the JointFragments pair frag
287         *
288         * @param ca1 the array of all atoms of structure1
289         * @param ca2 the array of all atoms of structure1
290         * @param frag the JointFragments object that contains the list of identical positions
291         * @return the rms
292         */
293        public static double getRMS(Atom[] ca1, Atom[]ca2,JointFragments frag) throws StructureException {
294                //      now svd ftmp and check if the rms is < X ...
295                AlternativeAlignment ali = new AlternativeAlignment();
296                ali.apairs_from_idxlst(frag);
297                double rms = 999;
298
299                int[] idx1 = ali.getIdx1();
300                int[] idx2 = ali.getIdx2();
301
302                Atom[] ca1subset = AlignTools.getFragmentFromIdxList(ca1, idx1);
303
304                Atom[] ca2subset = AlignTools.getFragmentFromIdxList(ca2,idx2);
305
306                ali.calculateSuperpositionByIdx(ca1,ca2);
307
308                Matrix rot = ali.getRotationMatrix();
309                Atom atom = ali.getShift();
310
311                for (Atom a : ca2subset) {
312                        Calc.rotate(a, rot);
313                        Calc.shift(a, atom);
314                }
315
316                rms = SVDSuperimposer.getRMS(ca1subset,ca2subset);
317
318                return rms;
319        }
320
321        public boolean angleCheckOk(FragmentPair a, FragmentPair b, float distcutoff){
322
323                double dist = -999;
324
325                Atom v1 = a.getUnitv();
326                Atom v2 = b.getUnitv();
327                dist = Calc.getDistance(v1,v2);
328
329                return dist <= distcutoff;
330        }
331
332        private boolean distanceCheckOk(FragmentPair a, FragmentPair b, float fragCompatDist){
333
334                double dd ;
335
336                Atom c1i = a.getCenter1();
337                Atom c1j = b.getCenter1();
338                Atom c2i = a.getCenter2();
339                Atom c2j = b.getCenter2();
340                dd = Calc.getDistance(c1i,c1j) - Calc.getDistance(c2i,c2j);
341
342
343                if ( dd < 0) dd = -dd;
344                return dd <= fragCompatDist;
345
346        }
347
348        /**
349         * Calculate the pairwise compatibility of fpairs.
350
351         Iterates through a list of fpairs and joins them if
352         they have compatible rotation and translation parameters.
353
354
355         * @param fraglst FragmentPair[] array
356         * @param angleDiff angle cutoff
357         * @param fragCompatDist distance cutoff
358         * @param maxRefine max number of solutions to keep
359         * @return JointFragments[]
360         */
361        public JointFragments[] frag_pairwise_compat(FragmentPair[] fraglst, int angleDiff, float fragCompatDist, int maxRefine  ){
362
363
364                FragmentPair[] tmpfidx = new FragmentPair[fraglst.length];
365                for ( int i=0 ; i < fraglst.length; i++){
366                        tmpfidx[i] = (FragmentPair)fraglst[i].clone();
367                }
368
369                int n  = tmpfidx.length;
370
371                //indicator if a fragment was already used
372                int[] used = new int[n];
373
374                //the final list of joined fragments stores as apairs
375                List<JointFragments> fll = new ArrayList<JointFragments>();
376
377                double adiff = angleDiff * Math.PI / 180d;
378                logger.finer("addiff" + adiff);
379                //distance between two unit vectors with angle adiff
380                double ddiff = Math.sqrt(2.0-2.0*Math.cos(adiff));
381                logger.finer("ddiff" + ddiff);
382
383                // the fpairs in the flist have to be sorted with respect to their positions
384
385                while (tmpfidx.length > 0){
386                        int i = 0;
387                        int j = 1;
388                        used[i]=1;
389                        JointFragments f = new JointFragments();
390
391                        int p1i  = tmpfidx[i].getPos1();
392                        int p1j  = tmpfidx[i].getPos2();
393
394                        int maxi = p1i+tmpfidx[i].getLength();
395
396                        f.add(p1i,p1j,0,tmpfidx[i].getLength());
397
398                        n = tmpfidx.length;
399
400                        while ( j < n) {
401                                int p2i = tmpfidx[j].getPos1();
402                                int p2j = tmpfidx[j].getPos2();
403                                int l2  = tmpfidx[j].getLength();
404                                if ( p2i > maxi) {
405
406                                        double dist = Calc.getDistance(tmpfidx[i].getUnitv(), tmpfidx[j].getUnitv());
407                                        if ( dist  < ddiff) {
408
409                                                // centroids have to be closer than fragCompatDist
410                                                double dd = Calc.getDistance(tmpfidx[i].getCenter1(),tmpfidx[j].getCenter1()) -
411                                                                  Calc.getDistance(tmpfidx[i].getCenter2(),tmpfidx[j].getCenter2());
412                                                if ( dd < 0)
413                                                        dd = -dd;
414                                                if ( dd < fragCompatDist){
415                                                        maxi = p2i+l2;
416                                                        used[j]=1;
417                                                        f.add(p2i,p2j,0,tmpfidx[j].getLength());
418                                                }
419                                        }
420
421
422                                }
423                                j+=1;
424                        }
425
426                        int red = 0;
427                        for (int k = 0 ; k < n ; k ++) {
428                                if (used[k] == 0) {
429                                        tmpfidx[red] = tmpfidx[k];
430                                        red+=1;
431                                }
432                        }
433
434
435                        used = new int[n];
436                        tmpfidx = (FragmentPair[])resizeArray(tmpfidx,red);
437
438                        fll.add(f);
439                }
440
441
442                Comparator<JointFragments> comp = new JointFragmentsComparator();
443                Collections.sort(fll,comp);
444                Collections.reverse(fll);
445                int m = Math.min(maxRefine,fll.size());
446                List<JointFragments> retlst = new ArrayList<JointFragments>();
447                for ( int i = 0 ; i < m ; i++){
448                        JointFragments jf = fll.get(i);
449                        retlst.add(jf);
450                }
451
452                return retlst.toArray(new JointFragments[retlst.size()]);
453
454        }
455
456
457        public void extendFragments(Atom[] ca1, Atom[] ca2 ,JointFragments[] fragments, StrucAligParameters params) throws StructureException {
458
459                for(JointFragments p : fragments){
460                        extendFragments(ca1, ca2, p, params);
461                }
462
463        }
464
465        public void extendFragments(Atom[] ca1, Atom[] ca2 , JointFragments fragments, StrucAligParameters params) throws StructureException {
466
467                List<int[]> pos = fragments.getIdxlist();
468
469                int[] firstP = pos.get(0);
470                int pstart1 = firstP[0];
471                int pstart2 = firstP[1];
472
473                int[] lastP = pos.get(pos.size()-1);
474                int pend1 = lastP[0];
475                int pend2 = lastP[1];
476
477                boolean startReached = false;
478                boolean endReached   = false;
479
480                while (! (startReached && endReached)){
481
482                        if ( ! (startReached) && ((pstart1 <= 0) || (pstart2 <= 0))) {
483                                startReached = true;
484
485                        } else {
486                                pstart1--;
487                                pstart2--;
488                        }
489
490                        if ( ! (endReached) && (( pend1 >= (ca1.length -1) ) || ( pend2 >= ca2.length -1  )) ){
491                                endReached = true;
492                        } else {
493                                pend1++;
494                                pend2++;
495                        }
496
497                        if ( ! startReached){
498                                double newRms1 = testAdd(ca1, ca2, fragments,pstart1,pstart2);
499
500                                if (newRms1 < 3.7 ) {
501                                        fragments.add(pstart1,pstart2,0,1);
502                                } else {
503                                        startReached = true;
504                                }
505                        }
506
507                        if( ! endReached){
508
509                                double newRms2 = testAdd(ca1, ca2, fragments, pend1, pend2);
510                                if ( newRms2 < 3.7) {
511                                        fragments.add(pend1,pend2,0,1);
512                                } else {
513                                        endReached = true;
514                                }
515                        }
516
517                }
518
519        }
520
521
522        private double testAdd(Atom[] ca1, Atom[] ca2, JointFragments fragments, int pstart1, int pstart2) throws StructureException {
523
524                JointFragments frag = new JointFragments();
525                frag.setIdxlist(fragments.getIdxlist());
526                frag.add(pstart1, pstart2, 0,1);
527
528                return FragmentJoiner.getRMS(ca1, ca2, frag);
529
530        }
531
532}
533
534
535
536
537class JointFragmentsComparator
538                  implements Comparator<JointFragments> {
539
540        @Override
541        public int compare(JointFragments one, JointFragments two) {
542
543
544                int s1 = one.getIdxlist().size();
545                int s2 = two.getIdxlist().size();
546
547                double rms1 = one.getRms();
548                double rms2 = two.getRms();
549                if ( s1 > s2 ) {
550                        return 1 ;
551                }
552
553                else if ( s1 < s2){
554                        return -1;
555                }
556                else{
557                        if ( rms1 < rms2)
558                                return 1;
559                        if ( rms1 > rms2)
560                                return -1;
561                        return 0;
562                }
563        }
564
565
566}