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;
22 import java.io.BufferedReader;
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;
31 import javax.swing.JOptionPane;
33 public class HMMSearch extends HmmerCommand
35 static final String HMMSEARCH = "hmmsearch";
38 * constants for i18n lookup of passed parameter names
40 static final String DATABASE_KEY = "label.database";
42 static final String THIS_ALIGNMENT_KEY = "label.this_alignment";
44 static final String USE_ACCESSIONS_KEY = "label.use_accessions";
46 static final String AUTO_ALIGN_SEQS_KEY = "label.auto_align_seqs";
48 static final String NUMBER_OF_RESULTS_KEY = "label.number_of_results";
50 static final String TRIM_TERMINI_KEY = "label.trim_termini";
52 static final String REPORTING_CUTOFF_KEY = "label.reporting_cutoff";
54 static final String CUTOFF_NONE = "None";
56 static final String CUTOFF_SCORE = "Score";
58 static final String CUTOFF_EVALUE = "E-Value";
60 static final String SEQ_EVALUE_KEY = "label.seq_evalue";
62 static final String DOM_EVALUE_KEY = "label.dom_evalue";
64 static final String SEQ_SCORE_KEY = "label.seq_score";
66 static final String DOM_SCORE_KEY = "label.dom_score";
68 boolean realign = false;
72 int seqsToReturn = Integer.MAX_VALUE;
76 private String databaseName;
79 * Constructor for the HMMSearchThread
83 public HMMSearch(AlignFrame af, List<ArgumentI> args)
89 * Runs the HMMSearchThread: the data on the alignment or group is exported,
90 * then the command is executed in the command line and then the data is
91 * imported and displayed in a new frame. Call this method directly to execute
92 * synchronously, or via start() in a new Thread for asynchronously.
97 HiddenMarkovModel hmm = getHmmProfile();
100 // shouldn't happen if we got this far
101 Cache.log.error("Error: no hmm for hmmsearch");
105 SequenceI hmmSeq = hmm.getConsensusSequence();
106 long msgId = System.currentTimeMillis();
107 af.setProgressBar(MessageManager.getString("status.running_hmmsearch"),
112 File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
113 File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
115 File searchOutputFile = FileUtils.createTempFile("searchOutput",
118 exportHmm(hmm, hmmFile.getAbsoluteFile());
120 boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
123 JvOptionPane.showInternalMessageDialog(af, MessageManager
124 .formatMessage("warn.command_failed", "hmmsearch"));
128 importData(hmmSeq, hitsAlignmentFile, hmmFile, searchOutputFile);
129 // TODO make realignment of search results a step at this level
130 // and make it conditional on this.realign
131 } catch (IOException | InterruptedException e)
137 af.setProgressBar("", msgId);
142 * Executes an hmmsearch with the given hmm as input. The database to be
143 * searched is a local file as specified by the 'Database' parameter, or the
144 * current alignment (written to file) if none is specified.
146 * @param searchOutputFile
147 * @param hitsAlignmentFile
151 * @throws IOException
153 private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
154 File hmmFile) throws IOException
156 String command = getCommandPath(HMMSEARCH);
162 List<String> args = new ArrayList<>();
164 buildArguments(args, searchOutputFile, hitsAlignmentFile, hmmFile);
166 return runCommand(args);
170 * Appends command line arguments to the given list, to specify input and
171 * output files for the search, and any additional options that may have been
172 * passed from the parameters dialog
175 * @param searchOutputFile
176 * @param hitsAlignmentFile
178 * @throws IOException
180 protected void buildArguments(List<String> args, File searchOutputFile,
181 File hitsAlignmentFile, File hmmFile) throws IOException
184 args.add(getFilePath(searchOutputFile, true));
186 args.add(getFilePath(hitsAlignmentFile, true));
188 boolean dbFound = false;
190 File databaseFile = null;
192 boolean useEvalueCutoff = false;
193 boolean useScoreCutoff = false;
194 String seqEvalueCutoff = null;
195 String domEvalueCutoff = null;
196 String seqScoreCutoff = null;
197 String domScoreCutoff = null;
198 databaseName = "Alignment";
199 boolean searchAlignment = false;
203 for (ArgumentI arg : params)
205 String name = arg.getName();
206 if (MessageManager.getString(NUMBER_OF_RESULTS_KEY)
209 seqsToReturn = Integer.parseInt(arg.getValue());
211 else if (MessageManager.getString("action.search").equals(name))
213 searchAlignment = arg.getValue().equals(
214 MessageManager.getString(HMMSearch.THIS_ALIGNMENT_KEY));
216 else if (MessageManager.getString(DATABASE_KEY).equals(name))
218 dbPath = arg.getValue();
219 int pos = dbPath.lastIndexOf(File.separator);
220 databaseName = dbPath.substring(pos + 1);
221 databaseFile = new File(dbPath);
223 else if (MessageManager.getString(AUTO_ALIGN_SEQS_KEY)
228 else if (MessageManager.getString(USE_ACCESSIONS_KEY)
233 else if (MessageManager.getString(REPORTING_CUTOFF_KEY)
236 if (CUTOFF_EVALUE.equals(arg.getValue()))
238 useEvalueCutoff = true;
240 else if (CUTOFF_SCORE.equals(arg.getValue()))
242 useScoreCutoff = true;
245 else if (MessageManager.getString(SEQ_EVALUE_KEY).equals(name))
247 seqEvalueCutoff = arg.getValue();
249 else if (MessageManager.getString(SEQ_SCORE_KEY).equals(name))
251 seqScoreCutoff = arg.getValue();
253 else if (MessageManager.getString(DOM_EVALUE_KEY)
256 domEvalueCutoff = arg.getValue();
258 else if (MessageManager.getString(DOM_SCORE_KEY).equals(name))
260 domScoreCutoff = arg.getValue();
262 else if (MessageManager.getString(TRIM_TERMINI_KEY)
267 else if (MessageManager.getString(DATABASE_KEY).equals(name))
270 dbPath = arg.getValue();
271 if (!MessageManager.getString(THIS_ALIGNMENT_KEY)
274 int pos = dbPath.lastIndexOf(File.separator);
275 databaseName = dbPath.substring(pos + 1);
276 databaseFile = new File(dbPath);
285 args.add(seqEvalueCutoff);
287 args.add(domEvalueCutoff);
289 else if (useScoreCutoff)
292 args.add(seqScoreCutoff);
294 args.add(domScoreCutoff);
297 // if (!dbFound || MessageManager.getString(THIS_ALIGNMENT_KEY)
302 * no external database specified for search, so
303 * export current alignment as 'database' to search,
304 * excluding any HMM consensus sequences it contains
306 databaseFile = FileUtils.createTempFile("database", ".sto");
307 AlignmentI al = af.getViewport().getAlignment();
308 AlignmentI copy = new Alignment(al);
309 List<SequenceI> hmms = copy.getHmmSequences();
310 for (SequenceI hmmSeq : hmms)
312 copy.deleteSequence(hmmSeq);
314 exportStockholm(copy.getSequencesArray(), databaseFile, null);
317 args.add(getFilePath(hmmFile, true));
318 args.add(getFilePath(databaseFile, true));
322 * Imports the data from the temporary file to which the output of hmmsearch
323 * was directed. The results are optionally realigned using hmmalign.
327 private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
328 File hmmTemp, File searchOutputFile)
329 throws IOException, InterruptedException
331 BufferedReader br = new BufferedReader(
332 new FileReader(inputAlignmentTemp));
335 if (br.readLine() == null)
337 JOptionPane.showMessageDialog(af,
338 MessageManager.getString("label.no_sequences_found"));
341 StockholmFile file = new StockholmFile(new FileParse(
342 inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
343 seqs = file.getSeqsAsArray();
344 // look for PP cons and ref seq in alignment only annotation
345 AlignmentAnnotation modelpos = null, ppcons = null;
346 for (AlignmentAnnotation aa : file.getAnnotations())
348 if (aa.sequenceRef == null)
350 if (aa.label.equals("Reference Positions")) // RF feature type in
355 if (aa.label.equals("Posterior Probability"))
361 readTable(searchOutputFile);
363 int seqCount = Math.min(seqs.length, seqsToReturn);
364 SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
365 hmmSeq = hmmSeq.deriveSequence(); // otherwise all bad things happen
366 hmmAndSeqs[0] = hmmSeq;
367 System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
368 if (modelpos != null)
370 // TODO need - get ungapped sequence method
372 hmmSeq.getDatasetSequence().getSequenceAsString());
373 Annotation[] refpos = modelpos.annotations;
374 // insert gaps to match with refseq positions
375 int gc = 0, lcol = 0;
376 for (int c = 0; c < refpos.length; c++)
378 if (refpos[c] != null && ("x".equals(refpos[c].displayCharacter)))
382 hmmSeq.insertCharAt(lcol + 1, gc, '-');
395 realignResults(hmmAndSeqs);
399 AlignmentI al = new Alignment(hmmAndSeqs);
402 al.addAnnotation(ppcons);
404 if (modelpos != null)
406 al.addAnnotation(modelpos);
408 AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
409 AlignFrame.DEFAULT_HEIGHT);
410 String ttl = "hmmSearch of " + databaseName + " using "
412 Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
413 AlignFrame.DEFAULT_HEIGHT);
417 inputAlignmentTemp.delete();
418 searchOutputFile.delete();
429 * Realigns the given sequences using hmmalign, to the HMM profile sequence
430 * which is the first in the array, and opens the results in a new frame
434 protected void realignResults(SequenceI[] hmmAndSeqs)
437 * and align the search results to the HMM profile
439 AlignmentI al = new Alignment(hmmAndSeqs);
440 AlignFrame frame = new AlignFrame(al, 1, 1);
441 List<ArgumentI> alignArgs = new ArrayList<>();
442 String alignTo = hmmAndSeqs[0].getName();
443 List<String> options = Collections.singletonList(alignTo);
444 Option option = new Option(MessageManager.getString("label.use_hmm"),
445 "", true, alignTo, alignTo, options, null);
446 alignArgs.add(option);
449 alignArgs.add(new BooleanOption(
450 MessageManager.getString(TRIM_TERMINI_KEY),
451 MessageManager.getString("label.trim_termini_desc"), true,
454 HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
459 * Reads in the scores table output by hmmsearch and adds annotation to
460 * sequences for E-value and bit score
462 * @param inputTableTemp
463 * @throws IOException
465 void readTable(File inputTableTemp) throws IOException
467 BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
469 while (!line.startsWith("Query:"))
471 line = br.readLine();
473 for (int i = 0; i < 5; i++)
475 line = br.readLine();
479 while (!" ------ inclusion threshold ------".equals(line)
482 SequenceI seq = seqs[index];
483 AlignmentAnnotation pp = seq
484 .getAlignmentAnnotations("", "Posterior Probability")
486 Scanner scanner = new Scanner(line);
487 String str = scanner.next();
488 addScoreAnnotation(str, seq, "hmmsearch E-value",
489 "Full sequence E-value", pp);
490 str = scanner.next();
491 addScoreAnnotation(str, seq, "hmmsearch Score",
492 "Full sequence bit score", pp);
493 seq.removeAlignmentAnnotation(pp);
495 line = br.readLine();
503 * A helper method that adds one score-only (non-positional) annotation to a
511 protected void addScoreAnnotation(String value, SequenceI seq,
512 String label, String description)
514 addScoreAnnotation(value, seq, label, description, null);
518 * A helper method that adds one score-only (non-positional) annotation to a
526 * existing posterior probability annotation - values copied to new
529 protected void addScoreAnnotation(String value, SequenceI seq,
530 String label, String description, AlignmentAnnotation pp)
534 AlignmentAnnotation annot = null;
537 new AlignmentAnnotation(label,
542 annot = new AlignmentAnnotation(pp);
544 annot.description = description;
546 annot.setCalcId(HMMSEARCH);
547 double eValue = Double.parseDouble(value);
548 annot.setScore(eValue);
549 annot.setSequenceRef(seq);
550 seq.addAlignmentAnnotation(annot);
551 } catch (NumberFormatException e)
553 System.err.println("Error parsing " + label + " from " + value);