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 */
021package org.biojava.nbio.structure.domain.pdp;
022
023import org.biojava.nbio.structure.Atom;
024
025import java.util.List;
026
027public class Cut {
028
029        static boolean verbose = CutDomain.verbose;
030
031        public int cut( Atom[] ca, Domain dom, CutValues val, int[][] dist, PDPDistanceMatrix pdpMatrix) {
032
033                int nclose = pdpMatrix.getNclose();
034
035                int[] iclose = pdpMatrix.getIclose();
036                int[] jclose = pdpMatrix.getJclose();
037
038                int[] contacts = new int[PDPParameters.MAXLEN];
039                double[] max_contacts = new double [PDPParameters.MAXLEN];
040                double[] contact_density = new double[PDPParameters.MAXLEN];
041                double average_density,x,y;
042
043                int endsf,endst;
044                int k,l,nc;
045                int size1,size2;
046                int size1t,size2t;
047                int size11,size22,size0;
048                int contactsd;
049                int iseg,jseg,kseg;
050                int from,to,from1,to1,from2,to2,lseg;
051
052
053                int site_min = -1;
054
055
056                // AP add sort here..
057                //qsort(dom.segment,dom.nseg,sizeof(struct Segment),segcmp);
058                // what is going on with the segments??
059
060                List<Segment> segments = dom.getSegments();
061
062                java.util.Collections.sort(segments, new SegmentComparator());
063
064                if ( verbose)
065                        System.out.println("  ---  Cut.cut " + dom + " ");
066                average_density = 0.0;
067                size0=0;
068                for(iseg=0;iseg<dom.nseg;iseg++) {
069                        contactsd=1;
070                        size1t=0;
071                        size2t=0;
072                        for(jseg=0;jseg<iseg;jseg++)
073                                size1t+=(dom.getSegmentAtPos(jseg).getFrom() - dom.getSegmentAtPos(jseg).getFrom() + 1);
074                        for(jseg=iseg+1;jseg<dom.nseg;jseg++)
075                                size2t+=(dom.getSegmentAtPos(jseg).getTo() - dom.getSegmentAtPos(jseg).getFrom() + 1);
076                        for(jseg=0;jseg<iseg;jseg++) {
077                                from1 = dom.getSegmentAtPos(jseg).getFrom();
078                                to1 = dom.getSegmentAtPos(jseg).getTo();
079                                for(int i=from1;i<to1;i++) {
080                                        for(kseg=iseg+1;kseg<dom.nseg;kseg++) {
081                                                from2 = dom.getSegmentAtPos(kseg).getFrom();
082                                                to2 = dom.getSegmentAtPos(kseg).getFrom();
083                                                for(int j=from2;j<to2;j++)
084                                                        if(Math.abs(i-j)>4) contactsd+=(dist[i][j]);
085                                        }
086                                }
087                        }
088                        from = dom.getSegmentAtPos(iseg).getFrom();
089                        to = dom.getSegmentAtPos(iseg).getTo();
090                        for(k=from;k<to;k++) {
091                                contacts[k] = contactsd;
092                                /*
093        if(k==392) printf("init contacts = %d\n",contacts[k]);
094                                 */
095                                size11=size1t+(k-from+1);
096                                size22=size2t+(to-k);
097                                for(int i=from;i<=k;i++) {
098                                        for(kseg=iseg+1;kseg<dom.nseg;kseg++) {
099                                                from2 = dom.getSegmentAtPos(kseg).getFrom();
100                                                to2 = dom.getSegmentAtPos(kseg).getTo();
101                                                for(int j=from2;j<=to2;j++)
102                                                        if(Math.abs(i-j)>4) contacts[k]+=(dist[i][j]);
103                                        }
104                                }
105                                /*
106        if(k==392) printf("[from,k]x]iseg,nseg[ = %d\n",contacts[k]);
107                                 */
108                                for(int i=from;i<=k;i++) {
109                                        for (int j=k+1;j<=to;j++)
110                                                if(Math.abs(i-j)>4) contacts[k]+=(dist[i][j]);
111                                }
112                                /*
113        if(k==392) printf("[from,k]x]k,to[ = %d\n",contacts[k]);
114                                 */
115                                for(int i=k+1;i<=to;i++) {
116                                        for(kseg=0;kseg<iseg;kseg++) {
117                                                from2 = dom.getSegmentAtPos(kseg).getFrom();
118                                                to2 = dom.getSegmentAtPos(kseg).getTo();
119                                                for(int j=from2;j<to2;j++)
120                                                        if(Math.abs(i-j)>4) contacts[k]+=(dist[j][i]);
121                                        }
122                                }
123                                /*
124        if(k==392) printf("]k,to]x]0,iseg[ = %d\n",contacts[k]);
125                                 */
126                                size1=Math.min(size11,size22);
127                                size2=Math.max(size11,size22);
128                                x=Math.min(PDPParameters.MAXSIZE,size1);
129                                y=Math.min(PDPParameters.MAXSIZE,size2);
130                                if(x>150&&y>1.5*x) y=1.5*x;
131                                else if(y>2*x) y=2*x;
132                                /*
133                                 */
134                                x=Math.min(Math.pow(x,1.3/3)+PDPParameters.RG,Math.pow(x,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG);
135                                y=Math.min(Math.pow(y,1.3/3)+PDPParameters.RG,Math.pow(y,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG);
136                                /* max_ contacts depend on the size of domains */
137                                /* stella wanted comments at this point */
138                                /*
139                        max_contacts[k] = min(MAXCONT,x*y);
140                                 */
141                                max_contacts[k] = 10*x*y;
142                                if(size1>150) max_contacts[k] = 9*x*y;
143                                contact_density[k]=contacts[k]/max_contacts[k];
144                                /*
145                        if(contact_density[k]>2.5) contact_density[k]=2.5;
146                                 */
147                                /*
148        if(first_cut)
149                if(k==277 && !first_cut)
150                printf("k=%d s1=%d s2=%d x=%f y=%f mc=%d c=%d cd=%f\n",k,size1,size2,x,y,max_contacts[k],contacts[k],contact_density[k]);
151                                 */
152                                //if(verbose) System.out.println(String.format("%d      %d      %d      %f      %f      %f      %d      %f",k,size1,size2,x,y,max_contacts[k],contacts[k],contact_density[k]));
153                                /*
154                                 */
155                                if(from==0) endsf = PDPParameters.ENDSEND;
156                                else endsf = PDPParameters.ENDS;
157                                if(to==ca.length-1) endst = PDPParameters.ENDSEND;
158                                else endst = PDPParameters.ENDS;
159                                if((contact_density[k])<val.s_min&&k>from+endsf&&k<to-endst) {
160                                        val.s_min = (contact_density[k]);
161                                        site_min=k+1;
162                                }
163                                if(k>from+endsf&&k<to-endst) {
164                                        average_density+=contact_density[k];
165                                        size0++;
166                                }
167                        }
168                }
169                average_density/=size0;
170                if(verbose) System.out.println(String.format("  --- Trying to cut domain of size %d having %d segments and  average cont_density %f\n",dom.size,dom.nseg,average_density));
171
172                if ( verbose )
173                        for(kseg=0;kseg<dom.nseg;kseg++)
174                                System.out.println(String.format("      --- segment %d from %d to %d",kseg,dom.getSegmentAtPos(kseg).getFrom(),dom.getSegmentAtPos(kseg).getTo()) + " av density: " + average_density);
175
176
177                if(val.first_cut) {
178                        /*
179                for(k=from+ENDS;k<to-ENDS;k++)
180                        printf("%d      %s      %d      %d      %d      %d      %f\n",k,protein.res[k].type,k-from,to-k,contacts[k],max_contacts[k],contact_density[k]);
181                         */
182                        val.AD = average_density;
183                        if(verbose)     System.out.println(String.format("  --- AD=%f", average_density));
184                        /*
185                         */
186                }
187                val.AD = average_density;
188
189                val.s_min/=val.AD;
190
191                if(verbose) System.out.println(String.format("  --- after single cut: s_min = %f site_min = %d",val.s_min,site_min));
192
193                ///
194                k=0;
195
196                /* check double cuts */
197                if ( verbose )
198                System.out.println("  --- checking double cuts up to: " + nclose);
199                nc=0;
200                for(l=0;l<nclose;l++) {
201
202                        /************ find iseg, jseg ****************/
203                        iseg=jseg=-1;
204                        for(kseg=0;kseg<dom.nseg;kseg++) {
205                                from=dom.getSegmentAtPos(kseg).getFrom();
206                                to=dom.getSegmentAtPos(kseg).getTo();
207                                if(from==0) endsf = PDPParameters.ENDSEND;
208                                else endsf = PDPParameters.ENDS;
209                                if(to==ca.length-1) endst = PDPParameters.ENDSEND;
210                                else endst = PDPParameters.ENDS;
211                                if(iclose[l]>from+endsf&&iclose[l]<to-endst)
212                                        iseg=kseg;
213                                if(jclose[l]>from+endsf&&jclose[l]<to-endst)
214                                        jseg=kseg;
215                        }
216                        if(iseg<0||jseg<0) continue;
217                        /*
218        printf("l = %d iclose[l] = %d jclose[l] = %d\n",l,iclose[l],jclose[l]);
219                         */
220                        /*********************************************/
221
222                        from=dom.getSegmentAtPos(iseg).getFrom();
223                        to=dom.getSegmentAtPos(iseg).getTo();
224                        from1=dom.getSegmentAtPos(jseg).getFrom();
225                        to1=dom.getSegmentAtPos(jseg).getTo();
226
227                        /************ count contacts *****************/
228                        contacts[nc] = 1;
229                        //no=0;
230
231                        /******* contacts between [0,iseg[ and ]iseg,jseg[ ********/
232                        for(kseg=0;kseg<iseg;kseg++)
233                                for(lseg=iseg+1;lseg<jseg;lseg++)
234                                        for( int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++)
235                                                for(int j=dom.getSegmentAtPos(lseg).getFrom();j<dom.getSegmentAtPos(lseg).getTo();j++) {
236                                                        contacts[nc]+=(dist[i][j]);
237                                                }
238
239                        //System.out.println(String.format("[0,iseg[ - ]iseg,jseg[ : %d\n",contacts[nc]-no));
240
241                        //no=contacts[nc];
242                        /**********************************************************/
243                        /**********************************************************/
244                        /******* contacts between ]jseg,nseg[ and ]iseg,jseg[ ********/
245                        for(kseg=jseg+1;kseg<dom.nseg;kseg++)
246                                for(lseg=iseg+1;lseg<jseg;lseg++)
247                                        for(int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++)
248                                                for(int j=dom.getSegmentAtPos(lseg).getFrom();j<dom.getSegmentAtPos(lseg).getTo();j++) {
249                                                        contacts[nc]+=(dist[j][i]);
250                                                }
251                        /*
252                printf("]jseg,nseg] - ]iseg,jseg[ : %d\n",contacts[nc]-no);
253                         */
254                        //no=contacts[nc];
255                        /*************************************************************/
256                        /*************************************************************/
257                        /**** contacts between [from,iclose] in iseg and ]iseg,jseg[ ****/
258                        if(iseg==jseg) {
259                                //System.out.println(" CONTACT:  " + from + " " + iclose[l] + " " + iseg + " " + jseg);
260                                for(int i=from;i<=iclose[l];i++) {
261                                        for (int j=iclose[l]+1;j<=jclose[l];j++) {
262                                                contacts[nc]+=(dist[i][j]);
263                                        }
264                                }
265                                for (int j=iclose[l]+1;j<jclose[l];j++) {
266                                        for(kseg=0;kseg<iseg;kseg++)
267                                                for(int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++) {
268                                                        contacts[nc]+=(dist[i][j]);
269                                                }
270                                        for(int i=jclose[l];i<to;i++) {
271                                                contacts[nc]+=(dist[j][i]);
272                                        }
273                                        for(kseg=iseg+1;kseg<dom.nseg;kseg++)
274                                                for(int i=dom.getSegmentAtPos(kseg).getFrom();i<dom.getSegmentAtPos(kseg).getTo();i++) {
275                                                        contacts[nc]+=(dist[j][i]);
276                                                }
277                                }
278                                /*
279                printf("iclose==jclose : %d\n",contacts[nc]-no);
280                                 */
281                                //no=contacts[nc];
282                        }
283                        else {
284                                //System.out.println(" ISEG!=JSEG " + " " + from + " " + iclose[l]);
285                                for(int i=from;i<=iclose[l];i++) {
286                                        for(kseg=iseg+1;kseg<jseg;kseg++)
287                                                for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) {
288                                                        contacts[nc]+=(dist[i][j]);
289                                                }
290                                        for(int j=from1;j<jclose[l];j++) {
291                                                contacts[nc]+=(dist[i][j]);
292                                        }
293                                        for(int j=iclose[l]+1;j<to;j++) {
294                                                contacts[nc]+=(dist[i][j]);
295                                        }
296                                }
297                                for(int i=iclose[l]+1;i<to;i++) {
298                                        for(kseg=0;kseg<iseg;kseg++)
299                                                for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) {
300                                                        contacts[nc]+=(dist[j][i]);
301                                                }
302                                        for(kseg=jseg+1;kseg<dom.nseg;kseg++)
303                                                for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) {
304                                                        contacts[nc]+=(dist[i][j]);
305                                                }
306                                        for(int j=jclose[l];j<=to1;j++) {
307                                                contacts[nc]+=(dist[i][j]);
308                                        }
309                                }
310                                for (int i=from1;i<jclose[l];i++) {
311                                        for(kseg=0;kseg<iseg;kseg++)
312                                                for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) {
313                                                        contacts[nc]+=(dist[j][i]);
314                                                }
315                                        for(kseg=jseg+1;kseg<dom.nseg;kseg++)  {
316                                                for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++)
317                                                        contacts[nc]+=(dist[i][j]);
318                                        }
319                                        for(int j=jclose[l];j<to1;j++) {
320                                                contacts[nc]+=(dist[i][j]);
321                                        }
322                                }
323                                for(int i=jclose[l];i<to1;i++)
324                                        for(kseg=iseg+1;kseg<jseg;kseg++)
325                                                for(int j=dom.getSegmentAtPos(kseg).getFrom();j<dom.getSegmentAtPos(kseg).getTo();j++) {
326                                                        /*
327                                                if(iclose[l]==33&&jclose[l]==69&&dist[i][j]) printf("%d %s %d %s %d\n",i,protein.res[i].type,j,protein.res[j].type,dist[i][j]);
328                                                         */
329                                                        contacts[nc]+=(dist[j][i]);
330                                                }
331                        }
332                        /*******************************************************************/
333                        /*******************************************************************/
334                        /*******************************************************************/
335                        size11=0;
336                        size22=0;
337                        for(kseg=0;kseg<iseg;kseg++)
338                                size11+=(dom.getSegmentAtPos(kseg).getTo()-dom.getSegmentAtPos(kseg).getFrom()+1);
339                        for(kseg=jseg+1;kseg<dom.nseg;kseg++)
340                                size11+=(dom.getSegmentAtPos(kseg).getTo()-dom.getSegmentAtPos(kseg).getFrom()+1);
341                        size11+=(iclose[l]-from+1);
342                        size11+=(to1-jclose[l]+1);
343                        /*
344        printf("size11 = %d from = %d to1 = %d \n",size11,from,to1);
345                         */
346                        for(kseg=iseg+1;kseg<jseg;kseg++)
347                                size22+=(dom.getSegmentAtPos(kseg).getTo()-dom.getSegmentAtPos(kseg).getFrom()+1);
348                        if(iseg==jseg)
349                                size22+=(jclose[l]-iclose[l]);
350                        else {
351                                size22+=(jclose[l]-from1);
352                                size22+=(to-iclose[l]);
353                        }
354                        /*
355        printf("size22 = %d from1 = %d jclose %d iclose %d \n",size22,from1,jclose[l],iclose[l]);
356                         */
357                        size1=Math.min(size11,size22);
358                        size2=Math.max(size11,size22);
359                        x=Math.min(PDPParameters.MAXSIZE,size1);
360                        y=Math.max(PDPParameters.MAXSIZE,size2);
361                        if(y>2*x) y=2*x;
362                        /*
363                         */
364                        x=Math.min(Math.pow(x,1.3/3)+PDPParameters.RG,Math.pow(x,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG);
365                        y=Math.min(Math.pow(y,1.3/3)+PDPParameters.RG,Math.pow(y,1.1/3)+Math.pow(PDPParameters.TD,1.3/3)+PDPParameters.RG);
366                        /*
367                max_contacts[nc] = min(MAXCONT,x*y);
368                         */
369                        max_contacts[nc] = x*y*10;
370                        if(size1>150) max_contacts[k] = 9*x*y;
371                        contact_density[nc]=contacts[nc]/max_contacts[nc];
372                        /*
373        if(first_cut)
374                         */
375                        //if ( verbose)
376                        //      System.out.println("double cut" + l + " " + ca[iclose[l]].getGroup().getType() + " " + iclose[l] + " " + jclose[l] + " " + contacts[nc]+
377                        //                      " " + max_contacts[nc] + " " + x+ " " + y + " " + size11 + " " + size22 + " " + contact_density[nc] + " " + contact_density[nc]/val.AD );
378                        //if(verbose) System.out.println(String.format(" double cut: %d %s %d %d c=%d mc=%d x=%f y=%f s1=%d s2=%d cd=%f cd/ad=%f\n",l,ca[iclose[l]].getGroup().getType(),iclose[l],jclose[l],contacts[nc],max_contacts[nc],x,y,size11,size22,contact_density[nc],contact_density[nc]/val.AD));
379                        if((contact_density[nc]/val.AD+PDPParameters.DBL)<val.s_min&&contact_density[nc]/val.AD+PDPParameters.DBL<PDPParameters.CUT_OFF_VALUE2) {
380                                /*
381        printf("________________\n");
382                                 */
383                                val.s_min = (contact_density[nc]/val.AD)+PDPParameters.DBL;
384                                site_min=iclose[l];
385                                val.site2=jclose[l];
386                        }
387
388                        nc++;
389                        if ( nc >= PDPParameters.MAXSIZE)
390                                nc = PDPParameters.MAXSIZE-1;
391                }
392                val.first_cut=false;
393                if(verbose)
394                        System.out.println(String.format("  --- E ... at the end of cut: s_min %f CUTOFF %f site_min %d *site2 %d",val.s_min,PDPParameters.CUT_OFF_VALUE,site_min,val.site2));
395                if(val.s_min> PDPParameters.CUT_OFF_VALUE) return -1;
396
397                return(site_min);
398        }
399}