JAL-629 Some refactoring to allow more flexible use of PAE files. Not working as...
[jalview.git] / 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 java.io.File;
24 import java.io.IOException;
25 import java.util.ArrayList;
26 import java.util.HashMap;
27 import java.util.List;
28 import java.util.Locale;
29 import java.util.Map;
30 import java.util.Vector;
31
32 import org.jmol.api.JmolStatusListener;
33 import org.jmol.api.JmolViewer;
34 import org.jmol.c.CBK;
35 import org.jmol.c.STR;
36 import org.jmol.modelset.ModelSet;
37 import org.jmol.viewer.Viewer;
38
39 import com.stevesoft.pat.Regex;
40
41 import jalview.bin.Console;
42 import jalview.datamodel.Alignment;
43 import jalview.datamodel.AlignmentAnnotation;
44 import jalview.datamodel.Annotation;
45 import jalview.datamodel.PDBEntry;
46 import jalview.datamodel.SequenceI;
47 import jalview.datamodel.annotations.AlphaFoldAnnotationRowBuilder;
48 import jalview.datamodel.annotations.AnnotationRowBuilder;
49 import jalview.io.DataSourceType;
50 import jalview.io.FileParse;
51 import jalview.io.StructureFile;
52 import jalview.schemes.ResidueProperties;
53 import jalview.structure.StructureImportSettings;
54 import jalview.util.Format;
55 import jalview.util.MessageManager;
56 import jalview.ws.dbsources.EBIAlfaFold;
57 import mc_view.Atom;
58 import mc_view.PDBChain;
59 import mc_view.Residue;
60
61 /**
62  * Import and process files with Jmol for file like PDB, mmCIF
63  * 
64  * @author jprocter
65  * 
66  */
67 public class JmolParser extends StructureFile implements JmolStatusListener
68 {
69   Viewer viewer = null;
70
71   public JmolParser(boolean immediate, Object inFile,
72           DataSourceType sourceType) throws IOException
73   {
74     // BH 2018 File or String for filename
75     super(immediate, inFile, sourceType);
76
77   }
78
79   public JmolParser(Object inFile, DataSourceType sourceType)
80           throws IOException
81   {
82     this(inFile, sourceType, null);
83   }
84
85   public JmolParser(Object inFile, DataSourceType sourceType,
86           StructureImportSettings.TFType tempfacType) throws IOException
87   {
88     super(inFile, sourceType, tempfacType);
89   }
90
91   public JmolParser(FileParse fp) throws IOException
92   {
93     super(fp);
94   }
95
96   public JmolParser()
97   {
98   }
99
100   /**
101    * Calls the Jmol library to parse the PDB/mmCIF file, and then inspects the
102    * resulting object model to generate Jalview-style sequences, with secondary
103    * structure annotation added where available (i.e. where it has been computed
104    * by Jmol using DSSP).
105    * 
106    * @see jalview.io.AlignFile#parse()
107    */
108   @Override
109   public void parse() throws IOException
110   {
111     setChains(new Vector<PDBChain>());
112     Viewer jmolModel = getJmolData();
113     jmolModel.openReader(getDataName(), getDataName(), getReader());
114     waitForScript(jmolModel);
115
116     /*
117      * Convert one or more Jmol Model objects to Jalview sequences
118      */
119     if (jmolModel.ms.mc > 0)
120     {
121       // ideally we do this
122       // try
123       // {
124       // setStructureFileType(jmolModel.evalString("show _fileType"));
125       // } catch (Exception q)
126       // {
127       // }
128       // ;
129       // instead, we distinguish .cif from non-.cif by filename
130       setStructureFileType(
131               getDataName().toLowerCase(Locale.ROOT).endsWith(".cif")
132                       ? PDBEntry.Type.MMCIF.toString()
133                       : "PDB");
134
135       transformJmolModelToJalview(jmolModel.ms);
136     }
137   }
138
139   /**
140    * create a headless jmol instance for dataprocessing
141    * 
142    * @return
143    */
144   private Viewer getJmolData()
145   {
146     if (viewer == null)
147     {
148       try
149       {
150         /*
151          * params -o (output to sysout) -n (nodisplay) -x (exit when finished)
152          * see http://wiki.jmol.org/index.php/Jmol_Application
153          */
154
155         viewer = JalviewJmolBinding.getJmolData(this);
156         // ensure the 'new' (DSSP) not 'old' (Ramachandran) SS method is used
157         viewer.setBooleanProperty("defaultStructureDSSP", true);
158       } catch (ClassCastException x)
159       {
160         throw new Error(MessageManager.formatMessage(
161                 "error.jmol_version_not_compatible_with_jalview_version",
162                 new String[]
163                 { JmolViewer.getJmolVersion() }), x);
164       }
165     }
166     return viewer;
167   }
168
169   public static Regex getNewAlphafoldValidator()
170   {
171     Regex validator = new Regex("(AF-[A-Z]+[0-9]+[A-Z0-9]+-F1)");
172     validator.setIgnoreCase(true);
173     return validator;
174   }
175
176   PDBEntry.Type jmolFiletype = null;
177
178   /**
179    * resolve a jmol filetype string and update the jmolFiletype field
180    * accordingly
181    * 
182    * @param jmolIdentifiedFileType
183    * @return true if filetype was identified as MMCIF, PDB
184    */
185   public boolean updateFileType(String jmolIdentifiedFileType)
186   {
187     if (jmolIdentifiedFileType == null
188             || jmolIdentifiedFileType.trim().equals(""))
189     {
190       return false;
191     }
192     if ("mmcif".equalsIgnoreCase(jmolIdentifiedFileType))
193     {
194       jmolFiletype = PDBEntry.Type.MMCIF;
195       return true;
196     }
197     if ("pdb".equalsIgnoreCase(jmolIdentifiedFileType))
198     {
199       jmolFiletype = PDBEntry.Type.PDB;
200       return true;
201     }
202     return false;
203   }
204
205   public void transformJmolModelToJalview(ModelSet ms) throws IOException
206   {
207     try
208     {
209       Regex alphaFold = getNewAlphafoldValidator();
210       String lastID = "";
211       List<SequenceI> rna = new ArrayList<SequenceI>();
212       List<SequenceI> prot = new ArrayList<SequenceI>();
213       PDBChain tmpchain;
214       String pdbId = (String) ms.getInfo(0, "title");
215       boolean isMMCIF = false;
216       String jmolFileType_String = (String) ms.getInfo(0, "fileType");
217       if (updateFileType(jmolFileType_String))
218       {
219         setStructureFileType(jmolFiletype.toString());
220       }
221
222       isMMCIF = PDBEntry.Type.MMCIF.equals(jmolFiletype);
223
224       if (pdbId == null)
225       {
226         setId(safeName(getDataName()));
227         setPDBIdAvailable(false);
228       }
229       else
230       {
231         setId(pdbId);
232         setPDBIdAvailable(true);
233         Console.debug("##### DATASOURCETYPE=" + getDataSourceType());
234         setAlphafoldModel(alphaFold.search(pdbId) && isMMCIF
235                 && getDataSourceType() == DataSourceType.URL);
236       }
237       List<Atom> significantAtoms = convertSignificantAtoms(ms);
238       for (Atom tmpatom : significantAtoms)
239       {
240         if (tmpatom.resNumIns.trim().equals(lastID))
241         {
242           // phosphorylated protein - seen both CA and P..
243           continue;
244         }
245         tmpchain = findChain(tmpatom.chain);
246         if (tmpchain != null)
247         {
248           tmpchain.atoms.addElement(tmpatom);
249         }
250         else
251         {
252           AnnotationRowBuilder builder = null;
253           if (isAlphafoldModel()
254                   || getTemperatureFactorType() == StructureImportSettings.TFType.PLDDT)
255           {
256             builder = new AlphaFoldAnnotationRowBuilder();
257           }
258
259           tmpchain = new PDBChain(getId(), tmpatom.chain, builder);
260           getChains().add(tmpchain);
261           tmpchain.atoms.addElement(tmpatom);
262         }
263         lastID = tmpatom.resNumIns.trim();
264       }
265       if (isParseImmediately())
266       {
267         // configure parsing settings from the static singleton
268         xferSettings();
269       }
270
271       makeResidueList();
272       makeCaBondList();
273
274       for (PDBChain chain : getChains())
275       {
276         SequenceI chainseq = postProcessChain(chain);
277         if (isRNA(chainseq))
278         {
279           rna.add(chainseq);
280         }
281         else
282         {
283           prot.add(chainseq);
284         }
285
286         // look at local setting for adding secondary tructure
287         if (predictSecondaryStructure)
288         {
289           createAnnotation(chainseq, chain, ms.at);
290         }
291       }
292       // if Alphafold, fetch the PAE matrix if doesn't already have one
293       if (isAlphafoldModel() && !hasPAEMatrix())
294       {
295         try
296         {
297           Console.info("retrieving pAE for " + pdbId);
298           File paeFile = EBIAlfaFold.fetchAlphaFoldPAE(pdbId, null);
299           this.setPAEMatrix(paeFile.getAbsolutePath());
300         } catch (Throwable t)
301         {
302           Console.error("Couldn't get the pAE for " + pdbId, t);
303         }
304       }
305       // add a PAEMatrix if set (either by above or otherwise)
306       Console.debug("##### hasPAEMatrix()=" + hasPAEMatrix());
307       Console.debug("##### isAlphafoldModel()=" + isAlphafoldModel());
308       Console.debug("##### getPAEMatrix()=" + getPAEMatrix());
309       if (hasPAEMatrix())
310       {
311         Alignment al = new Alignment(prot.toArray(new SequenceI[0]));
312         if (isAlphafoldModel())
313         {
314           EBIAlfaFold.addAlphaFoldPAE(al, new File(this.getPAEMatrix()), 0,
315                   null, false, false);
316         }
317         else
318         {
319           EBIAlfaFold.addPAEToStructure(null, this.getInFile(),
320                   new File(this.getPAEMatrix()));
321         }
322
323         Console.debug("##### al.getAlignmentAnnotation()="
324                 + al.getAlignmentAnnotation());
325         if (al.getAlignmentAnnotation() != null)
326         {
327           for (AlignmentAnnotation alann : al.getAlignmentAnnotation())
328           {
329             Console.debug("##### Adding to alann" + alann.annotationId + "("
330                     + alann.getCalcId() + ")");
331             annotations.add(alann);
332           }
333         }
334       }
335     } catch (OutOfMemoryError er)
336     {
337       System.out.println(
338               "OUT OF MEMORY LOADING TRANSFORMING JMOL MODEL TO JALVIEW MODEL");
339       throw new IOException(MessageManager
340               .getString("exception.outofmemory_loading_mmcif_file"));
341     }
342   }
343
344   private List<Atom> convertSignificantAtoms(ModelSet ms)
345   {
346     List<Atom> significantAtoms = new ArrayList<Atom>();
347     HashMap<String, org.jmol.modelset.Atom> chainTerMap = new HashMap<String, org.jmol.modelset.Atom>();
348     org.jmol.modelset.Atom prevAtom = null;
349     for (org.jmol.modelset.Atom atom : ms.at)
350     {
351       if (atom.getAtomName().equalsIgnoreCase("CA")
352               || atom.getAtomName().equalsIgnoreCase("P"))
353       {
354         if (!atomValidated(atom, prevAtom, chainTerMap))
355         {
356           continue;
357         }
358         Atom curAtom = new Atom(atom.x, atom.y, atom.z);
359         curAtom.atomIndex = atom.getIndex();
360         curAtom.chain = atom.getChainIDStr();
361         curAtom.insCode = atom.group.getInsertionCode() == '\000' ? ' '
362                 : atom.group.getInsertionCode();
363         curAtom.name = atom.getAtomName();
364         curAtom.number = atom.getAtomNumber();
365         curAtom.resName = atom.getGroup3(true);
366         curAtom.resNumber = atom.getResno();
367         curAtom.occupancy = ms.occupancies != null
368                 ? ms.occupancies[atom.getIndex()]
369                 : Float.valueOf(atom.getOccupancy100());
370         String fmt = new Format("%4i").form(curAtom.resNumber);
371         curAtom.resNumIns = (fmt + curAtom.insCode);
372         curAtom.tfactor = atom.getBfactor100() / 100f;
373         curAtom.type = 0;
374         // significantAtoms.add(curAtom);
375         // ignore atoms from subsequent models
376         if (!significantAtoms.contains(curAtom))
377         {
378           significantAtoms.add(curAtom);
379         }
380         prevAtom = atom;
381       }
382     }
383     return significantAtoms;
384   }
385
386   private boolean atomValidated(org.jmol.modelset.Atom curAtom,
387           org.jmol.modelset.Atom prevAtom,
388           HashMap<String, org.jmol.modelset.Atom> chainTerMap)
389   {
390     // System.out.println("Atom: " + curAtom.getAtomNumber()
391     // + " Last atom index " + curAtom.group.lastAtomIndex);
392     if (chainTerMap == null || prevAtom == null)
393     {
394       return true;
395     }
396     String curAtomChId = curAtom.getChainIDStr();
397     String prevAtomChId = prevAtom.getChainIDStr();
398     // new chain encoutered
399     if (!prevAtomChId.equals(curAtomChId))
400     {
401       // On chain switch add previous chain termination to xTerMap if not exists
402       if (!chainTerMap.containsKey(prevAtomChId))
403       {
404         chainTerMap.put(prevAtomChId, prevAtom);
405       }
406       // if current atom belongs to an already terminated chain and the resNum
407       // diff < 5 then mark as valid and update termination Atom
408       if (chainTerMap.containsKey(curAtomChId))
409       {
410         if (curAtom.getResno() < chainTerMap.get(curAtomChId).getResno())
411         {
412           return false;
413         }
414         if ((curAtom.getResno()
415                 - chainTerMap.get(curAtomChId).getResno()) < 5)
416         {
417           chainTerMap.put(curAtomChId, curAtom);
418           return true;
419         }
420         return false;
421       }
422     }
423     // atom with previously terminated chain encountered
424     else if (chainTerMap.containsKey(curAtomChId))
425     {
426       if (curAtom.getResno() < chainTerMap.get(curAtomChId).getResno())
427       {
428         return false;
429       }
430       if ((curAtom.getResno()
431               - chainTerMap.get(curAtomChId).getResno()) < 5)
432       {
433         chainTerMap.put(curAtomChId, curAtom);
434         return true;
435       }
436       return false;
437     }
438     // HETATM with resNum jump > 2
439     return !(curAtom.isHetero()
440             && ((curAtom.getResno() - prevAtom.getResno()) > 2));
441   }
442
443   private void createAnnotation(SequenceI sequence, PDBChain chain,
444           org.jmol.modelset.Atom[] jmolAtoms)
445   {
446     char[] secstr = new char[sequence.getLength()];
447     char[] secstrcode = new char[sequence.getLength()];
448
449     // Ensure Residue size equals Seq size
450     if (chain.residues.size() != sequence.getLength())
451     {
452       return;
453     }
454     int annotIndex = 0;
455     for (Residue residue : chain.residues)
456     {
457       Atom repAtom = residue.getAtoms().get(0);
458       STR proteinStructureSubType = jmolAtoms[repAtom.atomIndex].group
459               .getProteinStructureSubType();
460       setSecondaryStructure(proteinStructureSubType, annotIndex, secstr,
461               secstrcode);
462       ++annotIndex;
463     }
464     addSecondaryStructureAnnotation(chain.pdbid, sequence, secstr,
465             secstrcode, chain.id, sequence.getStart());
466   }
467
468   /**
469    * Helper method that adds an AlignmentAnnotation for secondary structure to
470    * the sequence, provided at least one secondary structure assignment has been
471    * made
472    * 
473    * @param modelTitle
474    * @param seq
475    * @param secstr
476    * @param secstrcode
477    * @param chainId
478    * @param firstResNum
479    * @return
480    */
481   protected void addSecondaryStructureAnnotation(String modelTitle,
482           SequenceI sq, char[] secstr, char[] secstrcode, String chainId,
483           int firstResNum)
484   {
485     int length = sq.getLength();
486     boolean ssFound = false;
487     Annotation asecstr[] = new Annotation[length + firstResNum - 1];
488     for (int p = 0; p < length; p++)
489     {
490       if (secstr[p] >= 'A' && secstr[p] <= 'z')
491       {
492         try
493         {
494           asecstr[p] = new Annotation(null, null, secstrcode[p], Float.NaN);
495           ssFound = true;
496         } catch (Exception e)
497         {
498           // e.printStackTrace();
499         }
500       }
501     }
502
503     if (ssFound)
504     {
505       String mt = modelTitle == null ? getDataName() : modelTitle;
506       mt += chainId;
507       AlignmentAnnotation ann = new AlignmentAnnotation(
508               "Secondary Structure", "Secondary Structure for " + mt,
509               asecstr);
510       ann.belowAlignment = true;
511       ann.visible = true;
512       ann.autoCalculated = false;
513       ann.setCalcId(getClass().getName());
514       ann.adjustForAlignment();
515       ann.validateRangeAndDisplay();
516       annotations.add(ann);
517       sq.addAlignmentAnnotation(ann);
518     }
519   }
520
521   private void waitForScript(Viewer jmd)
522   {
523     while (jmd.isScriptExecuting())
524     {
525       try
526       {
527         Thread.sleep(50);
528
529       } catch (InterruptedException x)
530       {
531       }
532     }
533   }
534
535   /**
536    * Convert Jmol's secondary structure code to Jalview's, and stored it in the
537    * secondary structure arrays at the given sequence position
538    * 
539    * @param proteinStructureSubType
540    * @param pos
541    * @param secstr
542    * @param secstrcode
543    */
544   protected void setSecondaryStructure(STR proteinStructureSubType, int pos,
545           char[] secstr, char[] secstrcode)
546   {
547     switch (proteinStructureSubType)
548     {
549     case HELIX310:
550       secstr[pos] = '3';
551       break;
552     case HELIX:
553     case HELIXALPHA:
554       secstr[pos] = 'H';
555       break;
556     case HELIXPI:
557       secstr[pos] = 'P';
558       break;
559     case SHEET:
560       secstr[pos] = 'E';
561       break;
562     default:
563       secstr[pos] = 0;
564     }
565
566     switch (proteinStructureSubType)
567     {
568     case HELIX310:
569     case HELIXALPHA:
570     case HELIXPI:
571     case HELIX:
572       secstrcode[pos] = 'H';
573       break;
574     case SHEET:
575       secstrcode[pos] = 'E';
576       break;
577     default:
578       secstrcode[pos] = 0;
579     }
580   }
581
582   /**
583    * Convert any non-standard peptide codes to their standard code table
584    * equivalent. (Initial version only does Selenomethionine MSE->MET.)
585    * 
586    * @param threeLetterCode
587    * @param seq
588    * @param pos
589    */
590   protected void replaceNonCanonicalResidue(String threeLetterCode,
591           char[] seq, int pos)
592   {
593     String canonical = ResidueProperties
594             .getCanonicalAminoAcid(threeLetterCode);
595     if (canonical != null && !canonical.equalsIgnoreCase(threeLetterCode))
596     {
597       seq[pos] = ResidueProperties.getSingleCharacterCode(canonical);
598     }
599   }
600
601   /**
602    * Not implemented - returns null
603    */
604   @Override
605   public String print(SequenceI[] seqs, boolean jvSuffix)
606   {
607     return null;
608   }
609
610   /**
611    * Not implemented
612    */
613   @Override
614   public void setCallbackFunction(String callbackType,
615           String callbackFunction)
616   {
617   }
618
619   @Override
620   public void notifyCallback(CBK cbType, Object[] data)
621   {
622     String strInfo = (data == null || data[1] == null ? null
623             : data[1].toString());
624     switch (cbType)
625     {
626     case ECHO:
627       sendConsoleEcho(strInfo);
628       break;
629     case SCRIPT:
630       notifyScriptTermination((String) data[2],
631               ((Integer) data[3]).intValue());
632       break;
633     case MEASURE:
634       String mystatus = (String) data[3];
635       if (mystatus.indexOf("Picked") >= 0
636               || mystatus.indexOf("Sequence") >= 0)
637       {
638         // Picking mode
639         sendConsoleMessage(strInfo);
640       }
641       else if (mystatus.indexOf("Completed") >= 0)
642       {
643         sendConsoleEcho(strInfo.substring(strInfo.lastIndexOf(",") + 2,
644                 strInfo.length() - 1));
645       }
646       break;
647     case MESSAGE:
648       sendConsoleMessage(data == null ? null : strInfo);
649       break;
650     case PICK:
651       sendConsoleMessage(strInfo);
652       break;
653     default:
654       break;
655     }
656   }
657
658   String lastConsoleEcho = "";
659
660   private void sendConsoleEcho(String string)
661   {
662     lastConsoleEcho += string;
663     lastConsoleEcho += "\n";
664   }
665
666   String lastConsoleMessage = "";
667
668   private void sendConsoleMessage(String string)
669   {
670     lastConsoleMessage += string;
671     lastConsoleMessage += "\n";
672   }
673
674   int lastScriptTermination = -1;
675
676   String lastScriptMessage = "";
677
678   private void notifyScriptTermination(String string, int intValue)
679   {
680     lastScriptMessage += string;
681     lastScriptMessage += "\n";
682     lastScriptTermination = intValue;
683   }
684
685   @Override
686   public boolean notifyEnabled(CBK callbackPick)
687   {
688     switch (callbackPick)
689     {
690     case MESSAGE:
691     case SCRIPT:
692     case ECHO:
693     case LOADSTRUCT:
694     case ERROR:
695       return true;
696     default:
697       return false;
698     }
699   }
700
701   /**
702    * Not implemented - returns null
703    */
704   @Override
705   public String eval(String strEval)
706   {
707     return null;
708   }
709
710   /**
711    * Not implemented - returns null
712    */
713   @Override
714   public float[][] functionXY(String functionName, int x, int y)
715   {
716     return null;
717   }
718
719   /**
720    * Not implemented - returns null
721    */
722   @Override
723   public float[][][] functionXYZ(String functionName, int nx, int ny,
724           int nz)
725   {
726     return null;
727   }
728
729   /**
730    * Not implemented - returns null
731    */
732   @Override
733   public String createImage(String fileName, String imageType,
734           Object text_or_bytes, int quality)
735   {
736     return null;
737   }
738
739   /**
740    * Not implemented - returns null
741    */
742   @Override
743   public Map<String, Object> getRegistryInfo()
744   {
745     return null;
746   }
747
748   /**
749    * Not implemented
750    */
751   @Override
752   public void showUrl(String url)
753   {
754   }
755
756   /**
757    * Not implemented - returns null
758    */
759   @Override
760   public int[] resizeInnerPanel(String data)
761   {
762     return null;
763   }
764
765   @Override
766   public Map<String, Object> getJSpecViewProperty(String arg0)
767   {
768     return null;
769   }
770
771   public boolean isPredictSecondaryStructure()
772   {
773     return predictSecondaryStructure;
774   }
775
776   public void setPredictSecondaryStructure(
777           boolean predictSecondaryStructure)
778   {
779     this.predictSecondaryStructure = predictSecondaryStructure;
780   }
781
782   public boolean isVisibleChainAnnotation()
783   {
784     return visibleChainAnnotation;
785   }
786
787   public void setVisibleChainAnnotation(boolean visibleChainAnnotation)
788   {
789     this.visibleChainAnnotation = visibleChainAnnotation;
790   }
791
792 }