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