src/jalview/analysis/AlignSeq.java

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