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     boolean dbFound = false;
 158     String dbPath = "";
 159     File databaseFile = null;
 160
 161     boolean useEvalueCutoff = false;
 162     boolean useScoreCutoff = false;
 163     String seqEvalueCutoff = null;
 164     String domEvalueCutoff = null;
 165     String seqScoreCutoff = null;
 166     String domScoreCutoff = null;
 167     databaseName = "Alignment";
 168
 169     if (params != null)
 170     {
 171       for (ArgumentI arg : params)
 172       {
 173         String name = arg.getName();
 174
 175         if (MessageManager.getString(REPORTING_CUTOFF_KEY)
 176                 .equals(name))
 177         {
 178           if (MessageManager.getString(CUTOFF_EVALUE)
 179                   .equals(arg.getValue()))
 180           {
 181             useEvalueCutoff = true;
 182           }
 183           else if (MessageManager.getString(CUTOFF_SCORE)
 184                   .equals(arg.getValue()))
 185           {
 186             useScoreCutoff = true;
 187           }
 188         }
 189         else if (MessageManager.getString(SEQ_EVALUE_KEY).equals(name))
 190         {
 191           seqEvalueCutoff = arg.getValue();
 192         }
 193         else if (MessageManager.getString(SEQ_SCORE_KEY).equals(name))
 194         {
 195           seqScoreCutoff = arg.getValue();
 196         }
 197         else if (MessageManager.getString(DOM_EVALUE_KEY).equals(name))
 198         {
 199           domEvalueCutoff = arg.getValue();
 200         }
 201         else if (MessageManager.getString(DOM_SCORE_KEY).equals(name))
 202         {
 203           domScoreCutoff = arg.getValue();
 204         }
 205         else if (MessageManager.getString(DATABASE_KEY).equals(name))
 206         {
 207           dbFound = true;
 208           dbPath = arg.getValue();
 209           if (!MessageManager.getString(THIS_ALIGNMENT_KEY).equals(dbPath))
 210           {
 211             int pos = dbPath.lastIndexOf(File.separator);
 212             databaseName = dbPath.substring(pos + 1);
 213             databaseFile = new File(dbPath);
 214           }
 215           searchAlignment = false;
 216         }
 217       }
 218     }
 219
 220     if (useEvalueCutoff)
 221     {
 222       args.add("-E");
 223       args.add(seqEvalueCutoff);
 224       args.add("--domE");
 225       args.add(domEvalueCutoff);
 226     }
 227     else if (useScoreCutoff)
 228     {
 229       args.add("-T");
 230       args.add(seqScoreCutoff);
 231       args.add("--domT");
 232       args.add(domScoreCutoff);
 233     }
 234
 235     // if (!dbFound || MessageManager.getString(THIS_ALIGNMENT_KEY)
 236     // .equals(dbPath))
 237     if (searchAlignment)
 238     {
 239       /*
 240        * no external database specified for search, so
 241        * export current alignment as 'database' to search
 242        */
 243       databaseFile = FileUtils.createTempFile("database", ".sto");
 244       AlignmentI al = af.getViewport().getAlignment();
 245       AlignmentI copy = new Alignment(al);
 246
 247       deleteHmmSequences(copy);
 248
 249       sequencesHash = stashSequences(copy.getSequencesArray());
 250
 251       exportStockholm(copy.getSequencesArray(), databaseFile, null);
 252     }
 253
 254     args.add(getFilePath(seqFile, true));
 255     args.add(getFilePath(databaseFile, true));
 256   }
 257
 258   /**
 259    * Imports the data from the temporary file to which the output of jackhmmer was
 260    * directed.
 261    */
 262   private void importData(File inputAlignmentTemp, File seqTemp,
 263           File searchOutputFile) throws IOException, InterruptedException
 264   {
 265     BufferedReader br = new BufferedReader(
 266             new FileReader(inputAlignmentTemp));
 267     try
 268     {
 269       if (br.readLine() == null)
 270       {
 271         JOptionPane.showMessageDialog(af,
 272                 MessageManager.getString("label.no_sequences_found"));
 273         return;
 274       }
 275       StockholmFile file = new StockholmFile(new FileParse(
 276               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
 277       seqs = file.getSeqsAsArray();
 278
 279       if (searchAlignment)
 280       {
 281         recoverSequences(sequencesHash, seqs);
 282       }
 283
 284       readTable(searchOutputFile);
 285
 286       int seqCount = Math.min(seqs.length, seqsToReturn);
 287
 288       AlignmentI al = new Alignment(seqs);
 289
 290       AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 291               AlignFrame.DEFAULT_HEIGHT);
 292       String ttl = "jackhmmer search of " + databaseName + " using "
 293               + seqs[0].getName();
 294       Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 295               AlignFrame.DEFAULT_HEIGHT);
 296
 297       seqTemp.delete();
 298       inputAlignmentTemp.delete();
 299       searchOutputFile.delete();
 300     } finally
 301     {
 302       if (br != null)
 303       {
 304         br.close();
 305       }
 306     }
 307   }
 308
 309   /**
 310    * Reads in the scores table output by jackhmmer and adds annotation to
 311    * sequences for E-value and bit score
 312    *
 313    * @param inputTableTemp
 314    * @throws IOException
 315    */
 316   void readTable(File inputTableTemp) throws IOException
 317   {
 318     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 319     String line = "";
 320     while (!line.startsWith("Query:"))
 321     {
 322       line = br.readLine();
 323     }
 324     while (!line.contains("-------"))
 325     {
 326       line = br.readLine();
 327     }
 328     line = br.readLine();
 329
 330     int index = 0;
 331     while (!"  ------ inclusion threshold ------".equals(line)
 332             && !"".equals(line))
 333     {
 334       SequenceI seq = seqs[index];
 335
 336       Scanner scanner = new Scanner(line);
 337       String evalue = scanner.next();
 338       evalue = scanner.next();
 339       String score = scanner.next();
 340       addScoreAnnotations(evalue, score, seq);
 341       scanner.close();
 342       line = br.readLine();
 343       index++;
 344     }
 345
 346     br.close();
 347   }
 348
 349   protected void addScoreAnnotations(String eValue, String bitScore,
 350           SequenceI seq)
 351   {
 352     String label = "Search Scores";
 353     String description = "Full sequence bit score and E-Value";
 354
 355     try
 356     {
 357       AlignmentAnnotation annot = new AlignmentAnnotation(label,
 358               description, null);
 359
 360       annot.label = label;
 361       annot.description = description;
 362
 363       annot.setCalcId(JACKHMMER);
 364
 365       double dEValue = Double.parseDouble(eValue);
 366       annot.setEValue(dEValue);
 367
 368       double dBitScore = Double.parseDouble(bitScore);
 369       annot.setBitScore(dBitScore);
 370
 371       annot.setSequenceRef(seq);
 372       seq.addAlignmentAnnotation(annot);
 373
 374     } catch (NumberFormatException e)
 375     {
 376       System.err.println("Error parsing " + label + " from " + eValue
 377               + " & " + bitScore);
 378     }
 379   }
 380
 381
 382
 383 }