JAL-2629 fix hmsearch and jackhmmer not being able to search alignment
[jalview.git] / src / jalview / hmmer / JackHMMER.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.SequenceI;
8 import jalview.gui.AlignFrame;
9 import jalview.gui.Desktop;
10 import jalview.gui.JvOptionPane;
11 import jalview.io.DataSourceType;
12 import jalview.io.FileParse;
13 import jalview.io.StockholmFile;
14 import jalview.util.FileUtils;
15 import jalview.util.MessageManager;
16 import jalview.ws.params.ArgumentI;
17
18 import java.io.BufferedReader;
19 import java.io.File;
20 import java.io.FileReader;
21 import java.io.IOException;
22 import java.util.ArrayList;
23 import java.util.Hashtable;
24 import java.util.List;
25 import java.util.Scanner;
26
27 import javax.swing.JOptionPane;
28
29 public class JackHMMER extends HmmerCommand
30 {
31   static final String JACKHMMER = "jackhmmer";
32
33   boolean realign = false;
34
35   boolean trim = false;
36
37   int seqsToReturn = Integer.MAX_VALUE;
38
39   SequenceI[] seqs;
40
41   private String databaseName;
42
43   boolean searchAlignment = true;
44
45   Hashtable sequencesHash;
46
47   /**
48    * Constructor for the JackhmmerThread
49    * 
50    * @param af
51    */
52   public JackHMMER(AlignFrame af, List<ArgumentI> args)
53   {
54     super(af, args);
55   }
56
57   /**
58    * Runs the JackhmmerThread: the data on the alignment or group is exported,
59    * then the command is executed in the command line and then the data is
60    * imported and displayed in a new frame. Call this method directly to execute
61    * synchronously, or via start() in a new Thread for asynchronously.
62    */
63   @Override
64   public void run()
65   {
66     SequenceI seq = getSequence();
67     if (seq == null)
68     {
69       // shouldn't happen if we got this far
70       Cache.log.error("Error: no sequence for jackhmmer");
71       return;
72     }
73
74     long msgId = System.currentTimeMillis();
75     af.setProgressBar(MessageManager.getString("status.running_search"),
76             msgId);
77
78     try
79     {
80       File seqFile = FileUtils.createTempFile("seq", ".sto");
81       File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
82               ".sto");
83       File searchOutputFile = FileUtils.createTempFile("searchOutput",
84               ".txt");
85
86       exportStockholm(new SequenceI[] { seq }, seqFile.getAbsoluteFile(),
87               null);
88
89       boolean ran = runCommand(searchOutputFile, hitsAlignmentFile,
90               seqFile);
91       if (!ran)
92       {
93         JvOptionPane.showInternalMessageDialog(af, MessageManager
94                 .formatMessage("warn.command_failed", "jackhmmer"));
95         return;
96       }
97
98       importData(hitsAlignmentFile, seqFile, searchOutputFile);
99       // TODO make realignment of search results a step at this level
100       // and make it conditional on this.realign
101     } catch (IOException | InterruptedException e)
102     {
103       e.printStackTrace();
104     } finally
105     {
106       af.setProgressBar("", msgId);
107     }
108   }
109
110   /**
111    * Executes an jackhmmer search with the given sequence as input. The database
112    * to be searched is a local file as specified by the 'Database' parameter, or
113    * the current alignment (written to file) if none is specified.
114    * 
115    * @param searchOutputFile
116    * @param hitsAlignmentFile
117    * @param seqFile
118    * 
119    * @return
120    * @throws IOException
121    */
122   private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
123           File seqFile) throws IOException
124   {
125     String command = getCommandPath(JACKHMMER);
126     if (command == null)
127     {
128       return false;
129     }
130
131     List<String> args = new ArrayList<>();
132     args.add(command);
133     buildArguments(args, searchOutputFile, hitsAlignmentFile, seqFile);
134
135     return runCommand(args);
136   }
137
138   /**
139    * Appends command line arguments to the given list, to specify input and output
140    * files for the search, and any additional options that may have been passed
141    * from the parameters dialog
142    * 
143    * @param args
144    * @param searchOutputFile
145    * @param hitsAlignmentFile
146    * @param seqFile
147    * @throws IOException
148    */
149   protected void buildArguments(List<String> args, File searchOutputFile,
150           File hitsAlignmentFile, File seqFile) throws IOException
151   {
152     args.add("-o");
153     args.add(getFilePath(searchOutputFile, true));
154     args.add("-A");
155     args.add(getFilePath(hitsAlignmentFile, true));
156
157     File databaseFile = null;
158
159     boolean useEvalueCutoff = false;
160     boolean useScoreCutoff = false;
161     String seqEvalueCutoff = null;
162     String domEvalueCutoff = null;
163     String seqScoreCutoff = null;
164     String domScoreCutoff = null;
165     databaseName = "Alignment";
166
167     if (params != null)
168     {
169       for (ArgumentI arg : params)
170       {
171         String name = arg.getName();
172
173         if (MessageManager.getString(REPORTING_CUTOFF_KEY)
174                 .equals(name))
175         {
176           if (MessageManager.getString(CUTOFF_EVALUE)
177                   .equals(arg.getValue()))
178           {
179             useEvalueCutoff = true;
180           }
181           else if (MessageManager.getString(CUTOFF_SCORE)
182                   .equals(arg.getValue()))
183           {
184             useScoreCutoff = true;
185           }
186         }
187         else if (MessageManager.getString(SEQ_EVALUE_KEY).equals(name))
188         {
189           seqEvalueCutoff = arg.getValue();
190         }
191         else if (MessageManager.getString(SEQ_SCORE_KEY).equals(name))
192         {
193           seqScoreCutoff = arg.getValue();
194         }
195         else if (MessageManager.getString(DOM_EVALUE_KEY).equals(name))
196         {
197           domEvalueCutoff = arg.getValue();
198         }
199         else if (MessageManager.getString(DOM_SCORE_KEY).equals(name))
200         {
201           domScoreCutoff = arg.getValue();
202         }
203         else if (MessageManager.getString(DATABASE_KEY).equals(name))
204         {
205           databaseFile = new File(arg.getValue());
206           if (!arg.getValue().isEmpty())
207           {
208             searchAlignment = false;
209           }
210         }
211       }
212     }
213
214     if (useEvalueCutoff)
215     {
216       args.add("-E");
217       args.add(seqEvalueCutoff);
218       args.add("--domE");
219       args.add(domEvalueCutoff);
220     }
221     else if (useScoreCutoff)
222     {
223       args.add("-T");
224       args.add(seqScoreCutoff);
225       args.add("--domT");
226       args.add(domScoreCutoff);
227     }
228
229     // if (!dbFound || MessageManager.getString(THIS_ALIGNMENT_KEY)
230     // .equals(dbPath))
231     if (searchAlignment)
232     {
233       /*
234        * no external database specified for search, so
235        * export current alignment as 'database' to search
236        */
237       databaseFile = FileUtils.createTempFile("database", ".sto");
238       AlignmentI al = af.getViewport().getAlignment();
239       AlignmentI copy = new Alignment(al);
240
241       deleteHmmSequences(copy);
242
243       if (searchAlignment)
244       {
245         sequencesHash = stashSequences(copy.getSequencesArray());
246       }
247
248       exportStockholm(copy.getSequencesArray(), databaseFile, null);
249     }
250
251     args.add(getFilePath(seqFile, true));
252     args.add(getFilePath(databaseFile, true));
253   }
254
255   /**
256    * Imports the data from the temporary file to which the output of jackhmmer was
257    * directed.
258    */
259   private void importData(File inputAlignmentTemp, File seqTemp,
260           File searchOutputFile) throws IOException, InterruptedException
261   {
262     BufferedReader br = new BufferedReader(
263             new FileReader(inputAlignmentTemp));
264     try
265     {
266       if (br.readLine() == null)
267       {
268         JOptionPane.showMessageDialog(af,
269                 MessageManager.getString("label.no_sequences_found"));
270         return;
271       }
272       StockholmFile file = new StockholmFile(new FileParse(
273               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
274       seqs = file.getSeqsAsArray();
275
276       if (searchAlignment)
277       {
278         recoverSequences(sequencesHash, seqs);
279       }
280
281       readTable(searchOutputFile);
282
283       int seqCount = Math.min(seqs.length, seqsToReturn);
284
285       AlignmentI al = new Alignment(seqs);
286
287       AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
288               AlignFrame.DEFAULT_HEIGHT);
289       String ttl = "jackhmmer search of " + databaseName + " using "
290               + seqs[0].getName();
291       Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
292               AlignFrame.DEFAULT_HEIGHT);
293
294       seqTemp.delete();
295       inputAlignmentTemp.delete();
296       searchOutputFile.delete();
297     } finally
298     {
299       if (br != null)
300       {
301         br.close();
302       }
303     }
304   }
305
306   /**
307    * Reads in the scores table output by jackhmmer and adds annotation to
308    * sequences for E-value and bit score
309    * 
310    * @param inputTableTemp
311    * @throws IOException
312    */
313   void readTable(File inputTableTemp) throws IOException
314   {
315     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
316     String line = "";
317     while (!line.startsWith("Query:"))
318     {
319       line = br.readLine();
320     }
321     while (!line.contains("-------"))
322     {
323       line = br.readLine();
324     }
325     line = br.readLine();
326
327     int index = 0;
328     while (!"  ------ inclusion threshold ------".equals(line)
329             && !"".equals(line))
330     {
331       SequenceI seq = seqs[index];
332
333       Scanner scanner = new Scanner(line);
334       String evalue = scanner.next();
335       evalue = scanner.next();
336       String score = scanner.next();
337       addScoreAnnotations(evalue, score, seq);
338       scanner.close();
339       line = br.readLine();
340       index++;
341     }
342
343     br.close();
344   }
345
346   protected void addScoreAnnotations(String eValue, String bitScore,
347           SequenceI seq)
348   {
349     String label = "Search Scores";
350     String description = "Full sequence bit score and E-Value";
351
352     try
353     {
354       AlignmentAnnotation annot = new AlignmentAnnotation(label,
355               description, null);
356
357       annot.label = label;
358       annot.description = description;
359
360       annot.setCalcId(JACKHMMER);
361
362       double dEValue = Double.parseDouble(eValue);
363       annot.setEValue(dEValue);
364
365       double dBitScore = Double.parseDouble(bitScore);
366       annot.setBitScore(dBitScore);
367
368       annot.setSequenceRef(seq);
369       seq.addAlignmentAnnotation(annot);
370
371     } catch (NumberFormatException e)
372     {
373       System.err.println("Error parsing " + label + " from " + eValue
374               + " & " + bitScore);
375     }
376   }
377
378
379
380 }