src/jalview/ext/jmol/PDBFileWithJmol.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.ext.jmol;
  22
  23 import jalview.datamodel.AlignmentAnnotation;
  24 import jalview.datamodel.Annotation;
  25 import jalview.datamodel.PDBEntry;
  26 import jalview.datamodel.Sequence;
  27 import jalview.datamodel.SequenceI;
  28 import jalview.io.AlignFile;
  29 import jalview.io.FileParse;
  30 import jalview.schemes.ResidueProperties;
  31 import jalview.util.Comparison;
  32 import jalview.util.MessageManager;
  33
  34 import java.io.IOException;
  35 import java.util.ArrayList;
  36 import java.util.Hashtable;
  37 import java.util.List;
  38 import java.util.Map;
  39
  40 import javajs.awt.Dimension;
  41
  42 import org.jmol.api.JmolStatusListener;
  43 import org.jmol.api.JmolViewer;
  44 import org.jmol.c.CBK;
  45 import org.jmol.c.STR;
  46 import org.jmol.modelset.Group;
  47 import org.jmol.modelset.Model;
  48 import org.jmol.modelset.ModelSet;
  49 import org.jmol.modelsetbio.BioModel;
  50 import org.jmol.modelsetbio.BioPolymer;
  51 import org.jmol.modelsetbio.Monomer;
  52 import org.jmol.viewer.Viewer;
  53
  54 /**
  55  * Import and process PDB files with Jmol
  56  *
  57  * @author jprocter
  58  *
  59  */
  60 public class PDBFileWithJmol extends AlignFile implements
  61         JmolStatusListener
  62 {
  63   Viewer viewer = null;
  64
  65   public PDBFileWithJmol(String inFile, String type) throws IOException
  66   {
  67     super(inFile, type);
  68   }
  69
  70   public PDBFileWithJmol(FileParse fp) throws IOException
  71   {
  72     super(fp);
  73   }
  74
  75   public PDBFileWithJmol()
  76   {
  77   }
  78
  79   /**
  80    * create a headless jmol instance for dataprocessing
  81    *
  82    * @return
  83    */
  84   private Viewer getJmolData()
  85   {
  86     if (viewer == null)
  87     {
  88       try
  89       {
  90         viewer = (Viewer) JmolViewer.allocateViewer(null, null, null, null,
  91                 null, "-x -o -n", this);
  92         // ensure the 'new' (DSSP) not 'old' (Ramachandran) SS method is used
  93         viewer.setBooleanProperty("defaultStructureDSSP", true);
  94       } catch (ClassCastException x)
  95       {
  96         throw new Error(MessageManager.formatMessage(
  97                 "error.jmol_version_not_compatible_with_jalview_version",
  98                 new String[] { JmolViewer.getJmolVersion() }), x);
  99       }
 100     }
 101     return viewer;
 102   }
 103
 104   private void waitForScript(Viewer jmd)
 105   {
 106     while (jmd.isScriptExecuting())
 107     {
 108       try
 109       {
 110         Thread.sleep(50);
 111
 112       } catch (InterruptedException x)
 113       {
 114       }
 115     }
 116   }
 117
 118   /**
 119    * Convert Jmol's secondary structure code to Jalview's, and stored it in the
 120    * secondary structure arrays at the given sequence position
 121    *
 122    * @param proteinStructureSubType
 123    * @param pos
 124    * @param secstr
 125    * @param secstrcode
 126    */
 127   protected void setSecondaryStructure(STR proteinStructureSubType,
 128           int pos, char[] secstr, char[] secstrcode)
 129   {
 130     switch (proteinStructureSubType)
 131     {
 132     case HELIX310:
 133       if (secstr[pos] == 0)
 134       {
 135         secstr[pos] = '3';
 136       }
 137     case HELIXALPHA:
 138       if (secstr[pos] == 0)
 139       {
 140         secstr[pos] = 'H';
 141       }
 142     case HELIXPI:
 143       if (secstr[pos] == 0)
 144       {
 145         secstr[pos] = 'P';
 146       }
 147     case HELIX:
 148       if (secstr[pos] == 0)
 149       {
 150         secstr[pos] = 'H';
 151       }
 152       secstrcode[pos] = 'H';
 153       break;
 154     case SHEET:
 155       secstr[pos] = 'E';
 156       secstrcode[pos] = 'E';
 157       break;
 158     default:
 159       secstr[pos] = 0;
 160       secstrcode[pos] = 0;
 161     }
 162   }
 163
 164   /**
 165    * Convert any non-standard peptide codes to their standard code table
 166    * equivalent. (Initial version only does Selenomethionine MSE->MET.)
 167    *
 168    * @param threeLetterCode
 169    * @param seq
 170    * @param pos
 171    */
 172   protected void replaceNonCanonicalResidue(String threeLetterCode,
 173           char[] seq, int pos)
 174   {
 175     String canonical = ResidueProperties
 176             .getCanonicalAminoAcid(threeLetterCode);
 177     if (canonical != null && !canonical.equalsIgnoreCase(threeLetterCode))
 178     {
 179       seq[pos] = ResidueProperties.getSingleCharacterCode(canonical);
 180     }
 181   }
 182
 183   /**
 184    * Not implemented - returns null
 185    */
 186   @Override
 187   public String print()
 188   {
 189     return null;
 190   }
 191
 192   /**
 193    * Not implemented
 194    */
 195   @Override
 196   public void setCallbackFunction(String callbackType,
 197           String callbackFunction)
 198   {
 199   }
 200
 201   @Override
 202   public void notifyCallback(CBK cbType, Object[] data)
 203   {
 204     String strInfo = (data == null || data[1] == null ? null : data[1]
 205             .toString());
 206     switch (cbType)
 207     {
 208     case ECHO:
 209       sendConsoleEcho(strInfo);
 210       break;
 211     case SCRIPT:
 212       notifyScriptTermination((String) data[2],
 213               ((Integer) data[3]).intValue());
 214       break;
 215     case MEASURE:
 216       String mystatus = (String) data[3];
 217       if (mystatus.indexOf("Picked") >= 0
 218               || mystatus.indexOf("Sequence") >= 0)
 219       {
 220         // Picking mode
 221         sendConsoleMessage(strInfo);
 222       }
 223       else if (mystatus.indexOf("Completed") >= 0)
 224       {
 225         sendConsoleEcho(strInfo.substring(strInfo.lastIndexOf(",") + 2,
 226                 strInfo.length() - 1));
 227       }
 228       break;
 229     case MESSAGE:
 230       sendConsoleMessage(data == null ? null : strInfo);
 231       break;
 232     case PICK:
 233       sendConsoleMessage(strInfo);
 234       break;
 235     default:
 236       break;
 237     }
 238   }
 239
 240   String lastConsoleEcho = "";
 241
 242   private void sendConsoleEcho(String string)
 243   {
 244     lastConsoleEcho += string;
 245     lastConsoleEcho += "\n";
 246   }
 247
 248   String lastConsoleMessage = "";
 249
 250   private void sendConsoleMessage(String string)
 251   {
 252     lastConsoleMessage += string;
 253     lastConsoleMessage += "\n";
 254   }
 255
 256   int lastScriptTermination = -1;
 257
 258   String lastScriptMessage = "";
 259
 260   private void notifyScriptTermination(String string, int intValue)
 261   {
 262     lastScriptMessage += string;
 263     lastScriptMessage += "\n";
 264     lastScriptTermination = intValue;
 265   }
 266
 267   @Override
 268   public boolean notifyEnabled(CBK callbackPick)
 269   {
 270     switch (callbackPick)
 271     {
 272     case MESSAGE:
 273     case SCRIPT:
 274     case ECHO:
 275     case LOADSTRUCT:
 276     case ERROR:
 277       return true;
 278     default:
 279       return false;
 280     }
 281   }
 282
 283   /**
 284    * Not implemented - returns null
 285    */
 286   @Override
 287   public String eval(String strEval)
 288   {
 289     return null;
 290   }
 291
 292   /**
 293    * Not implemented - returns null
 294    */
 295   @Override
 296   public float[][] functionXY(String functionName, int x, int y)
 297   {
 298     return null;
 299   }
 300
 301   /**
 302    * Not implemented - returns null
 303    */
 304   @Override
 305   public float[][][] functionXYZ(String functionName, int nx, int ny, int nz)
 306   {
 307     return null;
 308   }
 309
 310   /**
 311    * Not implemented - returns null
 312    */
 313   @Override
 314   public String createImage(String fileName, String imageType,
 315           Object text_or_bytes, int quality)
 316   {
 317     return null;
 318   }
 319
 320   /**
 321    * Not implemented - returns null
 322    */
 323   @Override
 324   public Map<String, Object> getRegistryInfo()
 325   {
 326     return null;
 327   }
 328
 329   /**
 330    * Not implemented
 331    */
 332   @Override
 333   public void showUrl(String url)
 334   {
 335   }
 336
 337   /**
 338    * Not implemented - returns null
 339    */
 340   @Override
 341   public Dimension resizeInnerPanel(String data)
 342   {
 343     return null;
 344   }
 345
 346   @Override
 347   public Map<String, Object> getJSpecViewProperty(String arg0)
 348   {
 349     return null;
 350   }
 351
 352   /**
 353    * Calls the Jmol library to parse the PDB file, and then inspects the
 354    * resulting object model to generate Jalview-style sequences, with secondary
 355    * structure annotation added where available (i.e. where it has been computed
 356    * by Jmol using DSSP).
 357    *
 358    * @see jalview.io.AlignFile#parse()
 359    */
 360   @Override
 361   public void parse() throws IOException
 362   {
 363     Viewer jmolModel = getJmolData();
 364     jmolModel.openReader(getDataName(), getDataName(), getReader());
 365     waitForScript(jmolModel);
 366
 367     /*
 368      * Convert one or more Jmol Model objects to Jalview objects.
 369      */
 370     if (jmolModel.ms.mc > 0)
 371     {
 372       parseBiopolymers(jmolModel.ms);
 373     }
 374   }
 375
 376   /**
 377    * Process the Jmol BioPolymer array and generate a Jalview sequence for each
 378    * chain found (including any secondary structure annotation from DSSP)
 379    *
 380    * @param ms
 381    * @throws IOException
 382    */
 383   public void parseBiopolymers(ModelSet ms) throws IOException
 384   {
 385     int modelIndex = -1;
 386     for (Model model : ms.am)
 387     {
 388       modelIndex++;
 389       String modelTitle = (String) ms.getInfo(modelIndex, "title");
 390
 391       /*
 392        * as chains can span BioPolymers, we first make a flattened list,
 393        * and then work out the lengths of chains present
 394        */
 395       List<Monomer> monomers = getMonomers(ms, (BioModel) model);
 396       List<Integer> chainLengths = getChainLengths(monomers);
 397
 398       /*
 399        * now chop up the Monomer list to make Jalview Sequences
 400        */
 401       int from = 0;
 402       for (int length : chainLengths)
 403       {
 404         buildSequenceFromChain(monomers.subList(from, from + length), modelTitle);
 405         from += length;
 406       }
 407     }
 408   }
 409
 410   /**
 411    * Helper method to construct a sequence for one chain and add it to the seqs
 412    * list
 413    *
 414    * @param monomers
 415    *          a list of all monomers in the chain
 416    * @param modelTitle
 417    */
 418   protected void buildSequenceFromChain(List<Monomer> monomers, String modelTitle)
 419   {
 420     final int length = monomers.size();
 421
 422     /*
 423      * arrays to hold sequence and secondary structure
 424      */
 425     char[] seq = new char[length];
 426     char[] secstr = new char[length];
 427     char[] secstrcode = new char[length];
 428
 429     /*
 430      * populate the sequence and secondary structure arrays
 431      */
 432     extractJmolChainData(monomers, seq, secstr, secstrcode);
 433
 434     /*
 435      * grab chain code and start position from first residue;
 436      */
 437     String chainId = monomers.get(0).chain.getIDStr();
 438     int firstResNum = monomers.get(0).getResno();
 439     if (firstResNum < 1)
 440     {
 441       // Jalview doesn't like residue < 1, so force this to 1
 442       System.err.println("Converting chain " + chainId + " first RESNUM ("
 443               + firstResNum + ") to 1");
 444       firstResNum = 1;
 445     }
 446
 447     /*
 448      * convert any non-gap unknown residues to 'X'
 449      */
 450     convertNonGapCharacters(seq);
 451
 452     /*
 453      * construct and add the Jalview sequence
 454      */
 455     SequenceI sq = new Sequence("" + getDataName() + "|" + modelTitle + "|"
 456             + chainId, seq, firstResNum, firstResNum + length - 1);
 457     seqs.add(sq);
 458
 459     /*
 460      * add secondary structure predictions (if any)
 461      */
 462     addSecondaryStructureAnnotation(modelTitle, sq, secstr, secstrcode,
 463             chainId, firstResNum);
 464
 465     /*
 466      * record the PDB id for the sequence
 467      */
 468     addPdbid(sq, chainId);
 469   }
 470
 471   /**
 472    * Scans the list of (Jmol) Monomer objects, and adds the residue for each to
 473    * the sequence array, and any converted secondary structure prediction to the
 474    * secondary structure arrays
 475    *
 476    * @param monomers
 477    * @param seq
 478    * @param secstr
 479    * @param secstrcode
 480    */
 481   protected void extractJmolChainData(List<Monomer> monomers, char[] seq,
 482           char[] secstr, char[] secstrcode)
 483   {
 484     int pos = 0;
 485     for (Monomer monomer : monomers)
 486     {
 487       seq[pos] = monomer.getGroup1();
 488
 489       /*
 490        * JAL-1828 replace a modified amino acid with its standard
 491        * equivalent (e.g. MSE with MET->M) to maximise sequence matching
 492        */
 493       replaceNonCanonicalResidue(monomer.getGroup3(), seq, pos);
 494
 495       /*
 496        * if Jmol has derived a secondary structure prediction for
 497        * this position, convert it to Jalview equivalent and save it
 498        */
 499       setSecondaryStructure(monomer.getProteinStructureSubType(), pos,
 500               secstr, secstrcode);
 501       pos++;
 502     }
 503   }
 504
 505   /**
 506    * Helper method that adds an AlignmentAnnotation for secondary structure to
 507    * the sequence, provided at least one secondary structure prediction has been
 508    * made
 509    *
 510    * @param modelTitle
 511    * @param seq
 512    * @param secstr
 513    * @param secstrcode
 514    * @param chainId
 515    * @param firstResNum
 516    * @return
 517    */
 518   protected void addSecondaryStructureAnnotation(String modelTitle,
 519           SequenceI sq, char[] secstr, char[] secstrcode,
 520           String chainId, int firstResNum)
 521   {
 522     char[] seq = sq.getSequence();
 523     boolean ssFound = false;
 524     Annotation asecstr[] = new Annotation[seq.length + firstResNum - 1];
 525     for (int p = 0; p < seq.length; p++)
 526     {
 527       if (secstr[p] >= 'A' && secstr[p] <= 'z')
 528       {
 529         asecstr[p] = new Annotation(String.valueOf(secstr[p]), null,
 530                 secstrcode[p], Float.NaN);
 531         ssFound = true;
 532       }
 533     }
 534
 535     if (ssFound)
 536     {
 537       String mt = modelTitle == null ? getDataName() : modelTitle;
 538       mt += chainId;
 539       AlignmentAnnotation ann = new AlignmentAnnotation(
 540               "Secondary Structure", "Secondary Structure for " + mt,
 541               asecstr);
 542       ann.belowAlignment = true;
 543       ann.visible = true;
 544       ann.autoCalculated = false;
 545       ann.setCalcId(getClass().getName());
 546       ann.adjustForAlignment();
 547       ann.validateRangeAndDisplay();
 548       annotations.add(ann);
 549       sq.addAlignmentAnnotation(ann);
 550     }
 551   }
 552
 553   /**
 554    * Replace any non-gap miscellaneous characters with 'X'
 555    *
 556    * @param seq
 557    * @return
 558    */
 559   protected void convertNonGapCharacters(char[] seq)
 560   {
 561     boolean isNa = Comparison.areNucleotide(new char[][] { seq });
 562     int[] cinds = isNa ? ResidueProperties.nucleotideIndex
 563             : ResidueProperties.aaIndex;
 564     int nonGap = isNa ? ResidueProperties.maxNucleotideIndex
 565             : ResidueProperties.maxProteinIndex;
 566
 567     for (int p = 0; p < seq.length; p++)
 568     {
 569       if (cinds[seq[p]] == nonGap)
 570       {
 571         seq[p] = 'X';
 572       }
 573     }
 574   }
 575
 576   /**
 577    * @param sq
 578    * @param chainId
 579    */
 580   protected void addPdbid(SequenceI sq, String chainId)
 581   {
 582     PDBEntry pdbe = new PDBEntry();
 583     pdbe.setFile(getDataName());
 584     pdbe.setId(getDataName());
 585     pdbe.setProperty(new Hashtable());
 586     // pdbe.getProperty().put("CHAIN", "" + _lastChainId);
 587     pdbe.setChainCode(chainId);
 588     sq.addPDBId(pdbe);
 589   }
 590
 591   /**
 592    * Scans the list of Monomers (residue models), inspecting the chain id for
 593    * each, and returns an array whose length is the number of chains, and values
 594    * the length of each chain
 595    *
 596    * @param monomers
 597    * @return
 598    */
 599   protected List<Integer> getChainLengths(List<Monomer> monomers)
 600   {
 601     List<Integer> chainLengths = new ArrayList<Integer>();
 602     int lastChainId = -1;
 603     int length = 0;
 604
 605     for (Monomer monomer : monomers)
 606     {
 607       int chainId = monomer.chain.chainID;
 608       if (chainId != lastChainId && length > 0)
 609       {
 610         /*
 611          * change of chain - record the length of the last one
 612          */
 613         chainLengths.add(length);
 614         length = 0;
 615       }
 616       lastChainId = chainId;
 617       length++;
 618     }
 619     if (length > 0)
 620     {
 621       /*
 622        * record the length of the final chain
 623        */
 624       chainLengths.add(length);
 625     }
 626
 627     return chainLengths;
 628   }
 629
 630   /**
 631    * Returns a flattened list of Monomer (residue) in order, across all
 632    * BioPolymers in the model. This simplifies assembling chains which span
 633    * BioPolymers. The result does not include any alternate residues reported
 634    * for the same sequence position (RESNUM value).
 635    *
 636    * @param ms
 637    * @param model
 638    * @return
 639    */
 640   protected List<Monomer> getMonomers(ModelSet ms, BioModel model)
 641   {
 642     List<Monomer> result = new ArrayList<Monomer>();
 643     String lastSeqCode = "";
 644
 645     for (BioPolymer bp : model.bioPolymers) {
 646       for (int groupLeadAtoms : bp.getLeadAtomIndices())
 647       {
 648         Group group = ms.at[groupLeadAtoms].group;
 649         if (group instanceof Monomer)
 650         {
 651           /*
 652            * ignore alternate residue at same position
 653            * example: 1ejg has residues A:LEU, B:ILE, C:ILE at RESNUM=25
 654            */
 655           String seqcodeString = group.getSeqcodeString();
 656           if (!lastSeqCode.equals(seqcodeString))
 657           {
 658             result.add((Monomer) group);
 659           }
 660           else
 661           {
 662             System.out.println("skipping");
 663           }
 664           lastSeqCode = seqcodeString;
 665         }
 666       }
 667     }
 668     return result;
 669   }
 670
 671 }