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