JAL-2629 hmmer searches now read domain rather than full scores
[jalview.git] / src / jalview / hmmer / HMMSearch.java
index caf0b37..f05823e 100644 (file)
@@ -1,5 +1,6 @@
 package jalview.hmmer;
 
+import jalview.bin.Cache;
 import jalview.datamodel.Alignment;
 import jalview.datamodel.AlignmentAnnotation;
 import jalview.datamodel.AlignmentI;
@@ -7,36 +8,38 @@ import jalview.datamodel.Annotation;
 import jalview.datamodel.HiddenMarkovModel;
 import jalview.datamodel.SequenceI;
 import jalview.gui.AlignFrame;
+import jalview.gui.Desktop;
 import jalview.gui.JvOptionPane;
 import jalview.io.DataSourceType;
 import jalview.io.FileParse;
 import jalview.io.StockholmFile;
+import jalview.util.FileUtils;
 import jalview.util.MessageManager;
 import jalview.ws.params.ArgumentI;
 import jalview.ws.params.simple.BooleanOption;
+import jalview.ws.params.simple.Option;
 
 import java.io.BufferedReader;
 import java.io.File;
 import java.io.FileReader;
 import java.io.IOException;
-import java.io.PrintWriter;
 import java.util.ArrayList;
+import java.util.Collections;
 import java.util.List;
-import java.util.Scanner;
 
 import javax.swing.JOptionPane;
 
