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