JAL-2629 hmmer searches now read domain rather than full scores
[jalview.git] / src / jalview / hmmer / HMMSearch.java
index 7413d14..f05823e 100644 (file)
@@ -8,6 +8,7 @@ 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;
@@ -25,23 +26,20 @@ import java.io.IOException;
 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
 {
-  private static final String PARAMNAME_NO_OF_RESULTS = MessageManager.getString("label.number_of_results");
-
-  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
@@ -70,9 +68,9 @@ public class HMMSearch extends HmmerCommand
       return;
     }
 
-    SequenceI hmmSeq = hmm.getConsensusSequence();// 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
@@ -129,108 +127,15 @@ public class HMMSearch extends HmmerCommand
 
     List<String> args = new ArrayList<>();
     args.add(command);
-    args.add("-o");
-    args.add(getFilePath(searchOutputFile));
-    args.add("-A");
-    args.add(getFilePath(hitsAlignmentFile));
-
-    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))
-    {
-      /*
-       * no external database specified for search, so
-       * export current alignment as 'database' to search,
-       * excluding any HMM consensus sequences it contains
-       */
-      databaseFile = FileUtils.createTempFile("database", ".sto");
-      AlignmentI al = af.getViewport().getAlignment();
-      AlignmentI copy = new Alignment(al);
-      List<SequenceI> hmms = copy.getHmmSequences();
-      for (SequenceI hmmSeq : hmms)
-      {
-        copy.deleteSequence(hmmSeq);
-      }
-      exportStockholm(copy.getSequencesArray(), databaseFile, null);
-      // StockholmFile stoFile = new StockholmFile(copy);
-      // stoFile.setSeqs(copy.getSequencesArray());
-      // String alignmentString = stoFile.print();
-      // PrintWriter writer = new PrintWriter(databaseFile);
-      // writer.print(alignmentString);
-      // writer.close();
-    }
-
-    args.add(getFilePath(hmmFile));
-    args.add(getFilePath(databaseFile));
+    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
    */
@@ -252,34 +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++;
+          }
+        }
+      }
 
-      /*
-       * 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 defSeq = hmmSeq.getName();
-      List<String> options = Collections.singletonList(defSeq);
-      Option option = new Option(MessageManager.getString("label.use_hmm"),
-              "", true, defSeq, defSeq, options, null);
-      alignArgs.add(option);
-      if (trim)
+      if (realign)
+      {
+        realignResults(hmmAndSeqs);
+      }
+      else
       {
-        alignArgs.add(new BooleanOption(
-                MessageManager.getString("label.trim_termini"),
-                MessageManager.getString("label.trim_termini_desc"), true,
-                true, true, null));
+        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"));
+        }
+
       }
-      HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
-      hmmalign.run();
-      frame = null;
+
+
       hmmTemp.delete();
       inputAlignmentTemp.delete();
       searchOutputFile.delete();
@@ -292,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();
   }
 
 }