-public class HMMSearch extends HmmerCommand
+public class HMMSearch extends Search
 {
-  static final String HMMSEARCH = "hmmsearch";
 
   boolean realign = false;
 
   boolean trim = false;
 
+  boolean returnNoOfNewSeqs = false;
+
   int seqsToReturn = Integer.MAX_VALUE;
 
-  SequenceI[] seqs;
 
   /**
    * Constructor for the HMMSearchThread
@@ -57,32 +60,34 @@ public class HMMSearch extends HmmerCommand
   @Override
   public void run()
   {
-    HiddenMarkovModel hmm = af.getSelectedHMM();
+    HiddenMarkovModel hmm = getHmmProfile();
     if (hmm == null)
     {
-      JOptionPane.showMessageDialog(af,
-              MessageManager.getString("warn.no_selected_hmm"));
+      // shouldn't happen if we got this far
+      Cache.log.error("Error: no hmm for hmmsearch");
       return;
     }
 
-    SequenceI hmmSeq = af.getSelectedHMMSequence();
+    SequenceI hmmSeq = hmm.getConsensusSequence();
     long msgId = System.currentTimeMillis();
-    af.setProgressBar(MessageManager.getString("status.running_hmmsearch"),
+    af.setProgressBar(MessageManager.getString("status.running_search"),
             msgId);
 
     try
     {
-      File hmmFile = createTempFile("hmm", ".hmm");
-      File hitsAlignmentFile = createTempFile("hitAlignment", ".sto");
-      File searchOutputFile = createTempFile("searchOutput", ".sto");
+      File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
+      File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
+              ".sto");
+      File searchOutputFile = FileUtils.createTempFile("searchOutput",
+              ".sto");
 
       exportHmm(hmm, hmmFile.getAbsoluteFile());
 
       boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
       if (!ran)
       {
-        JvOptionPane.showInternalMessageDialog(af,
-                MessageManager.getString("warn.hmmsearch_failed"));
+        JvOptionPane.showInternalMessageDialog(af, MessageManager
+                .formatMessage("warn.command_failed", "hmmsearch"));
         return;
       }
 
@@ -122,102 +127,15 @@ public class HMMSearch extends HmmerCommand
 
     List<String> args = new ArrayList<>();
     args.add(command);
-    args.add("-o");
-    args.add(searchOutputFile.getAbsolutePath());
-    args.add("-A");
-    args.add(hitsAlignmentFile.getAbsolutePath());
-
-    boolean dbFound = false;
-    String dbPath = "";
-    File databaseFile = null;
-
-    if (params != null)
-    {
-      for (ArgumentI arg : params)
-      {
-        String name = arg.getName();
-        if (MessageManager.getString("label.number_of_results")
-                .equals(name))
-        {
-          seqsToReturn = Integer.parseInt(arg.getValue());
-        }
-        else if (MessageManager.getString("label.auto_align_seqs")
-                .equals(name))
-        {
-          realign = true; // TODO: not used
-        }
-        else if (MessageManager.getString("label.use_accessions")
-                .equals(name))
-        {
-          args.add("--acc");
-        }
-        else if (MessageManager.getString("label.seq_e_value").equals(name))
-        {
-          args.add("--incE");
-          args.add(arg.getValue());
-        }
-        else if (MessageManager.getString("label.seq_score").equals(name))
-        {
-          args.add("-incT");
-          args.add(arg.getValue());
-        }
-        else if (MessageManager.getString("label.dom_e_value_desc")
-                .equals(name))
-        {
-          args.add("--incdomE");
-          args.add(arg.getValue());
-        }
-        else if (MessageManager.getString("label.dom_score").equals(name))
-        {
-          args.add("--incdomT");
-          args.add(arg.getValue());
-        }
-        else if (MessageManager.getString("label.trim_termini")
-                .equals(name))
-        {
-          trim = true;
-        }
-        else if (MessageManager.getString("label.database").equals(name))
-        {
-          dbFound = true;
-          dbPath = arg.getValue();
-          if (!MessageManager.getString("label.this_alignment")
-                  .equals(dbPath))
-          {
-            databaseFile = new File(dbPath);
-          }
-        }
-      }
-    }
-
-    if (!dbFound || MessageManager.getString("label.this_alignment")
-            .equals(dbPath))
-    {
-      AlignmentI alignment = af.getViewport().getAlignment();
-      AlignmentI copy = new Alignment(alignment);
-      SequenceI hmms = copy.getHmmConsensus();
-      if (hmms != null)
-      {
-        copy.deleteSequence(hmms);
-      }
-      StockholmFile stoFile = new StockholmFile(copy);
-      stoFile.setSeqs(copy.getSequencesArray());
-      String alignmentString = stoFile.print();
-      databaseFile = createTempFile("database", ".sto");
-      PrintWriter writer = new PrintWriter(databaseFile);
-      writer.print(alignmentString);
-      writer.close();
-    }
-
-    args.add(hmmFile.getAbsolutePath());
-    args.add(databaseFile.getAbsolutePath());
+    buildArguments(args, searchOutputFile, hitsAlignmentFile, hmmFile);
 
     return runCommand(args);
   }
 
+
   /**
    * Imports the data from the temporary file to which the output of hmmsearch
-   * is directed.
+   * was directed. The results are optionally realigned using hmmalign.
    * 
    * @param hmmSeq
    */
@@ -239,27 +157,95 @@ public class HMMSearch extends HmmerCommand
               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
       seqs = file.getSeqsAsArray();
 
-      readTable(searchOutputFile);
+      readDomainTable(searchOutputFile, false);
+
+      if (searchAlignment)
+      {
+        recoverSequences(sequencesHash, seqs);
+      }
+
+      // look for PP cons and ref seq in alignment only annotation
+      AlignmentAnnotation modelpos = null, ppcons = null;
+      for (AlignmentAnnotation aa : file.getAnnotations())
+      {
+        if (aa.sequenceRef == null)
+        {
+          if (aa.label.equals("Reference Positions")) // RF feature type in
+                                                      // stockholm parser
+          {
+            modelpos = aa;
+          }
+          if (aa.label.equals("Posterior Probability"))
+          {
+            ppcons = aa;
+          }
+        }
+      }
+
 
       int seqCount = Math.min(seqs.length, seqsToReturn);
       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
+      hmmSeq = hmmSeq.deriveSequence(); // otherwise all bad things happen
       hmmAndSeqs[0] = hmmSeq;
       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
+      if (modelpos != null)
+      {
+        // TODO need - get ungapped sequence method
+        hmmSeq.setSequence(
+                hmmSeq.getDatasetSequence().getSequenceAsString());
+        Annotation[] refpos = modelpos.annotations;
+        // insert gaps to match with refseq positions
+        int gc = 0, lcol = 0;
+        for (int c = 0; c < refpos.length; c++)
+        {
+          if (refpos[c] != null && ("x".equals(refpos[c].displayCharacter)))
+          {
+            if (gc > 0)
+            {
+              hmmSeq.insertCharAt(lcol + 1, gc, '-');
+            }
+            gc = 0;
+            lcol = c;
+          }
+          else
+          {
+            gc++;
+          }
+        }
+      }
 
-      AlignmentI alignment = new Alignment(hmmAndSeqs);
-      AlignFrame frame = new AlignFrame(alignment, 1, 1);
-      frame.setSelectedHMMSequence(hmmSeq);
-      List<ArgumentI> alignArgs = new ArrayList<>();
-      if (trim)
+      if (realign)
       {
-        alignArgs.add(new BooleanOption(
-                MessageManager.getString("label.trim_termini"),
-                MessageManager.getString("label.trim_termini_desc"), true,
-                true, true, null));
+        realignResults(hmmAndSeqs);
       }
-      HMMAlign hmmalign = new HMMAlign(frame, alignArgs);
-      hmmalign.run();
-      frame = null;
+      else
+      {
+        AlignmentI al = new Alignment(hmmAndSeqs);
+        if (ppcons != null)
+        {
+          al.addAnnotation(ppcons);
+        }
+        if (modelpos != null)
+        {
+          al.addAnnotation(modelpos);
+        }
+        AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
+                AlignFrame.DEFAULT_HEIGHT);
+        String ttl = "hmmSearch of " + databaseName + " using "
+                + hmmSeq.getName();
+        Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
+                AlignFrame.DEFAULT_HEIGHT);
+
+        if (returnNoOfNewSeqs)
+        {
+          int nNew = checkForNewSequences();
+          JvOptionPane.showMessageDialog(af.alignPanel, nNew + " "
+                  + MessageManager.getString("label.new_returned"));
+        }
+
+      }
+
+
       hmmTemp.delete();
       inputAlignmentTemp.delete();
       searchOutputFile.delete();
