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       secstr[pos] = '3';
 134       break;
 135     case HELIX:
 136     case HELIXALPHA:
 137       secstr[pos] = 'H';
 138       break;
 139     case HELIXPI:
 140       secstr[pos] = 'P';
 141       break;
 142     case SHEET:
 143       secstr[pos] = 'E';
 144       break;
 145     default:
 146       secstr[pos] = 0;
 147     }
 148
 149     switch (proteinStructureSubType)
 150     {
 151     case HELIX310:
 152     case HELIXALPHA:
 153     case HELIXPI:
 154     case HELIX:
 155       secstrcode[pos] = 'H';
 156       break;
 157     case SHEET:
 158       secstrcode[pos] = 'E';
 159       break;
 160     default:
 161       secstrcode[pos] = 0;
 162     }
 163   }
 164
 165   /**
 166    * Convert any non-standard peptide codes to their standard code table
 167    * equivalent. (Initial version only does Selenomethionine MSE->MET.)
 168    *
 169    * @param threeLetterCode
 170    * @param seq
 171    * @param pos
 172    */
 173   protected void replaceNonCanonicalResidue(String threeLetterCode,
 174           char[] seq, int pos)
 175   {
 176     String canonical = ResidueProperties
 177             .getCanonicalAminoAcid(threeLetterCode);
 178     if (canonical != null && !canonical.equalsIgnoreCase(threeLetterCode))
 179     {
 180       seq[pos] = ResidueProperties.getSingleCharacterCode(canonical);
 181     }
 182   }
 183
 184   /**
 185    * Not implemented - returns null
 186    */
 187   @Override
 188   public String print()
 189   {
 190     return null;
 191   }
 192
 193   /**
 194    * Not implemented
 195    */
 196   @Override
 197   public void setCallbackFunction(String callbackType,
 198           String callbackFunction)
 199   {
 200   }
 201
 202   @Override
 203   public void notifyCallback(CBK cbType, Object[] data)
 204   {
 205     String strInfo = (data == null || data[1] == null ? null : data[1]
 206             .toString());
 207     switch (cbType)
 208     {
 209     case ECHO:
 210       sendConsoleEcho(strInfo);
 211       break;
 212     case SCRIPT:
 213       notifyScriptTermination((String) data[2],
 214               ((Integer) data[3]).intValue());
 215       break;
 216     case MEASURE:
 217       String mystatus = (String) data[3];
 218       if (mystatus.indexOf("Picked") >= 0
 219               || mystatus.indexOf("Sequence") >= 0)
 220       {
 221         // Picking mode
 222         sendConsoleMessage(strInfo);
 223       }
 224       else if (mystatus.indexOf("Completed") >= 0)
 225       {
 226         sendConsoleEcho(strInfo.substring(strInfo.lastIndexOf(",") + 2,
 227                 strInfo.length() - 1));
 228       }
 229       break;
 230     case MESSAGE:
 231       sendConsoleMessage(data == null ? null : strInfo);
 232       break;
 233     case PICK:
 234       sendConsoleMessage(strInfo);
 235       break;
 236     default:
 237       break;
 238     }
 239   }
 240
 241   String lastConsoleEcho = "";
 242
 243   private void sendConsoleEcho(String string)
 244   {
 245     lastConsoleEcho += string;
 246     lastConsoleEcho += "\n";
 247   }
 248
 249   String lastConsoleMessage = "";
 250
 251   private void sendConsoleMessage(String string)
 252   {
 253     lastConsoleMessage += string;
 254     lastConsoleMessage += "\n";
 255   }
 256
 257   int lastScriptTermination = -1;
 258
 259   String lastScriptMessage = "";
 260
 261   private void notifyScriptTermination(String string, int intValue)
 262   {
 263     lastScriptMessage += string;
 264     lastScriptMessage += "\n";
 265     lastScriptTermination = intValue;
 266   }
 267
 268   @Override
 269   public boolean notifyEnabled(CBK callbackPick)
 270   {
 271     switch (callbackPick)
 272     {
 273     case MESSAGE:
 274     case SCRIPT:
 275     case ECHO:
 276     case LOADSTRUCT:
 277     case ERROR:
 278       return true;
 279     default:
 280       return false;
 281     }
 282   }
 283
 284   /**
 285    * Not implemented - returns null
 286    */
 287   @Override
 288   public String eval(String strEval)
 289   {
 290     return null;
 291   }
 292
 293   /**
 294    * Not implemented - returns null
 295    */
 296   @Override
 297   public float[][] functionXY(String functionName, int x, int y)
 298   {
 299     return null;
 300   }
 301
 302   /**
 303    * Not implemented - returns null
 304    */
 305   @Override
 306   public float[][][] functionXYZ(String functionName, int nx, int ny, int nz)
 307   {
 308     return null;
 309   }
 310
 311   /**
 312    * Not implemented - returns null
 313    */
 314   @Override
 315   public String createImage(String fileName, String imageType,
 316           Object text_or_bytes, int quality)
 317   {
 318     return null;
 319   }
 320
 321   /**
 322    * Not implemented - returns null
 323    */
 324   @Override
 325   public Map<String, Object> getRegistryInfo()
 326   {
 327     return null;
 328   }
 329
 330   /**
 331    * Not implemented
 332    */
 333   @Override
 334   public void showUrl(String url)
 335   {
 336   }
 337
 338   /**
 339    * Not implemented - returns null
 340    */
 341   @Override
 342   public Dimension resizeInnerPanel(String data)
 343   {
 344     return null;
 345   }
 346
 347   @Override
 348   public Map<String, Object> getJSpecViewProperty(String arg0)
 349   {
 350     return null;
 351   }
 352
 353   /**
 354    * Calls the Jmol library to parse the PDB file, and then inspects the
 355    * resulting object model to generate Jalview-style sequences, with secondary
 356    * structure annotation added where available (i.e. where it has been computed
 357    * by Jmol using DSSP).
 358    *
 359    * @see jalview.io.AlignFile#parse()
 360    */
 361   @Override
 362   public void parse() throws IOException
 363   {
 364     Viewer jmolModel = getJmolData();
 365     jmolModel.openReader(getDataName(), getDataName(), getReader());
 366     waitForScript(jmolModel);
 367
 368     /*
 369      * Convert one or more Jmol Model objects to Jalview sequences
 370      */
 371     if (jmolModel.ms.mc > 0)
 372     {
 373       parseBiopolymers(jmolModel.ms);
 374     }
 375   }
 376
 377   /**
 378    * Process the Jmol BioPolymer array and generate a Jalview sequence for each
 379    * chain found (including any secondary structure annotation from DSSP)
 380    *
 381    * @param ms
 382    * @throws IOException
 383    */
 384   public void parseBiopolymers(ModelSet ms) throws IOException
 385   {
 386     int modelIndex = -1;
 387     for (Model model : ms.am)
 388     {
 389       modelIndex++;
 390       String modelTitle = (String) ms.getInfo(modelIndex, "title");
 391
 392       /*
 393        * Chains can span BioPolymers, so first make a flattened list,
 394        * and then work out the lengths of chains present
 395        */
 396       List<Monomer> monomers = getMonomers(ms, (BioModel) model);
 397       List<Integer> chainLengths = getChainLengths(monomers);
 398
 399       /*
 400        * now chop up the Monomer list to make Jalview Sequences
 401        */
 402       int from = 0;
 403       for (int length : chainLengths)
 404       {
 405         buildSequenceFromChain(monomers.subList(from, from + length), modelTitle);
 406         from += length;
 407       }
 408     }
 409   }
 410
 411   /**
 412    * Helper method to construct a sequence for one chain and add it to the seqs
 413    * list
 414    *
 415    * @param monomers
 416    *          a list of all monomers in the chain
 417    * @param modelTitle
 418    */
 419   protected void buildSequenceFromChain(List<Monomer> monomers, String modelTitle)
 420   {
 421     final int length = monomers.size();
 422
 423     /*
 424      * arrays to hold sequence and secondary structure
 425      */
 426     char[] seq = new char[length];
 427     char[] secstr = new char[length];
 428     char[] secstrcode = new char[length];
 429
 430     /*
 431      * populate the sequence and secondary structure arrays
 432      */
 433     extractJmolChainData(monomers, seq, secstr, secstrcode);
 434
 435     /*
 436      * grab chain code and start position from first residue;
 437      */
 438     String chainId = monomers.get(0).chain.getIDStr();
 439     int firstResNum = monomers.get(0).getResno();
 440     if (firstResNum < 1)
 441     {
 442       // Jalview doesn't like residue < 1, so force this to 1
 443       System.err.println("Converting chain " + chainId + " first RESNUM ("
 444               + firstResNum + ") to 1");
 445       firstResNum = 1;
 446     }
 447
 448     /*
 449      * convert any non-gap unknown residues to 'X'
 450      */
 451     convertNonGapCharacters(seq);
 452
 453     /*
 454      * construct and add the Jalview sequence
 455      */
 456     String seqName = "" + getDataName() + "|" + modelTitle + "|"
 457             + chainId;
 458     SequenceI sq = new Sequence(seqName, seq, firstResNum, firstResNum + length - 1);
 459     seqs.add(sq);
 460
 461     /*
 462      * add secondary structure predictions (if any)
 463      */
 464     addSecondaryStructureAnnotation(modelTitle, sq, secstr, secstrcode,
 465             chainId, firstResNum);
 466
 467     /*
 468      * record the PDB id for the sequence
 469      */
 470     addPdbid(sq, chainId);
 471   }
 472
 473   /**
 474    * Scans the list of (Jmol) Monomer objects, and adds the residue for each to
 475    * the sequence array, and any converted secondary structure prediction to the
 476    * secondary structure arrays
 477    *
 478    * @param monomers
 479    * @param seq
 480    * @param secstr
 481    * @param secstrcode
 482    */
 483   protected void extractJmolChainData(List<Monomer> monomers, char[] seq,
 484           char[] secstr, char[] secstrcode)
 485   {
 486     int pos = 0;
 487     for (Monomer monomer : monomers)
 488     {
 489       seq[pos] = monomer.getGroup1();
 490
 491       /*
 492        * JAL-1828 replace a modified amino acid with its standard
 493        * equivalent (e.g. MSE with MET->M) to maximise sequence matching
 494        */
 495       replaceNonCanonicalResidue(monomer.getGroup3(), seq, pos);
 496
 497       /*
 498        * if Jmol has derived a secondary structure prediction for
 499        * this position, convert it to Jalview equivalent and save it
 500        */
 501       setSecondaryStructure(monomer.getProteinStructureSubType(), pos,
 502               secstr, secstrcode);
 503       pos++;
 504     }
 505   }
 506
 507   /**
 508    * Helper method that adds an AlignmentAnnotation for secondary structure to
 509    * the sequence, provided at least one secondary structure prediction has been
 510    * made
 511    *
 512    * @param modelTitle
 513    * @param seq
 514    * @param secstr
 515    * @param secstrcode
 516    * @param chainId
 517    * @param firstResNum
 518    * @return
 519    */
 520   protected void addSecondaryStructureAnnotation(String modelTitle,
 521           SequenceI sq, char[] secstr, char[] secstrcode,
 522           String chainId, int firstResNum)
 523   {
 524     char[] seq = sq.getSequence();
 525     boolean ssFound = false;
 526     Annotation asecstr[] = new Annotation[seq.length + firstResNum - 1];
 527     for (int p = 0; p < seq.length; p++)
 528     {
 529       if (secstr[p] >= 'A' && secstr[p] <= 'z')
 530       {
 531         asecstr[p] = new Annotation(String.valueOf(secstr[p]), null,
 532                 secstrcode[p], Float.NaN);
 533         ssFound = true;
 534       }
 535     }
 536
 537     if (ssFound)
 538     {
 539       String mt = modelTitle == null ? getDataName() : modelTitle;
 540       mt += chainId;
 541       AlignmentAnnotation ann = new AlignmentAnnotation(
 542               "Secondary Structure", "Secondary Structure for " + mt,
 543               asecstr);
 544       ann.belowAlignment = true;
 545       ann.visible = true;
 546       ann.autoCalculated = false;
 547       ann.setCalcId(getClass().getName());
 548       ann.adjustForAlignment();
 549       ann.validateRangeAndDisplay();
 550       annotations.add(ann);
 551       sq.addAlignmentAnnotation(ann);
 552     }
 553   }
 554
 555   /**
 556    * Replace any non-gap miscellaneous characters with 'X'
 557    *
 558    * @param seq
 559    * @return
 560    */
 561   protected void convertNonGapCharacters(char[] seq)
 562   {
 563     boolean isNa = Comparison.areNucleotide(new char[][] { seq });
 564     int[] cinds = isNa ? ResidueProperties.nucleotideIndex
 565             : ResidueProperties.aaIndex;
 566     int nonGap = isNa ? ResidueProperties.maxNucleotideIndex
 567             : ResidueProperties.maxProteinIndex;
 568
 569     for (int p = 0; p < seq.length; p++)
 570     {
 571       if (cinds[seq[p]] == nonGap)
 572       {
 573         seq[p] = 'X';
 574       }
 575     }
 576   }
 577
 578   /**
 579    * Add a PDBEntry giving the source of PDB data to the sequence
 580    *
 581    * @param sq
 582    * @param chainId
 583    */
 584   protected void addPdbid(SequenceI sq, String chainId)
 585   {
 586     PDBEntry pdbe = new PDBEntry();
 587     pdbe.setFile(getDataName());
 588     pdbe.setId(getDataName());
 589     pdbe.setProperty(new Hashtable());
 590     pdbe.setChainCode(chainId);
 591     sq.addPDBId(pdbe);
 592   }
 593
 594   /**
 595    * Scans the list of Monomers (residue models), inspecting the chain id for
 596    * each, and returns an array whose length is the number of chains, and values
 597    * the length of each chain
 598    *
 599    * @param monomers
 600    * @return
 601    */
 602   protected List<Integer> getChainLengths(List<Monomer> monomers)
 603   {
 604     List<Integer> chainLengths = new ArrayList<Integer>();
 605     int lastChainId = -1;
 606     int length = 0;
 607
 608     for (Monomer monomer : monomers)
 609     {
 610       int chainId = monomer.chain.chainID;
 611       if (chainId != lastChainId && length > 0)
 612       {
 613         /*
 614          * change of chain - record the length of the last one
 615          */
 616         chainLengths.add(length);
 617         length = 0;
 618       }
 619       lastChainId = chainId;
 620       length++;
 621     }
 622     if (length > 0)
 623     {
 624       /*
 625        * record the length of the final chain
 626        */
 627       chainLengths.add(length);
 628     }
 629
 630     return chainLengths;
 631   }
 632
 633   /**
 634    * Returns a flattened list of Monomer (residues) in order, across all
 635    * BioPolymers in the model. This simplifies assembling chains which span
 636    * BioPolymers. The result omits any alternate residues reported for the same
 637    * sequence position (RESNUM value).
 638    *
 639    * @param ms
 640    * @param model
 641    * @return
 642    */
 643   protected List<Monomer> getMonomers(ModelSet ms, BioModel model)
 644   {
 645     List<Monomer> result = new ArrayList<Monomer>();
 646     int lastResNo = Integer.MIN_VALUE;
 647
 648     for (BioPolymer bp : model.bioPolymers)
 649     {
 650       for (int groupLeadAtoms : bp.getLeadAtomIndices())
 651       {
 652         Group group = ms.at[groupLeadAtoms].group;
 653         if (group instanceof Monomer)
 654         {
 655           /*
 656            * ignore alternate residue at same position
 657            * example: 1ejg has residues A:LEU, B:ILE at RESNUM=25
 658            */
 659           int resNo = group.getResno();
 660           if (lastResNo != resNo)
 661           {
 662             result.add((Monomer) group);
 663           }
 664           lastResNo = resNo;
 665         }
 666       }
 667     }
 668     return result;
 669   }
 670
 671 }