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
 332       Scanner scanner = new Scanner(line);
 333       String evalue = scanner.next();
 334       evalue = scanner.next();
 335       checkSequenceOrder(index, scanner);
 336       SequenceI seq = seqs[index];
 337       String score = scanner.next();
 338       addScoreAnnotations(evalue, score, seq);
 339       scanner.close();
 340       line = br.readLine();
 341       index++;
 342     }
 343
 344     br.close();
 345   }
 346
 347   private void checkSequenceOrder(int index, Scanner scanner)
 348   {
 349     String seqName = null;
 350
 351     for (int i = 0; i < 8; i++)
 352     {
 353       seqName = scanner.next();
 354     }
 355
 356     if (!seqs[index + 1].getName().equals(seqName))
 357     {
 358       SequenceI temp = seqs[index + 1];
 359
 360       for (int j = 0; j < seqs.length; j++)
 361       {
 362         if (seqs[j].getName().equals(seqName))
 363         {
 364           seqs[index + 1] = seqs[j];
 365           seqs[j] = temp;
 366           break;
 367         }
 368       }
 369     }
 370   }
 371
 372   protected void addScoreAnnotations(String eValue, String bitScore,
 373           SequenceI seq)
 374   {
 375     String label = "Search Scores";
 376     String description = "Full sequence bit score and E-Value";
 377
 378     try
 379     {
 380       AlignmentAnnotation annot = new AlignmentAnnotation(label,
 381               description, null);
 382
 383       annot.label = label;
 384       annot.description = description;
 385
 386       annot.setCalcId(JACKHMMER);
 387
 388       double dEValue = Double.parseDouble(eValue);
 389       annot.setEValue(dEValue);
 390
 391       double dBitScore = Double.parseDouble(bitScore);
 392       annot.setBitScore(dBitScore);
 393
 394       annot.setSequenceRef(seq);
 395       seq.addAlignmentAnnotation(annot);
 396
 397     } catch (NumberFormatException e)
 398     {
 399       System.err.println("Error parsing " + label + " from " + eValue
 400               + " & " + bitScore);
 401     }
 402   }
 403
 404
 405
 406 }