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       deleteHmmSequences(copy);
 279       exportStockholm(copy.getSequencesArray(), databaseFile, null);
 280     }
 281
 282     args.add(getFilePath(hmmFile, true));
 283     args.add(getFilePath(databaseFile, true));
 284   }
 285
 286   /**
 287    * Imports the data from the temporary file to which the output of hmmsearch
 288    * was directed. The results are optionally realigned using hmmalign.
 289    *
 290    * @param hmmSeq
 291    */
 292   private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
 293           File hmmTemp, File searchOutputFile)
 294           throws IOException, InterruptedException
 295   {
 296     BufferedReader br = new BufferedReader(
 297             new FileReader(inputAlignmentTemp));
 298     try
 299     {
 300       if (br.readLine() == null)
 301       {
 302         JOptionPane.showMessageDialog(af,
 303                 MessageManager.getString("label.no_sequences_found"));
 304         return;
 305       }
 306       StockholmFile file = new StockholmFile(new FileParse(
 307               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
 308       seqs = file.getSeqsAsArray();
 309       // look for PP cons and ref seq in alignment only annotation
 310       AlignmentAnnotation modelpos = null, ppcons = null;
 311       for (AlignmentAnnotation aa : file.getAnnotations())
 312       {
 313         if (aa.sequenceRef == null)
 314         {
 315           if (aa.label.equals("Reference Positions")) // RF feature type in
 316                                                       // stockholm parser
 317           {
 318             modelpos = aa;
 319           }
 320           if (aa.label.equals("Posterior Probability"))
 321           {
 322             ppcons = aa;
 323           }
 324         }
 325       }
 326       readTable(searchOutputFile);
 327
 328       int seqCount = Math.min(seqs.length, seqsToReturn);
 329       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 330       hmmSeq = hmmSeq.deriveSequence(); // otherwise all bad things happen
 331       hmmAndSeqs[0] = hmmSeq;
 332       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 333       if (modelpos != null)
 334       {
 335         // TODO need - get ungapped sequence method
 336         hmmSeq.setSequence(
 337                 hmmSeq.getDatasetSequence().getSequenceAsString());
 338         Annotation[] refpos = modelpos.annotations;
 339         // insert gaps to match with refseq positions
 340         int gc = 0, lcol = 0;
 341         for (int c = 0; c < refpos.length; c++)
 342         {
 343           if (refpos[c] != null && ("x".equals(refpos[c].displayCharacter)))
 344           {
 345             if (gc > 0)
 346             {
 347               hmmSeq.insertCharAt(lcol + 1, gc, '-');
 348             }
 349             gc = 0;
 350             lcol = c;
 351           }
 352           else
 353           {
 354             gc++;
 355           }
 356         }
 357       }
 358       if (realign)
 359       {
 360         realignResults(hmmAndSeqs);
 361       }
 362       else
 363       {
 364         AlignmentI al = new Alignment(hmmAndSeqs);
 365         if (ppcons != null)
 366         {
 367           al.addAnnotation(ppcons);
 368         }
 369         if (modelpos != null)
 370         {
 371           al.addAnnotation(modelpos);
 372         }
 373         AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 374                 AlignFrame.DEFAULT_HEIGHT);
 375         String ttl = "hmmSearch of " + databaseName + " using "
 376                 + hmmSeq.getName();
 377         Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 378                 AlignFrame.DEFAULT_HEIGHT);
 379       }
 380
 381       hmmTemp.delete();
 382       inputAlignmentTemp.delete();
 383       searchOutputFile.delete();
 384     } finally
 385     {
 386       if (br != null)
 387       {
 388         br.close();
 389       }
 390     }
 391   }
 392
 393   /**
 394    * Realigns the given sequences using hmmalign, to the HMM profile sequence
 395    * which is the first in the array, and opens the results in a new frame
 396    *
 397    * @param hmmAndSeqs
 398    */
 399   protected void realignResults(SequenceI[] hmmAndSeqs)
 400   {
 401     /*
 402      * and align the search results to the HMM profile
 403      */
 404     AlignmentI al = new Alignment(hmmAndSeqs);
 405     AlignFrame frame = new AlignFrame(al, 1, 1);
 406     List<ArgumentI> alignArgs = new ArrayList<>();
 407     String alignTo = hmmAndSeqs[0].getName();
 408     List<String> options = Collections.singletonList(alignTo);
 409     Option option = new Option(MessageManager.getString("label.use_hmm"),
 410             "", true, alignTo, alignTo, options, null);
 411     alignArgs.add(option);
 412     if (trim)
 413     {
 414       alignArgs.add(new BooleanOption(
 415               MessageManager.getString(TRIM_TERMINI_KEY),
 416               MessageManager.getString("label.trim_termini_desc"), true,
 417               true, true, null));
 418     }
 419     HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 420     hmmalign.run();
 421   }
 422
 423   /**
 424    * Reads in the scores table output by hmmsearch and adds annotation to
 425    * sequences for E-value and bit score
 426    *
 427    * @param inputTableTemp
 428    * @throws IOException
 429    */
 430   void readTable(File inputTableTemp) throws IOException
 431   {
 432     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 433     String line = "";
 434     while (!line.startsWith("Query:"))
 435     {
 436       line = br.readLine();
 437     }
 438     while (!line.contains("-------"))
 439     {
 440       line = br.readLine();
 441     }
 442     line = br.readLine();
 443
 444     int index = 0;
 445     while (!"  ------ inclusion threshold ------".equals(line)
 446             && !"".equals(line))
 447     {
 448       SequenceI seq = seqs[index];
 449       Scanner scanner = new Scanner(line);
 450       String evalue = scanner.next();
 451       String score = scanner.next();
 452       addScoreAnnotations(evalue, score, seq);
 453       scanner.close();
 454       line = br.readLine();
 455       index++;
 456     }
 457
 458     br.close();
 459   }
 460
 461
 462   protected void addScoreAnnotations(String eValue, String bitScore,
 463           SequenceI seq)
 464   {
 465     String label = "Search Scores";
 466     String description = "Full sequence bit score and E-Value";
 467
 468     try
 469     {
 470       AlignmentAnnotation annot = new AlignmentAnnotation(label,
 471               description, null);
 472
 473       annot.label = label;
 474       annot.description = description;
 475
 476       annot.setCalcId(HMMSEARCH);
 477
 478       double dEValue = Double.parseDouble(eValue);
 479       annot.setEValue(dEValue);
 480
 481       double dBitScore = Double.parseDouble(bitScore);
 482       annot.setBitScore(dBitScore);
 483
 484       annot.setSequenceRef(seq);
 485       seq.addAlignmentAnnotation(annot);
 486     } catch (NumberFormatException e)
 487     {
 488       System.err.println("Error parsing " + label + " from " + eValue
 489               + " & " + bitScore);
 490     }
 491   }
 492
 493 }