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