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       AlignmentI copy = new Alignment(al);
 250       deleteHmmSequences(copy);
 251       exportFasta(copy, databaseFile);
 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       readTable(searchOutputFile);
 280
 281       int seqCount = Math.min(seqs.length, seqsToReturn);
 282
 283       AlignmentI al = new Alignment(seqs);
 284
 285       AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 286               AlignFrame.DEFAULT_HEIGHT);
 287       String ttl = "jackhmmer search of " + databaseName + " using "
 288               + seqs[0].getName();
 289       Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 290               AlignFrame.DEFAULT_HEIGHT);
 291
 292       seqTemp.delete();
 293       inputAlignmentTemp.delete();
 294       searchOutputFile.delete();
 295     } finally
 296     {
 297       if (br != null)
 298       {
 299         br.close();
 300       }
 301     }
 302   }
 303
 304   /**
 305    * Reads in the scores table output by jackhmmer and adds annotation to
 306    * sequences for E-value and bit score
 307    *
 308    * @param inputTableTemp
 309    * @throws IOException
 310    */
 311   void readTable(File inputTableTemp) throws IOException
 312   {
 313     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 314     String line = "";
 315     while (!line.startsWith("Query:"))
 316     {
 317       line = br.readLine();
 318     }
 319     while (!line.contains("-------"))
 320     {
 321       line = br.readLine();
 322     }
 323     line = br.readLine();
 324
 325     int index = 0;
 326     while (!"  ------ inclusion threshold ------".equals(line)
 327             && !"".equals(line))
 328     {
 329       SequenceI seq = seqs[index];
 330
 331       Scanner scanner = new Scanner(line);
 332       String evalue = scanner.next();
 333       evalue = scanner.next();
 334       String score = scanner.next();
 335       addScoreAnnotations(evalue, score, seq);
 336       scanner.close();
 337       line = br.readLine();
 338       index++;
 339     }
 340
 341     br.close();
 342   }
 343
 344   protected void addScoreAnnotations(String eValue, String bitScore,
 345           SequenceI seq)
 346   {
 347     String label = "Search Scores";
 348     String description = "Full sequence bit score and E-Value";
 349
 350     try
 351     {
 352       AlignmentAnnotation annot = new AlignmentAnnotation(label,
 353               description, null);
 354
 355       annot.label = label;
 356       annot.description = description;
 357
 358       annot.setCalcId(JACKHMMER);
 359
 360       double dEValue = Double.parseDouble(eValue);
 361       annot.setEValue(dEValue);
 362
 363       double dBitScore = Double.parseDouble(bitScore);
 364       annot.setEValue(dBitScore);
 365
 366       annot.setSequenceRef(seq);
 367       seq.addAlignmentAnnotation(annot);
 368     } catch (NumberFormatException e)
 369     {
 370       System.err.println("Error parsing " + label + " from " + eValue
 371               + " & " + bitScore);
 372     }
 373   }
 374
 375
 376
 377 }