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