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