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