@@ -272,45 +258,61 @@ public class HMMSearch extends HmmerCommand
     }
   }
 
-  void readTable(File inputTableTemp) throws IOException
+  private int checkForNewSequences()
   {
-    BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
-    String line = "";
-    while (!line.startsWith("Query:"))
+    int nNew = seqs.length;
+
+    for (SequenceI resultSeq : seqs)
     {
-      line = br.readLine();
+      for (SequenceI aliSeq : alignment.getSequencesArray())
+      {
+        if (resultSeq.getName().equals(aliSeq.getName()))
+        {
+          nNew--;
+          break;
+        }
+      }
     }
-    for (int i = 0; i < 5; i++)
+
+    return nNew;
+
+  }
+
+  /**
+   * Realigns the given sequences using hmmalign, to the HMM profile sequence
+   * which is the first in the array, and opens the results in a new frame
+   * 
+   * @param hmmAndSeqs
+   */
+  protected void realignResults(SequenceI[] hmmAndSeqs)
+  {
+    /*
+     * and align the search results to the HMM profile
+     */
+    AlignmentI al = new Alignment(hmmAndSeqs);
+    AlignFrame frame = new AlignFrame(al, 1, 1);
+    List<ArgumentI> alignArgs = new ArrayList<>();
+    String alignTo = hmmAndSeqs[0].getName();
+    List<String> options = Collections.singletonList(alignTo);
+    Option option = new Option(MessageManager.getString("label.use_hmm"),
+            "", true, alignTo, alignTo, options, null);
+    alignArgs.add(option);
+    if (trim)
     {
-      line = br.readLine();
+      alignArgs.add(new BooleanOption(
+              MessageManager.getString(TRIM_TERMINI_KEY),
+              MessageManager.getString("label.trim_termini_desc"), true,
+              true, true, null));
     }
+    HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
+    hmmalign.run();
 
-    int index = 0;
-    while (!"  ------ inclusion threshold ------".equals(line)
-            && !"".equals(line))
+    if (returnNoOfNewSeqs)
     {
-      Scanner scanner = new Scanner(line);
-
-      String str = scanner.next(); // full sequence eValue score
-      float eValue = Float.parseFloat(str);
-      int seqLength = seqs[index].getLength();
-      Annotation[] annots = new Annotation[seqLength];
-      for (int j = 0; j < seqLength; j++)
-      {
-        annots[j] = new Annotation(eValue);
-      }
-      AlignmentAnnotation annot = new AlignmentAnnotation("E-value",
-              "Score", annots);
-      annot.setScore(Double.parseDouble(str));
-      annot.setSequenceRef(seqs[index]);
-      seqs[index].addAlignmentAnnotation(annot);
-
-      scanner.close();
-      line = br.readLine();
-      index++;
+      int nNew = checkForNewSequences();
+      JvOptionPane.showMessageDialog(frame.alignPanel,
+              nNew + " " + MessageManager.getString("label.new_returned"));
     }
-
-    br.close();
   }
 
 }