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