src/jalview/hmmer/HMMSearch.java

   1 package jalview.hmmer;
   2
   3 import jalview.bin.Cache;
   4 import jalview.datamodel.Alignment;
   5 import jalview.datamodel.AlignmentAnnotation;
   6 import jalview.datamodel.AlignmentI;
   7 import jalview.datamodel.HiddenMarkovModel;
   8 import jalview.datamodel.SequenceI;
   9 import jalview.gui.AlignFrame;
  10 import jalview.gui.Desktop;
  11 import jalview.gui.JvOptionPane;
  12 import jalview.io.DataSourceType;
  13 import jalview.io.FileParse;
  14 import jalview.io.StockholmFile;
  15 import jalview.util.FileUtils;
  16 import jalview.util.MessageManager;
  17 import jalview.ws.params.ArgumentI;
  18 import jalview.ws.params.simple.BooleanOption;
  19 import jalview.ws.params.simple.Option;
  20
  21 import java.io.BufferedReader;
  22 import java.io.File;
  23 import java.io.FileReader;
  24 import java.io.IOException;
  25 import java.util.ArrayList;
  26 import java.util.Collections;
  27 import java.util.List;
  28 import java.util.Scanner;
  29
  30 import javax.swing.JOptionPane;
  31
  32 public class HMMSearch extends HmmerCommand
  33 {
  34   static final String HMMSEARCH = "hmmsearch";
  35
  36   /*
  37    * constants for i18n lookup of passed parameter names
  38    */
  39   static final String DATABASE_KEY = "label.database";
  40
  41   static final String THIS_ALIGNMENT_KEY = "label.this_alignment";
  42
  43   static final String USE_ACCESSIONS_KEY = "label.use_accessions";
  44
  45   static final String AUTO_ALIGN_SEQS_KEY = "label.auto_align_seqs";
  46
  47   static final String NUMBER_OF_RESULTS_KEY = "label.number_of_results";
  48
  49   static final String TRIM_TERMINI_KEY = "label.trim_termini";
  50
  51   static final String REPORTING_CUTOFF_KEY = "label.reporting_cutoff";
  52
  53   static final String CUTOFF_NONE = "None";
  54
  55   static final String CUTOFF_SCORE = "Score";
  56
  57   static final String CUTOFF_EVALUE = "E-Value";
  58
  59   static final String SEQ_EVALUE_KEY = "label.seq_evalue";
  60
  61   static final String DOM_EVALUE_KEY = "label.dom_evalue";
  62
  63   static final String SEQ_SCORE_KEY = "label.seq_score";
  64
  65   static final String DOM_SCORE_KEY = "label.dom_score";
  66
  67   boolean realign = false;
  68
  69   boolean trim = false;
  70
  71   int seqsToReturn = Integer.MAX_VALUE;
  72
  73   SequenceI[] seqs;
  74
  75   private String databaseName;
  76
  77   /**
  78    * Constructor for the HMMSearchThread
  79    *
  80    * @param af
  81    */
  82   public HMMSearch(AlignFrame af, List<ArgumentI> args)
  83   {
  84     super(af, args);
  85   }
  86
  87   /**
  88    * Runs the HMMSearchThread: the data on the alignment or group is exported,
  89    * then the command is executed in the command line and then the data is
  90    * imported and displayed in a new frame. Call this method directly to execute
  91    * synchronously, or via start() in a new Thread for asynchronously.
  92    */
  93   @Override
  94   public void run()
  95   {
  96     HiddenMarkovModel hmm = getHmmProfile();
  97     if (hmm == null)
  98     {
  99       // shouldn't happen if we got this far
 100       Cache.log.error("Error: no hmm for hmmsearch");
 101       return;
 102     }
 103
 104     SequenceI hmmSeq = hmm.getConsensusSequence();
 105     long msgId = System.currentTimeMillis();
 106     af.setProgressBar(MessageManager.getString("status.running_hmmsearch"),
 107             msgId);
 108
 109     try
 110     {
 111       File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
 112       File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
 113               ".sto");
 114       File searchOutputFile = FileUtils.createTempFile("searchOutput",
 115               ".sto");
 116
 117       exportHmm(hmm, hmmFile.getAbsoluteFile());
 118
 119       boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
 120       if (!ran)
 121       {
 122         JvOptionPane.showInternalMessageDialog(af, MessageManager
 123                 .formatMessage("warn.command_failed", "hmmsearch"));
 124         return;
 125       }
 126
 127       importData(hmmSeq, hitsAlignmentFile, hmmFile, searchOutputFile);
 128       // TODO make realignment of search results a step at this level
 129       // and make it conditional on this.realign
 130     } catch (IOException | InterruptedException e)
 131     {
 132       e.printStackTrace();
 133     }
 134     finally
 135     {
 136       af.setProgressBar("", msgId);
 137     }
 138   }
 139
 140   /**
 141    * Executes an hmmsearch with the given hmm as input. The database to be
 142    * searched is a local file as specified by the 'Database' parameter, or the
 143    * current alignment (written to file) if none is specified.
 144    *
 145    * @param searchOutputFile
 146    * @param hitsAlignmentFile
 147    * @param hmmFile
 148    *
 149    * @return
 150    * @throws IOException
 151    */
 152   private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
 153           File hmmFile) throws IOException
 154   {
 155     String command = getCommandPath(HMMSEARCH);
 156     if (command == null)
 157     {
 158       return false;
 159     }
 160
 161     List<String> args = new ArrayList<>();
 162     args.add(command);
 163     buildArguments(args, searchOutputFile, hitsAlignmentFile, hmmFile);
 164
 165     return runCommand(args);
 166   }
 167
 168   /**
 169    * Appends command line arguments to the given list, to specify input and
 170    * output files for the search, and any additional options that may have been
 171    * passed from the parameters dialog
 172    *
 173    * @param args
 174    * @param searchOutputFile
 175    * @param hitsAlignmentFile
 176    * @param hmmFile
 177    * @throws IOException
 178    */
 179   protected void buildArguments(List<String> args, File searchOutputFile,
 180           File hitsAlignmentFile, File hmmFile) throws IOException
 181   {
 182     args.add("-o");
 183     args.add(getFilePath(searchOutputFile));
 184     args.add("-A");
 185     args.add(getFilePath(hitsAlignmentFile));
 186
 187     boolean dbFound = false;
 188     String dbPath = "";
 189     File databaseFile = null;
 190
 191     boolean useEvalueCutoff = false;
 192     boolean useScoreCutoff = false;
 193     String seqEvalueCutoff = null;
 194     String domEvalueCutoff = null;
 195     String seqScoreCutoff = null;
 196     String domScoreCutoff = null;
 197     databaseName = "Alignment";
 198
 199     if (params != null)
 200     {
 201       for (ArgumentI arg : params)
 202       {
 203         String name = arg.getName();
 204         if (MessageManager.getString(NUMBER_OF_RESULTS_KEY)
 205                 .equals(name))
 206         {
 207           seqsToReturn = Integer.parseInt(arg.getValue());
 208         }
 209         else if (MessageManager.getString(AUTO_ALIGN_SEQS_KEY)
 210                 .equals(name))
 211         {
 212           realign = true;
 213         }
 214         else if (MessageManager.getString(USE_ACCESSIONS_KEY)
 215                 .equals(name))
 216         {
 217           args.add("--acc");
 218         }
 219         else if (MessageManager.getString(REPORTING_CUTOFF_KEY)
 220                 .equals(name))
 221         {
 222           if (CUTOFF_EVALUE.equals(arg.getValue()))
 223           {
 224             useEvalueCutoff = true;
 225           }
 226           else if (CUTOFF_SCORE.equals(arg.getValue()))
 227           {
 228             useScoreCutoff = true;
 229           }
 230         }
 231         else if (MessageManager.getString(SEQ_EVALUE_KEY).equals(name))
 232         {
 233           seqEvalueCutoff = arg.getValue();
 234         }
 235         else if (MessageManager.getString(SEQ_SCORE_KEY).equals(name))
 236         {
 237           seqScoreCutoff = arg.getValue();
 238         }
 239         else if (MessageManager.getString(DOM_EVALUE_KEY)
 240                 .equals(name))
 241         {
 242           domEvalueCutoff = arg.getValue();
 243         }
 244         else if (MessageManager.getString(DOM_SCORE_KEY).equals(name))
 245         {
 246           domScoreCutoff = arg.getValue();
 247         }
 248         else if (MessageManager.getString(TRIM_TERMINI_KEY)
 249                 .equals(name))
 250         {
 251           trim = true;
 252         }
 253         else if (MessageManager.getString(DATABASE_KEY).equals(name))
 254         {
 255           dbFound = true;
 256           dbPath = arg.getValue();
 257           if (!MessageManager.getString(THIS_ALIGNMENT_KEY)
 258                   .equals(dbPath))
 259           {
 260             int pos = dbPath.lastIndexOf(File.separator);
 261             databaseName = dbPath.substring(pos + 1);
 262             databaseFile = new File(dbPath);
 263           }
 264         }
 265       }
 266     }
 267
 268     if (useEvalueCutoff)
 269     {
 270       args.add("-E");
 271       args.add(seqEvalueCutoff);
 272       args.add("--domE");
 273       args.add(domEvalueCutoff);
 274     }
 275     else if (useScoreCutoff)
 276     {
 277       args.add("-T");
 278       args.add(seqScoreCutoff);
 279       args.add("--domT");
 280       args.add(domScoreCutoff);
 281     }
 282
 283     if (!dbFound || MessageManager.getString(THIS_ALIGNMENT_KEY)
 284             .equals(dbPath))
 285     {
 286       /*
 287        * no external database specified for search, so
 288        * export current alignment as 'database' to search,
 289        * excluding any HMM consensus sequences it contains
 290        */
 291       databaseFile = FileUtils.createTempFile("database", ".sto");
 292       AlignmentI al = af.getViewport().getAlignment();
 293       AlignmentI copy = new Alignment(al);
 294       List<SequenceI> hmms = copy.getHmmSequences();
 295       for (SequenceI hmmSeq : hmms)
 296       {
 297         copy.deleteSequence(hmmSeq);
 298       }
 299       exportStockholm(copy.getSequencesArray(), databaseFile, null);
 300     }
 301
 302     args.add(getFilePath(hmmFile));
 303     args.add(getFilePath(databaseFile));
 304   }
 305
 306   /**
 307    * Imports the data from the temporary file to which the output of hmmsearch
 308    * was directed. The results are optionally realigned using hmmalign.
 309    *
 310    * @param hmmSeq
 311    */
 312   private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
 313           File hmmTemp, File searchOutputFile)
 314           throws IOException, InterruptedException
 315   {
 316     BufferedReader br = new BufferedReader(
 317             new FileReader(inputAlignmentTemp));
 318     try
 319     {
 320       if (br.readLine() == null)
 321       {
 322         JOptionPane.showMessageDialog(af,
 323                 MessageManager.getString("label.no_sequences_found"));
 324         return;
 325       }
 326       StockholmFile file = new StockholmFile(new FileParse(
 327               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
 328       seqs = file.getSeqsAsArray();
 329
 330       readTable(searchOutputFile);
 331
 332       int seqCount = Math.min(seqs.length, seqsToReturn);
 333       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 334       hmmAndSeqs[0] = hmmSeq;
 335       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 336
 337       if (realign)
 338       {
 339         realignResults(hmmAndSeqs);
 340       }
 341       else
 342       {
 343         AlignmentI al = new Alignment(hmmAndSeqs);
 344         AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 345                 AlignFrame.DEFAULT_HEIGHT);
 346         String ttl = "hmmSearch of " + databaseName + " using "
 347                 + hmmSeq.getName();
 348         Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 349                 AlignFrame.DEFAULT_HEIGHT);
 350       }
 351
 352       hmmTemp.delete();
 353       inputAlignmentTemp.delete();
 354       searchOutputFile.delete();
 355     } finally
 356     {
 357       if (br != null)
 358       {
 359         br.close();
 360       }
 361     }
 362   }
 363
 364   /**
 365    * Realigns the given sequences using hmmalign, to the HMM profile sequence
 366    * which is the first in the array, and opens the results in a new frame
 367    *
 368    * @param hmmAndSeqs
 369    */
 370   protected void realignResults(SequenceI[] hmmAndSeqs)
 371   {
 372     /*
 373      * and align the search results to the HMM profile
 374      */
 375     AlignmentI al = new Alignment(hmmAndSeqs);
 376     AlignFrame frame = new AlignFrame(al, 1, 1);
 377     List<ArgumentI> alignArgs = new ArrayList<>();
 378     String alignTo = hmmAndSeqs[0].getName();
 379     List<String> options = Collections.singletonList(alignTo);
 380     Option option = new Option(MessageManager.getString("label.use_hmm"),
 381             "", true, alignTo, alignTo, options, null);
 382     alignArgs.add(option);
 383     if (trim)
 384     {
 385       alignArgs.add(new BooleanOption(
 386               MessageManager.getString(TRIM_TERMINI_KEY),
 387               MessageManager.getString("label.trim_termini_desc"), true,
 388               true, true, null));
 389     }
 390     HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 391     hmmalign.run();
 392   }
 393
 394   /**
 395    * Reads in the scores table output by hmmsearch and adds annotation to
 396    * sequences for E-value and bit score
 397    *
 398    * @param inputTableTemp
 399    * @throws IOException
 400    */
 401   void readTable(File inputTableTemp) throws IOException
 402   {
 403     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 404     String line = "";
 405     while (!line.startsWith("Query:"))
 406     {
 407       line = br.readLine();
 408     }
 409     for (int i = 0; i < 5; i++)
 410     {
 411       line = br.readLine();
 412     }
 413
 414     int index = 0;
 415     while (!"  ------ inclusion threshold ------".equals(line)
 416             && !"".equals(line))
 417     {
 418       SequenceI seq = seqs[index];
 419       Scanner scanner = new Scanner(line);
 420       String str = scanner.next();
 421       addScoreAnnotation(str, seq, "hmmsearch E-value",
 422               "Full sequence E-value");
 423       str = scanner.next();
 424       addScoreAnnotation(str, seq, "hmmsearch Score",
 425               "Full sequence bit score");
 426       scanner.close();
 427       line = br.readLine();
 428       index++;
 429     }
 430
 431     br.close();
 432   }
 433
 434   /**
 435    * A helper method that adds one score-only (non-positional) annotation to a
 436    * sequence
 437    *
 438    * @param value
 439    * @param seq
 440    * @param label
 441    * @param description
 442    */
 443   protected void addScoreAnnotation(String value, SequenceI seq,
 444           String label, String description)
 445   {
 446     try
 447     {
 448       AlignmentAnnotation annot = new AlignmentAnnotation(label,
 449               description, null);
 450       annot.setCalcId(HMMSEARCH);
 451       double eValue = Double.parseDouble(value);
 452       annot.setScore(eValue);
 453       annot.setSequenceRef(seq);
 454       seq.addAlignmentAnnotation(annot);
 455     } catch (NumberFormatException e)
 456     {
 457       System.err.println("Error parsing " + label + " from " + value);
 458     }
 459   }
 460
 461 }