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;
21 import java.io.BufferedReader;
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;
30 import javax.swing.JOptionPane;
32 public class HMMSearch extends HmmerCommand
34 private static final String PARAMNAME_NO_OF_RESULTS = MessageManager.getString("label.number_of_results");
36 static final String HMMSEARCH = "hmmsearch";
38 boolean realign = false;
42 int seqsToReturn = Integer.MAX_VALUE;
47 * Constructor for the HMMSearchThread
51 public HMMSearch(AlignFrame af, List<ArgumentI> args)
57 * Runs the HMMSearchThread: the data on the alignment or group is exported,
58 * then the command is executed in the command line and then the data is
59 * imported and displayed in a new frame. Call this method directly to execute
60 * synchronously, or via start() in a new Thread for asynchronously.
65 HiddenMarkovModel hmm = getHmmProfile();
68 // shouldn't happen if we got this far
69 Cache.log.error("Error: no hmm for hmmsearch");
73 SequenceI hmmSeq = hmm.getConsensusSequence();// af.getSelectedHMMSequence();
74 long msgId = System.currentTimeMillis();
75 af.setProgressBar(MessageManager.getString("status.running_hmmsearch"),
80 File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
81 File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
83 File searchOutputFile = FileUtils.createTempFile("searchOutput",
86 exportHmm(hmm, hmmFile.getAbsoluteFile());
88 boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
91 JvOptionPane.showInternalMessageDialog(af, MessageManager
92 .formatMessage("warn.command_failed", "hmmsearch"));
96 importData(hmmSeq, hitsAlignmentFile, hmmFile, searchOutputFile);
97 // TODO make realignment of search results a step at this level
98 // and make it conditional on this.realign
99 } catch (IOException | InterruptedException e)
105 af.setProgressBar("", msgId);
110 * Executes an hmmsearch with the given hmm as input. The database to be
111 * searched is a local file as specified by the 'Database' parameter, or the
112 * current alignment (written to file) if none is specified.
114 * @param searchOutputFile
115 * @param hitsAlignmentFile
119 * @throws IOException
121 private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
122 File hmmFile) throws IOException
124 String command = getCommandPath(HMMSEARCH);
130 List<String> args = new ArrayList<>();
133 args.add(getFilePath(searchOutputFile));
135 args.add(getFilePath(hitsAlignmentFile));
137 boolean dbFound = false;
139 File databaseFile = null;
143 for (ArgumentI arg : params)
145 String name = arg.getName();
146 if (MessageManager.getString("label.number_of_results")
149 seqsToReturn = Integer.parseInt(arg.getValue());
151 else if (MessageManager.getString("label.auto_align_seqs")
154 realign = true; // TODO: not used
156 else if (MessageManager.getString("label.use_accessions")
161 else if (MessageManager.getString("label.seq_e_value").equals(name))
164 args.add(arg.getValue());
166 else if (MessageManager.getString("label.seq_score").equals(name))
169 args.add(arg.getValue());
171 else if (MessageManager.getString("label.dom_e_value_desc")
174 args.add("--incdomE");
175 args.add(arg.getValue());
177 else if (MessageManager.getString("label.dom_score").equals(name))
179 args.add("--incdomT");
180 args.add(arg.getValue());
182 else if (MessageManager.getString("label.trim_termini")
187 else if (MessageManager.getString("label.database").equals(name))
190 dbPath = arg.getValue();
191 if (!MessageManager.getString("label.this_alignment")
194 databaseFile = new File(dbPath);
200 if (!dbFound || MessageManager.getString("label.this_alignment")
204 * no external database specified for search, so
205 * export current alignment as 'database' to search,
206 * excluding any HMM consensus sequences it contains
208 databaseFile = FileUtils.createTempFile("database", ".sto");
209 AlignmentI al = af.getViewport().getAlignment();
210 AlignmentI copy = new Alignment(al);
211 List<SequenceI> hmms = copy.getHmmSequences();
212 for (SequenceI hmmSeq : hmms)
214 copy.deleteSequence(hmmSeq);
216 exportStockholm(copy.getSequencesArray(), databaseFile, null);
217 // StockholmFile stoFile = new StockholmFile(copy);
218 // stoFile.setSeqs(copy.getSequencesArray());
219 // String alignmentString = stoFile.print();
220 // PrintWriter writer = new PrintWriter(databaseFile);
221 // writer.print(alignmentString);
225 args.add(getFilePath(hmmFile));
226 args.add(getFilePath(databaseFile));
228 return runCommand(args);
232 * Imports the data from the temporary file to which the output of hmmsearch
237 private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
238 File hmmTemp, File searchOutputFile)
239 throws IOException, InterruptedException
241 BufferedReader br = new BufferedReader(
242 new FileReader(inputAlignmentTemp));
245 if (br.readLine() == null)
247 JOptionPane.showMessageDialog(af,
248 MessageManager.getString("label.no_sequences_found"));
251 StockholmFile file = new StockholmFile(new FileParse(
252 inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
253 seqs = file.getSeqsAsArray();
255 readTable(searchOutputFile);
257 int seqCount = Math.min(seqs.length, seqsToReturn);
258 SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
259 hmmAndSeqs[0] = hmmSeq;
260 System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
263 * and align the search results to the HMM profile
265 AlignmentI al = new Alignment(hmmAndSeqs);
266 AlignFrame frame = new AlignFrame(al, 1, 1);
267 List<ArgumentI> alignArgs = new ArrayList<>();
268 String defSeq = hmmSeq.getName();
269 List<String> options = Collections.singletonList(defSeq);
270 Option option = new Option(MessageManager.getString("label.use_hmm"),
271 "", true, defSeq, defSeq, options, null);
272 alignArgs.add(option);
275 alignArgs.add(new BooleanOption(
276 MessageManager.getString("label.trim_termini"),
277 MessageManager.getString("label.trim_termini_desc"), true,
280 HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
284 inputAlignmentTemp.delete();
285 searchOutputFile.delete();
295 void readTable(File inputTableTemp) throws IOException
297 BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
299 while (!line.startsWith("Query:"))
301 line = br.readLine();
303 for (int i = 0; i < 5; i++)
305 line = br.readLine();
309 while (!" ------ inclusion threshold ------".equals(line)
312 Scanner scanner = new Scanner(line);
314 String str = scanner.next(); // full sequence eValue score
315 float eValue = Float.parseFloat(str);
316 int seqLength = seqs[index].getLength();
317 Annotation[] annots = new Annotation[seqLength];
318 for (int j = 0; j < seqLength; j++)
320 annots[j] = new Annotation(eValue);
322 AlignmentAnnotation annot = new AlignmentAnnotation("E-value",
324 annot.setScore(Double.parseDouble(str));
325 annot.setSequenceRef(seqs[index]);
326 seqs[index].addAlignmentAnnotation(annot);
329 line = br.readLine();