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