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.HiddenMarkovModel;
   8 import jalview.datamodel.SequenceI;
   9 import jalview.gui.AlignFrame;
  10 import jalview.gui.Desktop;
  11 import jalview.gui.JvOptionPane;
  12 import jalview.io.DataSourceType;
  13 import jalview.io.FileParse;
  14 import jalview.io.StockholmFile;
  15 import jalview.util.FileUtils;
  16 import jalview.util.MessageManager;
  17 import jalview.ws.params.ArgumentI;
  18 import jalview.ws.params.simple.BooleanOption;
  19 import jalview.ws.params.simple.Option;
  20
  21 import java.io.BufferedReader;
  22 import java.io.File;
  23 import java.io.FileReader;
  24 import java.io.IOException;
  25 import java.util.ArrayList;
  26 import java.util.Collections;
  27 import java.util.List;
  28 import java.util.Scanner;
  29
  30 import javax.swing.JOptionPane;
  31
  32 public class HMMSearch extends HmmerCommand
  33 {
  34   static final String HMMSEARCH = "hmmsearch";
  35
  36   /*
  37    * constants for i18n lookup of passed parameter names
  38    */
  39   static final String DATABASE_KEY = "label.database";
  40
  41   static final String THIS_ALIGNMENT_KEY = "label.this_alignment";
  42
  43   static final String USE_ACCESSIONS_KEY = "label.use_accessions";
  44
  45   static final String AUTO_ALIGN_SEQS_KEY = "label.auto_align_seqs";
  46
  47   static final String NUMBER_OF_RESULTS_KEY = "label.number_of_results";
  48
  49   static final String TRIM_TERMINI_KEY = "label.trim_termini";
  50
  51   static final String REPORTING_CUTOFF_KEY = "label.reporting_cutoff";
  52
  53   static final String CUTOFF_NONE = "None";
  54
  55   static final String CUTOFF_SCORE = "Score";
  56
  57   static final String CUTOFF_EVALUE = "E-Value";
  58
  59   static final String SEQ_EVALUE_KEY = "label.seq_evalue";
  60
  61   static final String DOM_EVALUE_KEY = "label.dom_evalue";
  62
  63   static final String SEQ_SCORE_KEY = "label.seq_score";
  64
  65   static final String DOM_SCORE_KEY = "label.dom_score";
  66
  67   boolean realign = false;
  68
  69   boolean trim = false;
  70
  71   int seqsToReturn = Integer.MAX_VALUE;
  72
  73   SequenceI[] seqs;
  74
  75   private String databaseName;
  76
  77   /**
  78    * Constructor for the HMMSearchThread
  79    *
  80    * @param af
  81    */
  82   public HMMSearch(AlignFrame af, List<ArgumentI> args)
  83   {
  84     super(af, args);
  85   }
  86
  87   /**
  88    * Runs the HMMSearchThread: the data on the alignment or group is exported,
  89    * then the command is executed in the command line and then the data is
  90    * imported and displayed in a new frame. Call this method directly to execute
  91    * synchronously, or via start() in a new Thread for asynchronously.
  92    */
  93   @Override
  94   public void run()
  95   {
  96     HiddenMarkovModel hmm = getHmmProfile();
  97     if (hmm == null)
  98     {
  99       // shouldn't happen if we got this far
 100       Cache.log.error("Error: no hmm for hmmsearch");
 101       return;
 102     }
 103
 104     SequenceI hmmSeq = hmm.getConsensusSequence();
 105     long msgId = System.currentTimeMillis();
 106     af.setProgressBar(MessageManager.getString("status.running_hmmsearch"),
 107             msgId);
 108
 109     try
 110     {
 111       File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
 112       File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
 113               ".sto");
 114       File searchOutputFile = FileUtils.createTempFile("searchOutput",
 115               ".sto");
 116
 117       exportHmm(hmm, hmmFile.getAbsoluteFile());
 118
 119       boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
 120       if (!ran)
 121       {
 122         JvOptionPane.showInternalMessageDialog(af, MessageManager
 123                 .formatMessage("warn.command_failed", "hmmsearch"));
 124         return;
 125       }
 126
 127       importData(hmmSeq, hitsAlignmentFile, hmmFile, searchOutputFile);
 128       // TODO make realignment of search results a step at this level
 129       // and make it conditional on this.realign
 130     } catch (IOException | InterruptedException e)
 131     {
 132       e.printStackTrace();
 133     }
 134     finally
 135     {
 136       af.setProgressBar("", msgId);
 137     }
 138   }
 139
 140   /**
 141    * Executes an hmmsearch with the given hmm as input. The database to be
 142    * searched is a local file as specified by the 'Database' parameter, or the
 143    * current alignment (written to file) if none is specified.
 144    *
 145    * @param searchOutputFile
 146    * @param hitsAlignmentFile
 147    * @param hmmFile
 148    *
 149    * @return
 150    * @throws IOException
 151    */
 152   private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
 153           File hmmFile) throws IOException
 154   {
 155     String command = getCommandPath(HMMSEARCH);
 156     if (command == null)
 157     {
 158       return false;
 159     }
 160
 161     List<String> args = new ArrayList<>();
 162     args.add(command);
 163     args.add("-o");
 164     args.add(getFilePath(searchOutputFile));
 165     args.add("-A");
 166     args.add(getFilePath(hitsAlignmentFile));
 167
 168     boolean dbFound = false;
 169     String dbPath = "";
 170     File databaseFile = null;
 171
 172     boolean useEvalueCutoff = false;
 173     boolean useScoreCutoff = false;
 174     String seqEvalueCutoff = null;
 175     String domEvalueCutoff = null;
 176     String seqScoreCutoff = null;
 177     String domScoreCutoff = null;
 178     databaseName = "Alignment";
 179
 180     if (params != null)
 181     {
 182       for (ArgumentI arg : params)
 183       {
 184         String name = arg.getName();
 185         if (MessageManager.getString(NUMBER_OF_RESULTS_KEY)
 186                 .equals(name))
 187         {
 188           seqsToReturn = Integer.parseInt(arg.getValue());
 189         }
 190         else if (MessageManager.getString(AUTO_ALIGN_SEQS_KEY)
 191                 .equals(name))
 192         {
 193           realign = true;
 194         }
 195         else if (MessageManager.getString(USE_ACCESSIONS_KEY)
 196                 .equals(name))
 197         {
 198           args.add("--acc");
 199         }
 200         else if (MessageManager.getString(REPORTING_CUTOFF_KEY)
 201                 .equals(name))
 202         {
 203           if (CUTOFF_EVALUE.equals(arg.getValue()))
 204           {
 205             useEvalueCutoff = true;
 206           }
 207           else if (CUTOFF_SCORE.equals(arg.getValue()))
 208           {
 209             useScoreCutoff = true;
 210           }
 211         }
 212         else if (MessageManager.getString(SEQ_EVALUE_KEY).equals(name))
 213         {
 214           seqEvalueCutoff = arg.getValue();
 215         }
 216         else if (MessageManager.getString(SEQ_SCORE_KEY).equals(name))
 217         {
 218           seqScoreCutoff = arg.getValue();
 219         }
 220         else if (MessageManager.getString(DOM_EVALUE_KEY)
 221                 .equals(name))
 222         {
 223           domEvalueCutoff = arg.getValue();
 224         }
 225         else if (MessageManager.getString(DOM_SCORE_KEY).equals(name))
 226         {
 227           domScoreCutoff = arg.getValue();
 228         }
 229         else if (MessageManager.getString(TRIM_TERMINI_KEY)
 230                 .equals(name))
 231         {
 232           trim = true;
 233         }
 234         else if (MessageManager.getString(DATABASE_KEY).equals(name))
 235         {
 236           dbFound = true;
 237           dbPath = arg.getValue();
 238           if (!MessageManager.getString(THIS_ALIGNMENT_KEY)
 239                   .equals(dbPath))
 240           {
 241             int pos = dbPath.lastIndexOf(File.separator);
 242             databaseName = dbPath.substring(pos + 1);
 243             databaseFile = new File(dbPath);
 244           }
 245         }
 246       }
 247     }
 248
 249     if (useEvalueCutoff)
 250     {
 251       args.add("-E");
 252       args.add(seqEvalueCutoff);
 253       args.add("--domE");
 254       args.add(domEvalueCutoff);
 255     }
 256     else if (useScoreCutoff)
 257     {
 258       args.add("-T");
 259       args.add(seqScoreCutoff);
 260       args.add("--domT");
 261       args.add(domScoreCutoff);
 262     }
 263
 264     if (!dbFound || MessageManager.getString(THIS_ALIGNMENT_KEY)
 265             .equals(dbPath))
 266     {
 267       /*
 268        * no external database specified for search, so
 269        * export current alignment as 'database' to search,
 270        * excluding any HMM consensus sequences it contains
 271        */
 272       databaseFile = FileUtils.createTempFile("database", ".sto");
 273       AlignmentI al = af.getViewport().getAlignment();
 274       AlignmentI copy = new Alignment(al);
 275       List<SequenceI> hmms = copy.getHmmSequences();
 276       for (SequenceI hmmSeq : hmms)
 277       {
 278         copy.deleteSequence(hmmSeq);
 279       }
 280       exportStockholm(copy.getSequencesArray(), databaseFile, null);
 281     }
 282
 283     args.add(getFilePath(hmmFile));
 284     args.add(getFilePath(databaseFile));
 285
 286     return runCommand(args);
 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       readTable(searchOutputFile);
 314
 315       int seqCount = Math.min(seqs.length, seqsToReturn);
 316       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 317       hmmAndSeqs[0] = hmmSeq;
 318       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 319
 320       if (realign)
 321       {
 322         realignResults(hmmAndSeqs);
 323       }
 324       else
 325       {
 326         AlignmentI al = new Alignment(hmmAndSeqs);
 327         AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 328                 AlignFrame.DEFAULT_HEIGHT);
 329         String ttl = "hmmSearch of " + databaseName + " using "
 330                 + hmmSeq.getName();
 331         Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 332                 AlignFrame.DEFAULT_HEIGHT);
 333       }
 334
 335       hmmTemp.delete();
 336       inputAlignmentTemp.delete();
 337       searchOutputFile.delete();
 338     } finally
 339     {
 340       if (br != null)
 341       {
 342         br.close();
 343       }
 344     }
 345   }
 346
 347   /**
 348    * Realigns the given sequences using hmmalign, to the HMM profile sequence
 349    * which is the first in the array, and opens the results in a new frame
 350    *
 351    * @param hmmAndSeqs
 352    */
 353   protected void realignResults(SequenceI[] hmmAndSeqs)
 354   {
 355     /*
 356      * and align the search results to the HMM profile
 357      */
 358     AlignmentI al = new Alignment(hmmAndSeqs);
 359     AlignFrame frame = new AlignFrame(al, 1, 1);
 360     List<ArgumentI> alignArgs = new ArrayList<>();
 361     String alignTo = hmmAndSeqs[0].getName();
 362     List<String> options = Collections.singletonList(alignTo);
 363     Option option = new Option(MessageManager.getString("label.use_hmm"),
 364             "", true, alignTo, alignTo, options, null);
 365     alignArgs.add(option);
 366     if (trim)
 367     {
 368       alignArgs.add(new BooleanOption(
 369               MessageManager.getString(TRIM_TERMINI_KEY),
 370               MessageManager.getString("label.trim_termini_desc"), true,
 371               true, true, null));
 372     }
 373     HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 374     hmmalign.run();
 375   }
 376
 377   /**
 378    * Reads in the scores table output by hmmsearch and adds annotation to
 379    * sequences for E-value and bit score
 380    *
 381    * @param inputTableTemp
 382    * @throws IOException
 383    */
 384   void readTable(File inputTableTemp) throws IOException
 385   {
 386     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 387     String line = "";
 388     while (!line.startsWith("Query:"))
 389     {
 390       line = br.readLine();
 391     }
 392     for (int i = 0; i < 5; i++)
 393     {
 394       line = br.readLine();
 395     }
 396
 397     int index = 0;
 398     while (!"  ------ inclusion threshold ------".equals(line)
 399             && !"".equals(line))
 400     {
 401       Scanner scanner = new Scanner(line);
 402
 403       String str = scanner.next(); // full sequence eValue score
 404       // float eValue = Float.parseFloat(str);
 405       // int seqLength = seqs[index].getLength();
 406       // Annotation[] annots = new Annotation[seqLength];
 407       // for (int j = 0; j < seqLength; j++)
 408       // {
 409       // annots[j] = new Annotation(eValue);
 410       // }
 411       AlignmentAnnotation annot = new AlignmentAnnotation("E-value",
 412               "Score", null);
 413       annot.setCalcId(HMMSEARCH);
 414       double eValue = Double.parseDouble(str);
 415       annot.setScore(eValue);
 416       annot.setSequenceRef(seqs[index]);
 417       seqs[index].addAlignmentAnnotation(annot);
 418
 419       scanner.close();
 420       line = br.readLine();
 421       index++;
 422     }
 423
 424     br.close();
 425   }
 426
 427 }