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 \r
  20 package jalview.analysis;\r
  21 \r
  22 import jalview.schemes.*;\r
  23 import jalview.datamodel.SequenceI;\r
  24 import jalview.util.*;\r
  25 import jalview.io.*;\r
  26 \r
  27 import java.util.*;\r
  28 import java.io.*;\r
  29 import java.awt.*;\r
  30 \r
  31 public class AlignSeq {\r
  32 \r
  33   int[][] score;\r
  34   int[][] E;\r
  35   int[][] F;\r
  36   int[][] traceback;\r
  37 \r
  38   int[] seq1;\r
  39   int[] seq2;\r
  40 \r
  41   SequenceI s1;\r
  42   SequenceI s2;\r
  43 \r
  44   String s1str;\r
  45   String s2str;\r
  46 \r
  47   int maxi;\r
  48   int maxj;\r
  49 \r
  50   int[] aseq1;\r
  51   int[] aseq2;\r
  52 \r
  53   String astr1 = "";\r
  54   String astr2 = "";\r
  55 \r
  56   public int seq1start;\r
  57   public int seq1end;\r
  58   public int seq2start;\r
  59   public int seq2end;\r
  60 \r
  61   int count;\r
  62 \r
  63   public int maxscore;\r
  64   float pid;\r
  65   int prev = 0;\r
  66 \r
  67   public static java.util.Hashtable dnaHash = new java.util.Hashtable();\r
  68 \r
  69   static  {\r
  70     dnaHash.put("C", new Integer(0));\r
  71     dnaHash.put("T", new Integer(1));\r
  72     dnaHash.put("A", new Integer(2));\r
  73     dnaHash.put("G", new Integer(3));\r
  74     dnaHash.put("-", new Integer(4));\r
  75   }\r
  76 \r
  77   static String dna[] = {"C","T","A","G","-"};\r
  78   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
  79 \r
  80   int gapOpen = 120;\r
  81   int gapExtend = 20;\r
  82 \r
  83   int lookup[][] = ResidueProperties.getBLOSUM62();\r
  84   String intToStr[] = pep;\r
  85   int defInt = 23;\r
  86 \r
  87   String output = "";\r
  88 \r
  89   String type;\r
  90   Runtime rt;\r
  91   public AlignSeq() {}\r
  92 \r
  93   public AlignSeq(SequenceI s1, SequenceI s2,String type) {\r
  94     rt  = Runtime.getRuntime();\r
  95     SeqInit(s1,s2,type);\r
  96   }\r
  97 \r
  98   public int getMaxScore() {\r
  99     return maxscore;\r
 100   }\r
 101 \r
 102   public int getSeq2Start() {\r
 103     return seq2start;\r
 104   }\r
 105 \r
 106   public int getSeq2End() {\r
 107     return seq2end;\r
 108   }\r
 109 \r
 110   public int getSeq1Start() {\r
 111     return seq1start;\r
 112   }\r
 113 \r
 114   public int getSeq1End() {\r
 115     return seq1end;\r
 116   }\r
 117 \r
 118   public String getOutput() {\r
 119     return output;\r
 120   }\r
 121 \r
 122   public String getAStr1() {\r
 123     return astr1;\r
 124   }\r
 125   public String getAStr2() {\r
 126     return astr2;\r
 127   }\r
 128   public int [] getASeq1() {\r
 129     return aseq1;\r
 130   }\r
 131   public int [] getASeq2() {\r
 132     return aseq2;\r
 133   }\r
 134   public SequenceI getS1() {\r
 135     return s1;\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 \r
 168     //System.out.println("lookuip " + rt.freeMemory() + " "+  rt.totalMemory());\r
 169     seq1 = new int[s1str.length()];\r
 170     //System.out.println("seq1 " + rt.freeMemory() +" "  + rt.totalMemory());\r
 171     seq2 = new int[s2str.length()];\r
 172     //System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());\r
 173     score = new int[s1str.length()][s2str.length()];\r
 174     //System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());\r
 175     E = new int[s1str.length()][s2str.length()];\r
 176     //System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());\r
 177     F = new int[s1str.length()][s2str.length()];\r
 178     traceback = new int[s1str.length()][s2str.length()];\r
 179     //System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());\r
 180     seq1 = stringToInt(s1str,type);\r
 181     //System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());\r
 182     seq2 = stringToInt(s2str,type);\r
 183     //System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());\r
 184 \r
 185     //   long tstart = System.currentTimeMillis();\r
 186     //    calcScoreMatrix();\r
 187     //long tend = System.currentTimeMillis();\r
 188 \r
 189     //System.out.println("Time take to calculate score matrix = " + (tend-tstart) + " ms");\r
 190 \r
 191 \r
 192     //   printScoreMatrix(score);\r
 193     //System.out.println();\r
 194 \r
 195     //printScoreMatrix(traceback);\r
 196     //System.out.println();\r
 197 \r
 198     //  printScoreMatrix(E);\r
 199     //System.out.println();\r
 200 \r
 201     ///printScoreMatrix(F);\r
 202     //System.out.println();\r
 203     // tstart = System.currentTimeMillis();\r
 204     //traceAlignment();\r
 205     //tend = System.currentTimeMillis();\r
 206     //System.out.println("Time take to traceback alignment = " + (tend-tstart) + " ms");\r
 207   }\r
 208 \r
 209   public void traceAlignment() {\r
 210 \r
 211     // Find the maximum score along the rhs or bottom row\r
 212     int max = -9999;\r
 213     for (int i = 0; i < seq1.length; i++) {\r
 214       if (score[i][seq2.length - 1] > max ) {\r
 215         max = score[i][seq2.length - 1];\r
 216         maxi = i;\r
 217         maxj = seq2.length-1;\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 \r
 230     int i = maxi;\r
 231     int j = maxj;\r
 232     int trace;\r
 233     maxscore = score[i][j] / 10;\r
 234 \r
 235     seq1end = maxi+1;\r
 236     seq2end = maxj+1;\r
 237 \r
 238     aseq1 = new int[seq1.length + seq2.length];\r
 239     aseq2 = new int[seq1.length + seq2.length];\r
 240 \r
 241     count = seq1.length + seq2.length - 1;\r
 242 \r
 243     while (i>0 && j >0) {\r
 244 \r
 245       if (aseq1[count] != defInt && i >=0) {\r
 246         aseq1[count] = seq1[i];\r
 247         astr1 = intToStr[seq1[i]] + astr1;\r
 248       }\r
 249 \r
 250       if (aseq2[count] != defInt && j > 0) {\r
 251         aseq2[count] = seq2[j];\r
 252         astr2 = intToStr[seq2[j]] + astr2;\r
 253       }\r
 254       trace = findTrace(i,j);\r
 255       if (trace == 0) {\r
 256         i--;\r
 257         j--;\r
 258 \r
 259       } else  if (trace == 1) {\r
 260         j--;\r
 261         aseq1[count] = defInt;\r
 262         astr1 = "-" + astr1.substring(1);\r
 263       } else  if (trace == -1) {\r
 264         i--;\r
 265         aseq2[count] = defInt;\r
 266         astr2 = "-" + astr2.substring(1);\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     int overlap = 0;\r
 297 \r
 298     output = output + ("Score = " + score[maxi][maxj] + "\n");\r
 299     output = output + ("Length of alignment = " + (aseq1.length-count) + "\n");\r
 300     output = output + ("Sequence ");\r
 301     output = output + (new Format("%" + maxid + "s").form(s1.getName()));\r
 302     output = output + (" :  " + seq1start + " - " + seq1end + " (Sequence length = " + s1str.length() + ")\n");\r
 303     output = output + ("Sequence ");\r
 304     output = output + (new Format("%" + maxid + "s").form(s2.getName()));\r
 305     output = output + (" :  " + seq2start + " - " + seq2end + " (Sequence length = " + s2str.length() + ")\n\n");\r
 306 \r
 307     for (int j = 0; j < nochunks; j++) {\r
 308       // Print the first aligned sequence\r
 309       output = output + (new Format("%" + (maxid) + "s").form(s1.getName()) + " ");\r
 310       for (int i = 0; i < len ; i++) {\r
 311 \r
 312         if ((count + i + j*len) < aseq1.length) {\r
 313           output = output + (new Format("%s").form(intToStr[aseq1[count + i + j*len]]));\r
 314         }\r
 315       }\r
 316 \r
 317       output = output + ("\n");\r
 318       output = output + (new Format("%" + (maxid) + "s").form(" ") + " ");\r
 319       // Print out the matching chars\r
 320       for (int i = 0; i < len ; i++) {\r
 321 \r
 322         if ((count + i + j*len) < aseq1.length) {\r
 323           if ( intToStr[aseq1[count+i+j*len]].equals(intToStr[aseq2[count+i+j*len]]) && !intToStr[aseq1[count+i+j*len]].equals("-")) {\r
 324             pid++;\r
 325             output = output + ("|");\r
 326           } else if (type.equals("pep")) {\r
 327             if (ResidueProperties.getPAM250(intToStr[aseq1[count+i+j*len]],intToStr[aseq2[count+i+j*len]]) > 0) {\r
 328               output  = output + (".");\r
 329             } else {\r
 330               output = output + (" ");\r
 331             }\r
 332           } else {\r
 333             output = output + (" ");\r
 334           }\r
 335 \r
 336         }\r
 337       }\r
 338       // Now print the second aligned sequence\r
 339       output = output + ("\n");\r
 340       output = output + (new Format("%" + (maxid) + "s").form(s2.getName()) + " " );\r
 341       for (int i = 0; i < len ; i++) {\r
 342         if ((count + i + j*len) < aseq1.length) {\r
 343           output = output + (new Format("%s").form(intToStr[aseq2[count + i + j*len]]));\r
 344         }\r
 345       }\r
 346       output = output + ("\n\n");\r
 347     }\r
 348     pid = pid/(float)(aseq1.length-count)*100;\r
 349     output = output + (new Format("Percentage ID = %2.2f\n\n").form(pid));\r
 350 \r
 351   }\r
 352 \r
 353   public void printScoreMatrix(int[][] mat) {\r
 354     int n = seq1.length;\r
 355     int m = seq2.length;\r
 356 \r
 357     for (int i = 0; i < n;i++) {\r
 358       // Print the top sequence\r
 359       if (i == 0) {\r
 360         Format.print(System.out,"%8s",s2str.substring(0,1));\r
 361         for (int jj = 1;jj < m; jj++) {\r
 362           Format.print(System.out,"%5s",s2str.substring(jj,jj+1));\r
 363         }\r
 364         System.out.println();\r
 365       }\r
 366 \r
 367       for (int j = 0;j < m; j++) {\r
 368         if (j == 0) {\r
 369           Format.print(System.out,"%3s",s1str.substring(i,i+1));\r
 370         }\r
 371         Format.print(System.out,"%3d ",mat[i][j]/10);\r
 372       }\r
 373       System.out.println();\r
 374     }\r
 375   }\r
 376 \r
 377   public int findTrace(int i,int j) {\r
 378     int t = 0;\r
 379     int max = score[i-1][j-1] + lookup[seq1[i]][seq2[j]] * 10;\r
 380 \r
 381     if (F[i][j] > max) {\r
 382       max = F[i][j];\r
 383       t = -1;\r
 384     } else if (F[i][j] == max) {\r
 385       if (prev == -1) {\r
 386         max = F[i][j];\r
 387         t = -1;\r
 388       }\r
 389     }\r
 390     if (E[i][j] >= max) {\r
 391       max = E[i][j];\r
 392       t = 1;\r
 393     } else if (E[i][j] == max) {\r
 394       if (prev == 1) {\r
 395         max = E[i][j];\r
 396         t = 1;\r
 397       }\r
 398     }\r
 399     prev = t;\r
 400     return t;\r
 401   }\r
 402 \r
 403   public void calcScoreMatrix() {\r
 404 \r
 405 \r
 406     int n = seq1.length;\r
 407     int m = seq2.length;\r
 408 \r
 409 \r
 410     // top left hand element\r
 411     score[0][0] =  lookup[seq1[0]][seq2[0]] * 10;\r
 412     E[0][0] = -gapExtend;\r
 413     F[0][0] = 0;\r
 414 \r
 415     // Calculate the top row first\r
 416     for (int j=1; j < m; j++) {\r
 417       // What should these values be? 0 maybe\r
 418       E[0][j] = max(score[0][j-1] - gapOpen,E[0][j-1] - gapExtend);\r
 419       F[0][j] = -gapExtend;\r
 420 \r
 421       score[0][j] = max( lookup[seq1[0]][seq2[j]] * 10 ,-gapOpen,-gapExtend);\r
 422 \r
 423       traceback[0][j] = 1;\r
 424     }\r
 425 \r
 426     // Now do the left hand column\r
 427     for (int i=1; i < n; i++) {\r
 428       E[i][0] =  -gapOpen;\r
 429       F[i][0] =  max(score[i-1][0]-gapOpen,F[i-1][0]-gapExtend);\r
 430 \r
 431       score[i][0] = max( lookup[seq1[i]][seq2[0]] * 10 ,E[i][0],F[i][0]);\r
 432       traceback[i][0] = -1;\r
 433     }\r
 434 \r
 435     // Now do all the other rows\r
 436     for (int i = 1; i < n; i++) {\r
 437       for (int j = 1; j < m; j++) {\r
 438 \r
 439         E[i][j] = max(score[i][j-1] - gapOpen,  E[i][j-1] - gapExtend);\r
 440         F[i][j] = max(score[i-1][j] - gapOpen,  F[i-1][j] - gapExtend);\r
 441 \r
 442         score[i][j] = max(score[i-1][j-1] + lookup[seq1[i]][seq2[j]]*10,\r
 443                           E[i][j],\r
 444                           F[i][j]);\r
 445         traceback[i][j] = findTrace(i,j);\r
 446       }\r
 447     }\r
 448 \r
 449   }\r
 450   public static String extractChars(String chars, String seq) {\r
 451     String out = seq;\r
 452     for (int i=0; i < chars.length(); i++) {\r
 453       String gap = chars.substring(i,i+1);\r
 454       out = extractGaps(gap,out);\r
 455     }\r
 456     return out;\r
 457   }\r
 458   public static String extractGaps(String gapChar, String seq) {\r
 459     StringTokenizer str = new StringTokenizer(seq,gapChar);\r
 460     String newString = "";\r
 461 \r
 462     while (str.hasMoreTokens()) {\r
 463       newString  = newString + str.nextToken();\r
 464     }\r
 465     return newString;\r
 466   }\r
 467 \r
 468 \r
 469   public int max(int i1, int i2, int i3) {\r
 470     int max = i1;\r
 471     if (i2 > i1)  {\r
 472       max = i2;\r
 473     }\r
 474     if (i3 > max) {\r
 475       max = i3;\r
 476     }\r
 477     return max;\r
 478   }\r
 479 \r
 480   public int max(int i1, int i2) {\r
 481     int max = i1;\r
 482     if (i2 > i1)  {\r
 483       max = i2;\r
 484     }\r
 485     return max;\r
 486   }\r
 487 \r
 488   public int[] stringToInt(String s,String type) {\r
 489     int[] seq1 = new int[s.length()];\r
 490 \r
 491     for (int i = 0;i < s.length(); i++) {\r
 492       String ss = s.substring(i,i+1).toUpperCase();\r
 493       try {\r
 494         if (type.equals("pep")) {\r
 495           seq1[i] = ((Integer)ResidueProperties.aaHash.get(ss)).intValue();\r
 496         } else if (type.equals("dna")) {\r
 497           seq1[i] = ((Integer)dnaHash.get(ss)).intValue();\r
 498         }\r
 499         if (seq1[i] > 23) {\r
 500           seq1[i] = 23;\r
 501         }\r
 502       } catch (Exception e) {\r
 503         if (type.equals("dna")) {\r
 504           seq1[i] = 4;\r
 505         } else {\r
 506           seq1[i] = 23;\r
 507         }\r
 508       }\r
 509     }\r
 510     return seq1;\r
 511   }\r
 512 \r
 513   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,int psize) {\r
 514 \r
 515     int max = -1000;\r
 516     int min = 1000;\r
 517 \r
 518     for (int i=0; i < n; i++) {\r
 519       for (int j=0; j < m; j++) {\r
 520         if (mat[i][j] >= max) {\r
 521           max = mat[i][j];\r
 522         }\r
 523         if (mat[i][j] <= min) {\r
 524           min = mat[i][j];\r
 525         }\r
 526       }\r
 527     }\r
 528     System.out.println(max + " " + min);\r
 529     for (int i=0; i < n; i++) {\r
 530       for (int j=0; j < m; j++) {\r
 531         int x = psize*i;\r
 532         int y = psize*j;\r
 533         //      System.out.println(mat[i][j]);\r
 534         float score = (float)(mat[i][j] - min)/(float)(max-min);\r
 535         g.setColor(new Color(score,0,0));\r
 536         g.fillRect(x,y,psize,psize);\r
 537         //      System.out.println(x + " " + y + " " + score);\r
 538       }\r
 539 \r
 540     }\r
 541   }\r
 542 \r
 543 }\r