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.*;
024
025
026public class GetDistanceMatrix {
027
028
029        /** A set of Calpha atoms that are representing the protein
030         *
031         * @param protein
032         */
033        public  PDPDistanceMatrix getDistanceMatrix(Atom[] protein) throws StructureException{
034                int[][] dist = new int[protein.length+3][protein.length+3];
035                int i,j;
036                double d,dt1,dt2,dt3,dt4;
037                int nclose=0;
038                int[] iclose = new int[protein.length*protein.length];
039                int[] jclose= new int[protein.length*protein.length];
040
041                if(protein.length >= PDPParameters.MAXLEN) {
042                        System.err.println(String.format("%d protein.len > MAXLEN %d\n",protein.length,PDPParameters.MAXLEN));
043                        return null;
044                }
045                for(i=0;i<protein.length;i++) {
046                        for(j=i;j<protein.length;j++) {
047                                dist[i][j]=0;
048                                dist[j][i]=0;
049
050                                d=0;
051
052                                Atom ca1 = protein[i];
053                                Atom ca2 = protein[j];
054                                Group g1 = ca1.getGroup();
055                                Group g2 = ca2.getGroup();
056
057                                Atom cb1 = getCBeta(g1);
058                                Atom cb2 = getCBeta(g2);
059                                boolean hasCbeta1 = cb1 != null;
060                                boolean hasCbeta2 = cb2 != null;
061
062                                dt1=81;
063                                dt2=64;
064                                dt3=49;
065                                dt4=36;
066
067                                if(hasCbeta1 && hasCbeta2) {
068                                        double distance = Calc.getDistance(cb1,cb2);
069                                        d+= distance*distance;
070                                }
071                                else if(hasCbeta1 && ! hasCbeta2) {
072                                        double distance = 999;
073
074                                        distance = Calc.getDistance(cb1, ca2);
075                                        d += distance * distance;
076                                }
077                                else if(!hasCbeta1&&hasCbeta2) {
078                                        double distance = Calc.getDistance(ca1, cb2);
079                                        d += distance * distance;
080                                }
081                                else if( ! hasCbeta1&&!hasCbeta2) {
082                                        double distance = Calc.getDistance(ca1, ca2);
083                                        d += distance * distance;
084                                }
085
086                                if(d<dt1) {
087                                        dist[i][j]=1;
088                                        dist[j][i]=1;
089                                        if(d<dt2) {
090                                                dist[i][j]=2;
091                                                dist[j][i]=2;
092                                                if(j-i>35) {
093                                                        iclose[nclose]=i;
094                                                        jclose[nclose]=j;
095                                                        nclose++;
096                                                }
097                                                if(d<dt3) {
098                                                        dist[i][j]=4;
099                                                        dist[j][i]=4;
100                                                        if(d<dt4) {
101                                                                dist[i][j]=6;
102                                                                dist[j][i]=6;
103                                                        }
104                                                }
105                                        }
106                                }
107                        }
108                }
109                /* secondary structure interaction */
110                for(i=1;i<protein.length;i++) {
111                        for(j=i;j<protein.length-1;j++) {
112                                /* beta-sheet */
113                                if(dist[i][j]>=2&&j-i>5) {
114                                        if(dist[i-1][j-1]>=2&&dist[i+1][j+1]>=2||dist[i-1][j+1]>=2&&dist[i+1][j-1]>=2) {
115                                                dist[i][j]+=4;
116                                                dist[j][i]+=4;
117                                                /*
118                                        printf("1: %d %d %d\n",i,j,dist[i][j]);
119                                                 */
120                                        }
121                                        /* alpha-helices */
122                                        else if(i>2&&j<protein.length-2) {
123                                                if(dist[i-3][j-3]>=1&&dist[i+3][j+3]>=1||dist[i-3][j+3]>=1&&dist[i+3][j-3]>=1) {
124                                                        dist[i][j]+=4;
125                                                        dist[j][i]+=4;
126                                                        /*
127                                                printf("3: %d %d %d\n",i,j,dist[i][j]);
128                                                         */
129                                                }
130                                                else if(i>3&&j<protein.length-3) {
131                                                        if((dist[i-3][j-3]>=1||dist[i-3][j-4]>=1||dist[i-4][j-3]>=1||dist[i-4][j-4]>=1)&&
132                                                                        (dist[i+4][j+4]>=1||dist[i+4][j+3]>=1||dist[i+3][j+3]>=1||dist[i+3][j+4]>=1)
133                                                                        ||(dist[i-4][j+4]>=1||dist[i-4][j+3]>=1||dist[i-3][j+4]>=1||dist[i-3][j+3]>=1)&&
134                                                                        (dist[i+4][j-4]>=1||dist[i+4][j-3]>=1||dist[i+3][j-4]>=1||dist[i+3][j-3]>=1)) {
135                                                                dist[i][j]+=4;
136                                                                dist[j][i]+=4;
137                                                                /*
138                                                        printf("4: %d %d %d\n",i,j,dist[i][j]);
139                                                                 */
140                                                        }
141                                                }
142                                        }
143                                }
144                        }
145                }
146
147                PDPDistanceMatrix matrix = new PDPDistanceMatrix();
148
149                matrix.setNclose(nclose);
150                matrix.setIclose(iclose);
151                matrix.setJclose(jclose);
152                matrix.setDist(dist);
153                return matrix;
154
155        }
156
157
158
159        private Atom getCBeta(Group g1) {
160                Atom cb = null;
161
162
163                cb = g1.getAtom("CB");
164                if ( cb == null)
165                        {
166                        if ( g1 instanceof AminoAcid) {
167                                AminoAcid aa = (AminoAcid) g1;
168
169                                try {
170                                        cb = Calc.createVirtualCBAtom(aa);
171                                } catch (StructureException e1) {
172                                        // TODO Auto-generated catch block
173                                        e1.printStackTrace();
174                                }
175                        }
176                }
177                return cb;
178        }
179
180
181
182
183
184
185
186
187
188}