src/jalview/hmmer/HMMSearch.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.Annotation;
   8 import jalview.datamodel.HiddenMarkovModel;
   9 import jalview.datamodel.SequenceI;
  10 import jalview.gui.AlignFrame;
  11 import jalview.gui.JvOptionPane;
  12 import jalview.io.DataSourceType;
  13 import jalview.io.FileParse;
  14 import jalview.io.StockholmFile;
  15 import jalview.util.FileUtils;
  16 import jalview.util.MessageManager;
  17 import jalview.ws.params.ArgumentI;
  18 import jalview.ws.params.simple.BooleanOption;
  19 import jalview.ws.params.simple.Option;
  20
  21 import java.io.BufferedReader;
  22 import java.io.File;
  23 import java.io.FileReader;
  24 import java.io.IOException;
  25 import java.util.ArrayList;
  26 import java.util.Collections;
  27 import java.util.List;
  28 import java.util.Scanner;
  29
  30 import javax.swing.JOptionPane;
  31
  32 public class HMMSearch extends HmmerCommand
  33 {
  34   private static final String PARAMNAME_NO_OF_RESULTS = MessageManager.getString("label.number_of_results");
  35
  36   static final String HMMSEARCH = "hmmsearch";
  37
  38   boolean realign = false;
  39
  40   boolean trim = false;
  41
  42   int seqsToReturn = Integer.MAX_VALUE;
  43
  44   SequenceI[] seqs;
  45
  46   /**
  47    * Constructor for the HMMSearchThread
  48    *
  49    * @param af
  50    */
  51   public HMMSearch(AlignFrame af, List<ArgumentI> args)
  52   {
  53     super(af, args);
  54   }
  55
  56   /**
  57    * Runs the HMMSearchThread: the data on the alignment or group is exported,
  58    * then the command is executed in the command line and then the data is
  59    * imported and displayed in a new frame. Call this method directly to execute
  60    * synchronously, or via start() in a new Thread for asynchronously.
  61    */
  62   @Override
  63   public void run()
  64   {
  65     HiddenMarkovModel hmm = getHmmProfile();
  66     if (hmm == null)
  67     {
  68       // shouldn't happen if we got this far
  69       Cache.log.error("Error: no hmm for hmmsearch");
  70       return;
  71     }
  72
  73     SequenceI hmmSeq = hmm.getConsensusSequence();// af.getSelectedHMMSequence();
  74     long msgId = System.currentTimeMillis();
  75     af.setProgressBar(MessageManager.getString("status.running_hmmsearch"),
  76             msgId);
  77
  78     try
  79     {
  80       File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
  81       File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
  82               ".sto");
  83       File searchOutputFile = FileUtils.createTempFile("searchOutput",
  84               ".sto");
  85
  86       exportHmm(hmm, hmmFile.getAbsoluteFile());
  87
  88       boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
  89       if (!ran)
  90       {
  91         JvOptionPane.showInternalMessageDialog(af, MessageManager
  92                 .formatMessage("warn.command_failed", "hmmsearch"));
  93         return;
  94       }
  95
  96       importData(hmmSeq, hitsAlignmentFile, hmmFile, 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     }
 103     finally
 104     {
 105       af.setProgressBar("", msgId);
 106     }
 107   }
 108
 109   /**
 110    * Executes an hmmsearch with the given hmm as input. The database to be
 111    * searched is a local file as specified by the 'Database' parameter, or the
 112    * current alignment (written to file) if none is specified.
 113    *
 114    * @param searchOutputFile
 115    * @param hitsAlignmentFile
 116    * @param hmmFile
 117    *
 118    * @return
 119    * @throws IOException
 120    */
 121   private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
 122           File hmmFile) throws IOException
 123   {
 124     String command = getCommandPath(HMMSEARCH);
 125     if (command == null)
 126     {
 127       return false;
 128     }
 129
 130     List<String> args = new ArrayList<>();
 131     args.add(command);
 132     args.add("-o");
 133     args.add(getFilePath(searchOutputFile));
 134     args.add("-A");
 135     args.add(getFilePath(hitsAlignmentFile));
 136
 137     boolean dbFound = false;
 138     String dbPath = "";
 139     File databaseFile = null;
 140
 141     if (params != null)
 142     {
 143       for (ArgumentI arg : params)
 144       {
 145         String name = arg.getName();
 146         if (MessageManager.getString("label.number_of_results")
 147                 .equals(name))
 148         {
 149           seqsToReturn = Integer.parseInt(arg.getValue());
 150         }
 151         else if (MessageManager.getString("label.auto_align_seqs")
 152                 .equals(name))
 153         {
 154           realign = true; // TODO: not used
 155         }
 156         else if (MessageManager.getString("label.use_accessions")
 157                 .equals(name))
 158         {
 159           args.add("--acc");
 160         }
 161         else if (MessageManager.getString("label.seq_e_value").equals(name))
 162         {
 163           args.add("-E");
 164           args.add(arg.getValue());
 165         }
 166         else if (MessageManager.getString("label.seq_score").equals(name))
 167         {
 168           args.add("-incT");
 169           args.add(arg.getValue());
 170         }
 171         else if (MessageManager.getString("label.dom_e_value")
 172                 .equals(name))
 173         {
 174           args.add("--domE");
 175           args.add(arg.getValue());
 176         }
 177         else if (MessageManager.getString("label.dom_score").equals(name))
 178         {
 179           args.add("--incdomT");
 180           args.add(arg.getValue());
 181         }
 182         else if (MessageManager.getString("label.trim_termini")
 183                 .equals(name))
 184         {
 185           trim = true;
 186         }
 187         else if (MessageManager.getString("label.database").equals(name))
 188         {
 189           dbFound = true;
 190           dbPath = arg.getValue();
 191           if (!MessageManager.getString("label.this_alignment")
 192                   .equals(dbPath))
 193           {
 194             databaseFile = new File(dbPath);
 195           }
 196         }
 197       }
 198     }
 199
 200     if (!dbFound || MessageManager.getString("label.this_alignment")
 201             .equals(dbPath))
 202     {
 203       /*
 204        * no external database specified for search, so
 205        * export current alignment as 'database' to search,
 206        * excluding any HMM consensus sequences it contains
 207        */
 208       databaseFile = FileUtils.createTempFile("database", ".sto");
 209       AlignmentI al = af.getViewport().getAlignment();
 210       AlignmentI copy = new Alignment(al);
 211       List<SequenceI> hmms = copy.getHmmSequences();
 212       for (SequenceI hmmSeq : hmms)
 213       {
 214         copy.deleteSequence(hmmSeq);
 215       }
 216       exportStockholm(copy.getSequencesArray(), databaseFile, null);
 217       // StockholmFile stoFile = new StockholmFile(copy);
 218       // stoFile.setSeqs(copy.getSequencesArray());
 219       // String alignmentString = stoFile.print();
 220       // PrintWriter writer = new PrintWriter(databaseFile);
 221       // writer.print(alignmentString);
 222       // writer.close();
 223     }
 224
 225     args.add(getFilePath(hmmFile));
 226     args.add(getFilePath(databaseFile));
 227
 228     return runCommand(args);
 229   }
 230
 231   /**
 232    * Imports the data from the temporary file to which the output of hmmsearch
 233    * is directed.
 234    *
 235    * @param hmmSeq
 236    */
 237   private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
 238           File hmmTemp, File searchOutputFile)
 239           throws IOException, InterruptedException
 240   {
 241     BufferedReader br = new BufferedReader(
 242             new FileReader(inputAlignmentTemp));
 243     try
 244     {
 245       if (br.readLine() == null)
 246       {
 247         JOptionPane.showMessageDialog(af,
 248                 MessageManager.getString("label.no_sequences_found"));
 249         return;
 250       }
 251       StockholmFile file = new StockholmFile(new FileParse(
 252               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
 253       seqs = file.getSeqsAsArray();
 254
 255       readTable(searchOutputFile);
 256
 257       int seqCount = Math.min(seqs.length, seqsToReturn);
 258       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 259       hmmAndSeqs[0] = hmmSeq;
 260       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 261
 262       /*
 263        * and align the search results to the HMM profile
 264        */
 265       AlignmentI al = new Alignment(hmmAndSeqs);
 266       AlignFrame frame = new AlignFrame(al, 1, 1);
 267       List<ArgumentI> alignArgs = new ArrayList<>();
 268       String defSeq = hmmSeq.getName();
 269       List<String> options = Collections.singletonList(defSeq);
 270       Option option = new Option(MessageManager.getString("label.use_hmm"),
 271               "", true, defSeq, defSeq, options, null);
 272       alignArgs.add(option);
 273       if (trim)
 274       {
 275         alignArgs.add(new BooleanOption(
 276                 MessageManager.getString("label.trim_termini"),
 277                 MessageManager.getString("label.trim_termini_desc"), true,
 278                 true, true, null));
 279       }
 280       HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 281       hmmalign.run();
 282       frame = null;
 283       hmmTemp.delete();
 284       inputAlignmentTemp.delete();
 285       searchOutputFile.delete();
 286     } finally
 287     {
 288       if (br != null)
 289       {
 290         br.close();
 291       }
 292     }
 293   }
 294
 295   void readTable(File inputTableTemp) throws IOException
 296   {
 297     BufferedReader br = new BufferedReader(new FileReader(inputTableTemp));
 298     String line = "";
 299     while (!line.startsWith("Query:"))
 300     {
 301       line = br.readLine();
 302     }
 303     for (int i = 0; i < 5; i++)
 304     {
 305       line = br.readLine();
 306     }
 307
 308     int index = 0;
 309     while (!"  ------ inclusion threshold ------".equals(line)
 310             && !"".equals(line))
 311     {
 312       Scanner scanner = new Scanner(line);
 313
 314       String str = scanner.next(); // full sequence eValue score
 315       float eValue = Float.parseFloat(str);
 316       int seqLength = seqs[index].getLength();
 317       Annotation[] annots = new Annotation[seqLength];
 318       for (int j = 0; j < seqLength; j++)
 319       {
 320         annots[j] = new Annotation(eValue);
 321       }
 322       AlignmentAnnotation annot = new AlignmentAnnotation("E-value",
 323               "Score", annots);
 324       annot.setScore(Double.parseDouble(str));
 325       annot.setSequenceRef(seqs[index]);
 326       seqs[index].addAlignmentAnnotation(annot);
 327
 328       scanner.close();
 329       line = br.readLine();
 330       index++;
 331     }
 332
 333     br.close();
 334   }
 335
 336 }