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.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   /**
  76    * Constructor for the HMMSearchThread
  77    *
  78    * @param af
  79    */
  80   public HMMSearch(AlignFrame af, List<ArgumentI> args)
  81   {
  82     super(af, args);
  83   }
  84
  85   /**
  86    * Runs the HMMSearchThread: the data on the alignment or group is exported,
  87    * then the command is executed in the command line and then the data is
  88    * imported and displayed in a new frame. Call this method directly to execute
  89    * synchronously, or via start() in a new Thread for asynchronously.
  90    */
  91   @Override
  92   public void run()
  93   {
  94     HiddenMarkovModel hmm = getHmmProfile();
  95     if (hmm == null)
  96     {
  97       // shouldn't happen if we got this far
  98       Cache.log.error("Error: no hmm for hmmsearch");
  99       return;
 100     }
 101
 102     SequenceI hmmSeq = hmm.getConsensusSequence();// af.getSelectedHMMSequence();
 103     long msgId = System.currentTimeMillis();
 104     af.setProgressBar(MessageManager.getString("status.running_hmmsearch"),
 105             msgId);
 106
 107     try
 108     {
 109       File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
 110       File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
 111               ".sto");
 112       File searchOutputFile = FileUtils.createTempFile("searchOutput",
 113               ".sto");
 114
 115       exportHmm(hmm, hmmFile.getAbsoluteFile());
 116
 117       boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
 118       if (!ran)
 119       {
 120         JvOptionPane.showInternalMessageDialog(af, MessageManager
 121                 .formatMessage("warn.command_failed", "hmmsearch"));
 122         return;
 123       }
 124
 125       importData(hmmSeq, hitsAlignmentFile, hmmFile, searchOutputFile);
 126       // TODO make realignment of search results a step at this level
 127       // and make it conditional on this.realign
 128     } catch (IOException | InterruptedException e)
 129     {
 130       e.printStackTrace();
 131     }
 132     finally
 133     {
 134       af.setProgressBar("", msgId);
 135     }
 136   }
 137
 138   /**
 139    * Executes an hmmsearch with the given hmm as input. The database to be
 140    * searched is a local file as specified by the 'Database' parameter, or the
 141    * current alignment (written to file) if none is specified.
 142    *
 143    * @param searchOutputFile
 144    * @param hitsAlignmentFile
 145    * @param hmmFile
 146    *
 147    * @return
 148    * @throws IOException
 149    */
 150   private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
 151           File hmmFile) throws IOException
 152   {
 153     String command = getCommandPath(HMMSEARCH);
 154     if (command == null)
 155     {
 156       return false;
 157     }
 158
 159     List<String> args = new ArrayList<>();
 160     args.add(command);
 161     args.add("-o");
 162     args.add(getFilePath(searchOutputFile));
 163     args.add("-A");
 164     args.add(getFilePath(hitsAlignmentFile));
 165
 166     boolean dbFound = false;
 167     String dbPath = "";
 168     File databaseFile = null;
 169
 170     boolean useEvalueCutoff = false;
 171     boolean useScoreCutoff = false;
 172     String seqEvalueCutoff = null;
 173     String domEvalueCutoff = null;
 174     String seqScoreCutoff = null;
 175     String domScoreCutoff = null;
 176
 177     if (params != null)
 178     {
 179       for (ArgumentI arg : params)
 180       {
 181         String name = arg.getName();
 182         if (MessageManager.getString(NUMBER_OF_RESULTS_KEY)
 183                 .equals(name))
 184         {
 185           seqsToReturn = Integer.parseInt(arg.getValue());
 186         }
 187         else if (MessageManager.getString(AUTO_ALIGN_SEQS_KEY)
 188                 .equals(name))
 189         {
 190           realign = true; // TODO: not used
 191         }
 192         else if (MessageManager.getString(USE_ACCESSIONS_KEY)
 193                 .equals(name))
 194         {
 195           args.add("--acc");
 196         }
 197         else if (MessageManager.getString(REPORTING_CUTOFF_KEY)
 198                 .equals(name))
 199         {
 200           if (CUTOFF_EVALUE.equals(arg.getValue()))
 201           {
 202             useEvalueCutoff = true;
 203           }
 204           else if (CUTOFF_SCORE.equals(arg.getValue()))
 205           {
 206             useScoreCutoff = true;
 207           }
 208         }
 209         else if (MessageManager.getString(SEQ_EVALUE_KEY).equals(name))
 210         {
 211           seqEvalueCutoff = arg.getValue();
 212         }
 213         else if (MessageManager.getString(SEQ_SCORE_KEY).equals(name))
 214         {
 215           seqScoreCutoff = arg.getValue();
 216         }
 217         else if (MessageManager.getString(DOM_EVALUE_KEY)
 218                 .equals(name))
 219         {
 220           domEvalueCutoff = arg.getValue();
 221         }
 222         else if (MessageManager.getString(DOM_SCORE_KEY).equals(name))
 223         {
 224           domScoreCutoff = arg.getValue();
 225         }
 226         else if (MessageManager.getString(TRIM_TERMINI_KEY)
 227                 .equals(name))
 228         {
 229           trim = true;
 230         }
 231         else if (MessageManager.getString(DATABASE_KEY).equals(name))
 232         {
 233           dbFound = true;
 234           dbPath = arg.getValue();
 235           if (!MessageManager.getString(THIS_ALIGNMENT_KEY)
 236                   .equals(dbPath))
 237           {
 238             databaseFile = new File(dbPath);
 239           }
 240         }
 241       }
 242     }
 243
 244     if (useEvalueCutoff)
 245     {
 246       args.add("-E");
 247       args.add(seqEvalueCutoff);
 248       args.add("--domE");
 249       args.add(domEvalueCutoff);
 250     }
 251     else if (useScoreCutoff)
 252     {
 253       args.add("-T");
 254       args.add(seqScoreCutoff);
 255       args.add("--domT");
 256       args.add(domScoreCutoff);
 257     }
 258
 259     if (!dbFound || MessageManager.getString(THIS_ALIGNMENT_KEY)
 260             .equals(dbPath))
 261     {
 262       /*
 263        * no external database specified for search, so
 264        * export current alignment as 'database' to search,
 265        * excluding any HMM consensus sequences it contains
 266        */
 267       databaseFile = FileUtils.createTempFile("database", ".sto");
 268       AlignmentI al = af.getViewport().getAlignment();
 269       AlignmentI copy = new Alignment(al);
 270       List<SequenceI> hmms = copy.getHmmSequences();
 271       for (SequenceI hmmSeq : hmms)
 272       {
 273         copy.deleteSequence(hmmSeq);
 274       }
 275       exportStockholm(copy.getSequencesArray(), databaseFile, null);
 276     }
 277
 278     args.add(getFilePath(hmmFile));
 279     args.add(getFilePath(databaseFile));
 280
 281     return runCommand(args);
 282   }
 283
 284   /**
 285    * Imports the data from the temporary file to which the output of hmmsearch
 286    * is directed.
 287    *
 288    * @param hmmSeq
 289    */
 290   private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
 291           File hmmTemp, File searchOutputFile)
 292           throws IOException, InterruptedException
 293   {
 294     BufferedReader br = new BufferedReader(
 295             new FileReader(inputAlignmentTemp));
 296     try
 297     {
 298       if (br.readLine() == null)
 299       {
 300         JOptionPane.showMessageDialog(af,
 301                 MessageManager.getString("label.no_sequences_found"));
 302         return;
 303       }
 304       StockholmFile file = new StockholmFile(new FileParse(
 305               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
 306       seqs = file.getSeqsAsArray();
 307
 308       readTable(searchOutputFile);
 309
 310       int seqCount = Math.min(seqs.length, seqsToReturn);
 311       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 312       hmmAndSeqs[0] = hmmSeq;
 313       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 314
 315       /*
 316        * and align the search results to the HMM profile
 317        */
 318       AlignmentI al = new Alignment(hmmAndSeqs);
 319       AlignFrame frame = new AlignFrame(al, 1, 1);
 320       List<ArgumentI> alignArgs = new ArrayList<>();
 321       String defSeq = hmmSeq.getName();
 322       List<String> options = Collections.singletonList(defSeq);
 323       Option option = new Option(MessageManager.getString("label.use_hmm"),
 324               "", true, defSeq, defSeq, options, null);
 325       alignArgs.add(option);
 326       if (trim)
 327       {
 328         alignArgs.add(new BooleanOption(
 329                 MessageManager.getString(TRIM_TERMINI_KEY),
 330                 MessageManager.getString("label.trim_termini_desc"), true,
 331                 true, true, null));
 332       }
 333       HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 334       hmmalign.run();
 335       frame = null;
 336       hmmTemp.delete();
 337       inputAlignmentTemp.delete();
 338       searchOutputFile.delete();
 339     } finally
 340     {
 341       if (br != null)
 342       {
 343         br.close();
 344       }
 345     }
 346   }
 347
 348   void readTable(File inputTableTemp) throws IOException
 349   {
 350     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 351     String line = "";
 352     while (!line.startsWith("Query:"))
 353     {
 354       line = br.readLine();
 355     }
 356     for (int i = 0; i < 5; i++)
 357     {
 358       line = br.readLine();
 359     }
 360
 361     int index = 0;
 362     while (!"  ------ inclusion threshold ------".equals(line)
 363             && !"".equals(line))
 364     {
 365       Scanner scanner = new Scanner(line);
 366
 367       String str = scanner.next(); // full sequence eValue score
 368       float eValue = Float.parseFloat(str);
 369       int seqLength = seqs[index].getLength();
 370       Annotation[] annots = new Annotation[seqLength];
 371       for (int j = 0; j < seqLength; j++)
 372       {
 373         annots[j] = new Annotation(eValue);
 374       }
 375       AlignmentAnnotation annot = new AlignmentAnnotation("E-value",
 376               "Score", annots);
 377       annot.setScore(Double.parseDouble(str));
 378       annot.setSequenceRef(seqs[index]);
 379       seqs[index].addAlignmentAnnotation(annot);
 380
 381       scanner.close();
 382       line = br.readLine();
 383       index++;
 384     }
 385
 386     br.close();
 387   }
 388
 389 }