src/jalview/analysis/AlignSeq.java

   1 /* Jalview - a java multiple alignment editor\r
   2  * Copyright (C) 1998  Michele Clamp\r
   3  *\r
   4  * This program is free software; you can redistribute it and/or\r
   5  * modify it under the terms of the GNU General Public License\r
   6  * as published by the Free Software Foundation; either version 2\r
   7  * of the License, or (at your option) any later version.\r
   8  *\r
   9  * This program is distributed in the hope that it will be useful,\r
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  12  * GNU General Public License for more details.\r
  13  *\r
  14  * You should have received a copy of the GNU General Public License\r
  15  * along with this program; if not, write to the Free Software\r
  16  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.\r
  17  */\r
  18 package jalview.analysis;\r
  19 \r
  20 import jalview.schemes.*;\r
  21 import jalview.datamodel.SequenceI;\r
  22 import jalview.util.*;\r
  23 import jalview.io.*;\r
  24 \r
  25 import java.util.*;\r
  26 import java.io.*;\r
  27 import java.awt.*;\r
  28 \r
  29 public class AlignSeq {\r
  30 \r
  31   int[][] score;\r
  32   int[][] E;\r
  33   int[][] F;\r
  34   int[][] traceback;\r
  35 \r
  36   int[] seq1;\r
  37   int[] seq2;\r
  38 \r
  39   SequenceI s1;\r
  40   SequenceI s2;\r
  41 \r
  42   String s1str;\r
  43   String s2str;\r
  44 \r
  45   int maxi;\r
  46   int maxj;\r
  47 \r
  48   int[] aseq1;\r
  49   int[] aseq2;\r
  50 \r
  51   String astr1 = "";\r
  52   String astr2 = "";\r
  53 \r
  54   public int seq1start;\r
  55   public int seq1end;\r
  56   public int seq2start;\r
  57   public int seq2end;\r
  58 \r
  59   int count;\r
  60 \r
  61   public int maxscore;\r
  62   float pid;\r
  63   int prev = 0;\r
  64 \r
  65   public static java.util.Hashtable dnaHash = new java.util.Hashtable();\r
  66 \r
  67   static  {\r
  68     dnaHash.put("C", new Integer(0));\r
  69     dnaHash.put("T", new Integer(1));\r
  70     dnaHash.put("A", new Integer(2));\r
  71     dnaHash.put("G", new Integer(3));\r
  72     dnaHash.put("-", new Integer(4));\r
  73   }\r
  74 \r
  75   static String dna[] = {"C","T","A","G","-"};\r
  76   static String pep[] = {"A","R","N","D","C","Q","E","G","H","I","L","K","M","F","P","S","T","W","Y","V","B","Z","X","-"};\r
  77 \r
  78   int gapOpen = 120;\r
  79   int gapExtend = 20;\r
  80 \r
  81   int lookup[][] = ResidueProperties.getBLOSUM62();\r
  82   String intToStr[] = pep;\r
  83   int defInt = 23;\r
  84 \r
  85   String output = "";\r
  86 \r
  87   String type;\r
  88   Runtime rt;\r
  89   public AlignSeq() {}\r
  90 \r
  91   public AlignSeq(SequenceI s1, SequenceI s2,String type) {\r
  92     rt  = Runtime.getRuntime();\r
  93     SeqInit(s1,s2,type);\r
  94   }\r
  95 \r
  96   public int getMaxScore() {\r
  97     return maxscore;\r
  98   }\r
  99 \r
 100   public int getSeq2Start() {\r
 101     return seq2start;\r
 102   }\r
 103 \r
 104   public int getSeq2End() {\r
 105     return seq2end;\r
 106   }\r
 107 \r
 108   public int getSeq1Start() {\r
 109     return seq1start;\r
 110   }\r
 111 \r
 112   public int getSeq1End() {\r
 113     return seq1end;\r
 114   }\r
 115 \r
 116   public String getOutput() {\r
 117     return output;\r
 118   }\r
 119 \r
 120   public String getAStr1() {\r
 121     return astr1;\r
 122   }\r
 123   public String getAStr2() {\r
 124     return astr2;\r
 125   }\r
 126   public int [] getASeq1() {\r
 127     return aseq1;\r
 128   }\r
 129   public int [] getASeq2() {\r
 130     return aseq2;\r
 131   }\r
 132   public SequenceI getS1() {\r
 133     return s1;\r
 134   }\r
 135   public SequenceI getS2() {\r
 136     return s2;\r
 137   }\r
 138 \r
 139   public void SeqInit(SequenceI s1, SequenceI s2,String type) {\r
 140     s1str = extractGaps(".",s1.getSequence());\r
 141     s2str = extractGaps(".",s2.getSequence());\r
 142     s1str = extractGaps("-",s1str);\r
 143     s2str = extractGaps("-",s2str);\r
 144     s1str = extractGaps(" ",s1str);\r
 145     s2str = extractGaps(" ",s2str);\r
 146 \r
 147     this.s1 = s1;\r
 148     this.s2 = s2;\r
 149 \r
 150     this.type = type;\r
 151 \r
 152     if (type.equals("pep")) {\r
 153       lookup = ResidueProperties.getBLOSUM62();\r
 154       intToStr = pep;\r
 155       defInt = 23;\r
 156     } else if (type.equals("dna")) {\r
 157       lookup = ResidueProperties.getDNA();\r
 158       intToStr = dna;\r
 159       defInt = 4;\r
 160     } else {\r
 161       output = output + ("Wrong type = dna or pep only");\r
 162       System.exit(0);\r
 163     }\r
 164 \r
 165 \r
 166     //System.out.println("lookuip " + rt.freeMemory() + " "+  rt.totalMemory());\r
 167     seq1 = new int[s1str.length()];\r
 168     //System.out.println("seq1 " + rt.freeMemory() +" "  + rt.totalMemory());\r
 169     seq2 = new int[s2str.length()];\r
 170     //System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());\r
 171     score = new int[s1str.length()][s2str.length()];\r
 172     //System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());\r
 173     E = new int[s1str.length()][s2str.length()];\r
 174     //System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());\r
 175     F = new int[s1str.length()][s2str.length()];\r
 176     traceback = new int[s1str.length()][s2str.length()];\r
 177     //System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());\r
 178     seq1 = stringToInt(s1str,type);\r
 179     //System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());\r
 180     seq2 = stringToInt(s2str,type);\r
 181     //System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());\r
 182 \r
 183     //   long tstart = System.currentTimeMillis();\r
 184     //    calcScoreMatrix();\r
 185     //long tend = System.currentTimeMillis();\r
 186 \r
 187     //System.out.println("Time take to calculate score matrix = " + (tend-tstart) + " ms");\r
 188 \r
 189 \r
 190     //   printScoreMatrix(score);\r
 191     //System.out.println();\r
 192 \r
 193     //printScoreMatrix(traceback);\r
 194     //System.out.println();\r
 195 \r
 196     //  printScoreMatrix(E);\r
 197     //System.out.println();\r
 198 \r
 199     ///printScoreMatrix(F);\r
 200     //System.out.println();\r
 201     // tstart = System.currentTimeMillis();\r
 202     //traceAlignment();\r
 203     //tend = System.currentTimeMillis();\r
 204     //System.out.println("Time take to traceback alignment = " + (tend-tstart) + " ms");\r
 205   }\r
 206 \r
 207   public void traceAlignment() {\r
 208 \r
 209     // Find the maximum score along the rhs or bottom row\r
 210     int max = -9999;\r
 211     for (int i = 0; i < seq1.length; i++) {\r
 212       if (score[i][seq2.length - 1] > max ) {\r
 213         max = score[i][seq2.length - 1];\r
 214         maxi = i;\r
 215         maxj = seq2.length-1;\r
 216       }\r
 217     }\r
 218     for (int j = 0; j < seq2.length; j++) {\r
 219       if (score[seq1.length - 1][j] > max ) {\r
 220         max = score[seq1.length - 1][j];\r
 221         maxi = seq1.length-1;\r
 222         maxj = j;\r
 223       }\r
 224     }\r
 225 \r
 226     //  System.out.println(maxi + " " + maxj + " " + score[maxi][maxj]);\r
 227 \r
 228     int i = maxi;\r
 229     int j = maxj;\r
 230     int trace;\r
 231     maxscore = score[i][j] / 10;\r
 232 \r
 233     seq1end = maxi+1;\r
 234     seq2end = maxj+1;\r
 235 \r
 236     aseq1 = new int[seq1.length + seq2.length];\r
 237     aseq2 = new int[seq1.length + seq2.length];\r
 238 \r
 239     count = seq1.length + seq2.length - 1;\r
 240 \r
 241     while (i>0 && j >0) {\r
 242 \r
 243       if (aseq1[count] != defInt && i >=0) {\r
 244         aseq1[count] = seq1[i];\r
 245         astr1 = intToStr[seq1[i]] + astr1;\r
 246       }\r
 247 \r
 248       if (aseq2[count] != defInt && j > 0) {\r
 249         aseq2[count] = seq2[j];\r
 250         astr2 = intToStr[seq2[j]] + astr2;\r
 251       }\r
 252       trace = findTrace(i,j);\r
 253       if (trace == 0) {\r
 254         i--;\r
 255         j--;\r
 256 \r
 257       } else  if (trace == 1) {\r
 258         j--;\r
 259         aseq1[count] = defInt;\r
 260         astr1 = "-" + astr1.substring(1);\r
 261       } else  if (trace == -1) {\r
 262         i--;\r
 263         aseq2[count] = defInt;\r
 264         astr2 = "-" + astr2.substring(1);\r
 265       }\r
 266       count--;\r
 267     }\r
 268 \r
 269     seq1start = i+1;\r
 270     seq2start = j+1;\r
 271 \r
 272     if (aseq1[count] != defInt) {\r
 273       aseq1[count] = seq1[i];\r
 274       astr1 = intToStr[seq1[i]] + astr1;\r
 275     }\r
 276 \r
 277     if (aseq2[count] != defInt) {\r
 278       aseq2[count] = seq2[j];\r
 279       astr2 = intToStr[seq2[j]] + astr2;\r
 280     }\r
 281   }\r
 282 \r
 283   public void printAlignment() {\r
 284     // Find the biggest id length for formatting purposes\r
 285     int maxid = s1.getName().length();\r
 286 \r
 287     if (s2.getName().length() > maxid) {\r
 288       maxid = s2.getName().length();\r
 289     }\r
 290 \r
 291     int len = 72 - maxid - 1;\r
 292     int nochunks = ((aseq1.length - count) / len) + 1;\r
 293     pid = 0;\r
 294     int overlap = 0;\r
 295 \r
 296     output = output + ("Score = " + score[maxi][maxj] + "\n");\r
 297     output = output + ("Length of alignment = " + (aseq1.length-count) + "\n");\r
 298     output = output + ("Sequence ");\r
 299     output = output + (new Format("%" + maxid + "s").form(s1.getName()));\r
 300     output = output + (" :  " + seq1start + " - " + seq1end + " (Sequence length = " + s1str.length() + ")\n");\r
 301     output = output + ("Sequence ");\r
 302     output = output + (new Format("%" + maxid + "s").form(s2.getName()));\r
 303     output = output + (" :  " + seq2start + " - " + seq2end + " (Sequence length = " + s2str.length() + ")\n\n");\r
 304 \r
 305     for (int j = 0; j < nochunks; j++) {\r
 306       // Print the first aligned sequence\r
 307       output = output + (new Format("%" + (maxid) + "s").form(s1.getName()) + " ");\r
 308       for (int i = 0; i < len ; i++) {\r
 309 \r
 310         if ((count + i + j*len) < aseq1.length) {\r
 311           output = output + (new Format("%s").form(intToStr[aseq1[count + i + j*len]]));\r
 312         }\r
 313       }\r
 314 \r
 315       output = output + ("\n");\r
 316       output = output + (new Format("%" + (maxid) + "s").form(" ") + " ");\r
 317       // Print out the matching chars\r
 318       for (int i = 0; i < len ; i++) {\r
 319 \r
 320         if ((count + i + j*len) < aseq1.length) {\r
 321           if ( intToStr[aseq1[count+i+j*len]].equals(intToStr[aseq2[count+i+j*len]]) && !intToStr[aseq1[count+i+j*len]].equals("-")) {\r
 322             pid++;\r
 323             output = output + ("|");\r
 324           } else if (type.equals("pep")) {\r
 325             if (ResidueProperties.getPAM250(intToStr[aseq1[count+i+j*len]],intToStr[aseq2[count+i+j*len]]) > 0) {\r
 326               output  = output + (".");\r
 327             } else {\r
 328               output = output + (" ");\r
 329             }\r
 330           } else {\r
 331             output = output + (" ");\r
 332           }\r
 333 \r
 334         }\r
 335       }\r
 336       // Now print the second aligned sequence\r
 337       output = output + ("\n");\r
 338       output = output + (new Format("%" + (maxid) + "s").form(s2.getName()) + " " );\r
 339       for (int i = 0; i < len ; i++) {\r
 340         if ((count + i + j*len) < aseq1.length) {\r
 341           output = output + (new Format("%s").form(intToStr[aseq2[count + i + j*len]]));\r
 342         }\r
 343       }\r
 344       output = output + ("\n\n");\r
 345     }\r
 346     pid = pid/(float)(aseq1.length-count)*100;\r
 347     output = output + (new Format("Percentage ID = %2.2f\n\n").form(pid));\r
 348 \r
 349   }\r
 350 \r
 351   public void printScoreMatrix(int[][] mat) {\r
 352     int n = seq1.length;\r
 353     int m = seq2.length;\r
 354 \r
 355     for (int i = 0; i < n;i++) {\r
 356       // Print the top sequence\r
 357       if (i == 0) {\r
 358         Format.print(System.out,"%8s",s2str.substring(0,1));\r
 359         for (int jj = 1;jj < m; jj++) {\r
 360           Format.print(System.out,"%5s",s2str.substring(jj,jj+1));\r
 361         }\r
 362         System.out.println();\r
 363       }\r
 364 \r
 365       for (int j = 0;j < m; j++) {\r
 366         if (j == 0) {\r
 367           Format.print(System.out,"%3s",s1str.substring(i,i+1));\r
 368         }\r
 369         Format.print(System.out,"%3d ",mat[i][j]/10);\r
 370       }\r
 371       System.out.println();\r
 372     }\r
 373   }\r
 374 \r
 375   public int findTrace(int i,int j) {\r
 376     int t = 0;\r
 377     int max = score[i-1][j-1] + lookup[seq1[i]][seq2[j]] * 10;\r
 378 \r
 379     if (F[i][j] > max) {\r
 380       max = F[i][j];\r
 381       t = -1;\r
 382     } else if (F[i][j] == max) {\r
 383       if (prev == -1) {\r
 384         max = F[i][j];\r
 385         t = -1;\r
 386       }\r
 387     }\r
 388     if (E[i][j] >= max) {\r
 389       max = E[i][j];\r
 390       t = 1;\r
 391     } else if (E[i][j] == max) {\r
 392       if (prev == 1) {\r
 393         max = E[i][j];\r
 394         t = 1;\r
 395       }\r
 396     }\r
 397     prev = t;\r
 398     return t;\r
 399   }\r
 400 \r
 401   public void calcScoreMatrix() {\r
 402 \r
 403 \r
 404     int n = seq1.length;\r
 405     int m = seq2.length;\r
 406 \r
 407 \r
 408     // top left hand element\r
 409     score[0][0] =  lookup[seq1[0]][seq2[0]] * 10;\r
 410     E[0][0] = -gapExtend;\r
 411     F[0][0] = 0;\r
 412 \r
 413     // Calculate the top row first\r
 414     for (int j=1; j < m; j++) {\r
 415       // What should these values be? 0 maybe\r
 416       E[0][j] = max(score[0][j-1] - gapOpen,E[0][j-1] - gapExtend);\r
 417       F[0][j] = -gapExtend;\r
 418 \r
 419       score[0][j] = max( lookup[seq1[0]][seq2[j]] * 10 ,-gapOpen,-gapExtend);\r
 420 \r
 421       traceback[0][j] = 1;\r
 422     }\r
 423 \r
 424     // Now do the left hand column\r
 425     for (int i=1; i < n; i++) {\r
 426       E[i][0] =  -gapOpen;\r
 427       F[i][0] =  max(score[i-1][0]-gapOpen,F[i-1][0]-gapExtend);\r
 428 \r
 429       score[i][0] = max( lookup[seq1[i]][seq2[0]] * 10 ,E[i][0],F[i][0]);\r
 430       traceback[i][0] = -1;\r
 431     }\r
 432 \r
 433     // Now do all the other rows\r
 434     for (int i = 1; i < n; i++) {\r
 435       for (int j = 1; j < m; j++) {\r
 436 \r
 437         E[i][j] = max(score[i][j-1] - gapOpen,  E[i][j-1] - gapExtend);\r
 438         F[i][j] = max(score[i-1][j] - gapOpen,  F[i-1][j] - gapExtend);\r
 439 \r
 440         score[i][j] = max(score[i-1][j-1] + lookup[seq1[i]][seq2[j]]*10,\r
 441                           E[i][j],\r
 442                           F[i][j]);\r
 443         traceback[i][j] = findTrace(i,j);\r
 444       }\r
 445     }\r
 446 \r
 447   }\r
 448   public static String extractChars(String chars, String seq) {\r
 449     String out = seq;\r
 450     for (int i=0; i < chars.length(); i++) {\r
 451       String gap = chars.substring(i,i+1);\r
 452       out = extractGaps(gap,out);\r
 453     }\r
 454     return out;\r
 455   }\r
 456   public static String extractGaps(String gapChar, String seq) {\r
 457     StringTokenizer str = new StringTokenizer(seq,gapChar);\r
 458     String newString = "";\r
 459 \r
 460     while (str.hasMoreTokens()) {\r
 461       newString  = newString + str.nextToken();\r
 462     }\r
 463     return newString;\r
 464   }\r
 465 \r
 466 \r
 467   public int max(int i1, int i2, int i3) {\r
 468     int max = i1;\r
 469     if (i2 > i1)  {\r
 470       max = i2;\r
 471     }\r
 472     if (i3 > max) {\r
 473       max = i3;\r
 474     }\r
 475     return max;\r
 476   }\r
 477 \r
 478   public int max(int i1, int i2) {\r
 479     int max = i1;\r
 480     if (i2 > i1)  {\r
 481       max = i2;\r
 482     }\r
 483     return max;\r
 484   }\r
 485 \r
 486   public int[] stringToInt(String s,String type) {\r
 487     int[] seq1 = new int[s.length()];\r
 488 \r
 489     for (int i = 0;i < s.length(); i++) {\r
 490       String ss = s.substring(i,i+1).toUpperCase();\r
 491       try {\r
 492         if (type.equals("pep")) {\r
 493           seq1[i] = ((Integer)ResidueProperties.aaHash.get(ss)).intValue();\r
 494         } else if (type.equals("dna")) {\r
 495           seq1[i] = ((Integer)dnaHash.get(ss)).intValue();\r
 496         }\r
 497         if (seq1[i] > 23) {\r
 498           seq1[i] = 23;\r
 499         }\r
 500       } catch (Exception e) {\r
 501         if (type.equals("dna")) {\r
 502           seq1[i] = 4;\r
 503         } else {\r
 504           seq1[i] = 23;\r
 505         }\r
 506       }\r
 507     }\r
 508     return seq1;\r
 509   }\r
 510 \r
 511   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,int psize) {\r
 512 \r
 513     int max = -1000;\r
 514     int min = 1000;\r
 515 \r
 516     for (int i=0; i < n; i++) {\r
 517       for (int j=0; j < m; j++) {\r
 518         if (mat[i][j] >= max) {\r
 519           max = mat[i][j];\r
 520         }\r
 521         if (mat[i][j] <= min) {\r
 522           min = mat[i][j];\r
 523         }\r
 524       }\r
 525     }\r
 526     System.out.println(max + " " + min);\r
 527     for (int i=0; i < n; i++) {\r
 528       for (int j=0; j < m; j++) {\r
 529         int x = psize*i;\r
 530         int y = psize*j;\r
 531         //      System.out.println(mat[i][j]);\r
 532         float score = (float)(mat[i][j] - min)/(float)(max-min);\r
 533         g.setColor(new Color(score,0,0));\r
 534         g.fillRect(x,y,psize,psize);\r
 535         //      System.out.println(x + " " + y + " " + score);\r
 536       }\r
 537 \r
 538     }\r
 539   }\r
 540 \r
 541 }\r