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     boolean searchAlignment = false;
 199
 200     if (params != null)
 201     {
 202       for (ArgumentI arg : params)
 203       {
 204         String name = arg.getName();
 205         if (MessageManager.getString(NUMBER_OF_RESULTS_KEY)
 206                 .equals(name))
 207         {
 208           seqsToReturn = Integer.parseInt(arg.getValue());
 209         }
 210         else if (MessageManager.getString("action.search").equals(name))
 211         {
 212           searchAlignment = arg.getValue().equals(
 213                   MessageManager.getString(HMMSearch.THIS_ALIGNMENT_KEY));
 214         }
 215         else if (MessageManager.getString(DATABASE_KEY).equals(name))
 216         {
 217           dbPath = arg.getValue();
 218           int pos = dbPath.lastIndexOf(File.separator);
 219           databaseName = dbPath.substring(pos + 1);
 220           databaseFile = new File(dbPath);
 221         }
 222         else if (MessageManager.getString(AUTO_ALIGN_SEQS_KEY)
 223                 .equals(name))
 224         {
 225           realign = true;
 226         }
 227         else if (MessageManager.getString(USE_ACCESSIONS_KEY)
 228                 .equals(name))
 229         {
 230           args.add("--acc");
 231         }
 232         else if (MessageManager.getString(REPORTING_CUTOFF_KEY)
 233                 .equals(name))
 234         {
 235           if (CUTOFF_EVALUE.equals(arg.getValue()))
 236           {
 237             useEvalueCutoff = true;
 238           }
 239           else if (CUTOFF_SCORE.equals(arg.getValue()))
 240           {
 241             useScoreCutoff = true;
 242           }
 243         }
 244         else if (MessageManager.getString(SEQ_EVALUE_KEY).equals(name))
 245         {
 246           seqEvalueCutoff = arg.getValue();
 247         }
 248         else if (MessageManager.getString(SEQ_SCORE_KEY).equals(name))
 249         {
 250           seqScoreCutoff = arg.getValue();
 251         }
 252         else if (MessageManager.getString(DOM_EVALUE_KEY)
 253                 .equals(name))
 254         {
 255           domEvalueCutoff = arg.getValue();
 256         }
 257         else if (MessageManager.getString(DOM_SCORE_KEY).equals(name))
 258         {
 259           domScoreCutoff = arg.getValue();
 260         }
 261         else if (MessageManager.getString(TRIM_TERMINI_KEY)
 262                 .equals(name))
 263         {
 264           trim = true;
 265         }
 266         else if (MessageManager.getString(DATABASE_KEY).equals(name))
 267         {
 268           dbFound = true;
 269           dbPath = arg.getValue();
 270           if (!MessageManager.getString(THIS_ALIGNMENT_KEY)
 271                   .equals(dbPath))
 272           {
 273             int pos = dbPath.lastIndexOf(File.separator);
 274             databaseName = dbPath.substring(pos + 1);
 275             databaseFile = new File(dbPath);
 276           }
 277         }
 278       }
 279     }
 280
 281     if (useEvalueCutoff)
 282     {
 283       args.add("-E");
 284       args.add(seqEvalueCutoff);
 285       args.add("--domE");
 286       args.add(domEvalueCutoff);
 287     }
 288     else if (useScoreCutoff)
 289     {
 290       args.add("-T");
 291       args.add(seqScoreCutoff);
 292       args.add("--domT");
 293       args.add(domScoreCutoff);
 294     }
 295
 296 //    if (!dbFound || MessageManager.getString(THIS_ALIGNMENT_KEY)
 297 //            .equals(dbPath))
 298       if (searchAlignment)
 299     {
 300       /*
 301        * no external database specified for search, so
 302        * export current alignment as 'database' to search,
 303        * excluding any HMM consensus sequences it contains
 304        */
 305       databaseFile = FileUtils.createTempFile("database", ".sto");
 306       AlignmentI al = af.getViewport().getAlignment();
 307       AlignmentI copy = new Alignment(al);
 308       List<SequenceI> hmms = copy.getHmmSequences();
 309       for (SequenceI hmmSeq : hmms)
 310       {
 311         copy.deleteSequence(hmmSeq);
 312       }
 313       exportStockholm(copy.getSequencesArray(), databaseFile, null);
 314     }
 315
 316     args.add(getFilePath(hmmFile));
 317     args.add(getFilePath(databaseFile));
 318   }
 319
 320   /**
 321    * Imports the data from the temporary file to which the output of hmmsearch
 322    * was directed. The results are optionally realigned using hmmalign.
 323    *
 324    * @param hmmSeq
 325    */
 326   private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
 327           File hmmTemp, File searchOutputFile)
 328           throws IOException, InterruptedException
 329   {
 330     BufferedReader br = new BufferedReader(
 331             new FileReader(inputAlignmentTemp));
 332     try
 333     {
 334       if (br.readLine() == null)
 335       {
 336         JOptionPane.showMessageDialog(af,
 337                 MessageManager.getString("label.no_sequences_found"));
 338         return;
 339       }
 340       StockholmFile file = new StockholmFile(new FileParse(
 341               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
 342       seqs = file.getSeqsAsArray();
 343
 344       readTable(searchOutputFile);
 345
 346       int seqCount = Math.min(seqs.length, seqsToReturn);
 347       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 348       hmmAndSeqs[0] = hmmSeq;
 349       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 350
 351       if (realign)
 352       {
 353         realignResults(hmmAndSeqs);
 354       }
 355       else
 356       {
 357         AlignmentI al = new Alignment(hmmAndSeqs);
 358         AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 359                 AlignFrame.DEFAULT_HEIGHT);
 360         String ttl = "hmmSearch of " + databaseName + " using "
 361                 + hmmSeq.getName();
 362         Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 363                 AlignFrame.DEFAULT_HEIGHT);
 364       }
 365
 366       hmmTemp.delete();
 367       inputAlignmentTemp.delete();
 368       searchOutputFile.delete();
 369     } finally
 370     {
 371       if (br != null)
 372       {
 373         br.close();
 374       }
 375     }
 376   }
 377
 378   /**
 379    * Realigns the given sequences using hmmalign, to the HMM profile sequence
 380    * which is the first in the array, and opens the results in a new frame
 381    *
 382    * @param hmmAndSeqs
 383    */
 384   protected void realignResults(SequenceI[] hmmAndSeqs)
 385   {
 386     /*
 387      * and align the search results to the HMM profile
 388      */
 389     AlignmentI al = new Alignment(hmmAndSeqs);
 390     AlignFrame frame = new AlignFrame(al, 1, 1);
 391     List<ArgumentI> alignArgs = new ArrayList<>();
 392     String alignTo = hmmAndSeqs[0].getName();
 393     List<String> options = Collections.singletonList(alignTo);
 394     Option option = new Option(MessageManager.getString("label.use_hmm"),
 395             "", true, alignTo, alignTo, options, null);
 396     alignArgs.add(option);
 397     if (trim)
 398     {
 399       alignArgs.add(new BooleanOption(
 400               MessageManager.getString(TRIM_TERMINI_KEY),
 401               MessageManager.getString("label.trim_termini_desc"), true,
 402               true, true, null));
 403     }
 404     HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 405     hmmalign.run();
 406   }
 407
 408   /**
 409    * Reads in the scores table output by hmmsearch and adds annotation to
 410    * sequences for E-value and bit score
 411    *
 412    * @param inputTableTemp
 413    * @throws IOException
 414    */
 415   void readTable(File inputTableTemp) throws IOException
 416   {
 417     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 418     String line = "";
 419     while (!line.startsWith("Query:"))
 420     {
 421       line = br.readLine();
 422     }
 423     for (int i = 0; i < 5; i++)
 424     {
 425       line = br.readLine();
 426     }
 427
 428     int index = 0;
 429     while (!"  ------ inclusion threshold ------".equals(line)
 430             && !"".equals(line))
 431     {
 432       SequenceI seq = seqs[index];
 433       Scanner scanner = new Scanner(line);
 434       String str = scanner.next();
 435       addScoreAnnotation(str, seq, "hmmsearch E-value",
 436               "Full sequence E-value");
 437       str = scanner.next();
 438       addScoreAnnotation(str, seq, "hmmsearch Score",
 439               "Full sequence bit score");
 440       scanner.close();
 441       line = br.readLine();
 442       index++;
 443     }
 444
 445     br.close();
 446   }
 447
 448   /**
 449    * A helper method that adds one score-only (non-positional) annotation to a
 450    * sequence
 451    *
 452    * @param value
 453    * @param seq
 454    * @param label
 455    * @param description
 456    */
 457   protected void addScoreAnnotation(String value, SequenceI seq,
 458           String label, String description)
 459   {
 460     try
 461     {
 462       AlignmentAnnotation annot = new AlignmentAnnotation(label,
 463               description, null);
 464       annot.setCalcId(HMMSEARCH);
 465       double eValue = Double.parseDouble(value);
 466       annot.setScore(eValue);
 467       annot.setSequenceRef(seq);
 468       seq.addAlignmentAnnotation(annot);
 469     } catch (NumberFormatException e)
 470     {
 471       System.err.println("Error parsing " + label + " from " + value);
 472     }
 473   }
 474
 475 }