src/jalview/hmmer/HMMSearch.java

   1 package jalview.hmmer;
   2
   3 import jalview.bin.Cache;
   4 import jalview.bin.Console;
   5 import jalview.datamodel.Alignment;
   6 import jalview.datamodel.AlignmentAnnotation;
   7 import jalview.datamodel.AlignmentI;
   8 import jalview.datamodel.Annotation;
   9 import jalview.datamodel.HiddenMarkovModel;
  10 import jalview.datamodel.SequenceI;
  11 import jalview.gui.AlignFrame;
  12 import jalview.gui.Desktop;
  13 import jalview.gui.JvOptionPane;
  14 import jalview.io.DataSourceType;
  15 import jalview.io.FileParse;
  16 import jalview.io.StockholmFile;
  17 import jalview.util.FileUtils;
  18 import jalview.util.MessageManager;
  19 import jalview.ws.params.ArgumentI;
  20 import jalview.ws.params.simple.BooleanOption;
  21 import jalview.ws.params.simple.Option;
  22
  23 import java.io.BufferedReader;
  24 import java.io.File;
  25 import java.io.FileReader;
  26 import java.io.IOException;
  27 import java.util.ArrayList;
  28 import java.util.Collections;
  29 import java.util.List;
  30
  31 import javax.swing.JOptionPane;
  32
  33 public class HMMSearch extends Search
  34 {
  35
  36   boolean realign = false;
  37
  38   boolean trim = false;
  39
  40   boolean returnNoOfNewSeqs = false;
  41
  42   int seqsToReturn = Integer.MAX_VALUE;
  43
  44
  45   /**
  46    * Constructor for the HMMSearchThread
  47    *
  48    * @param af
  49    */
  50   public HMMSearch(AlignFrame af, List<ArgumentI> args)
  51   {
  52     super(af, args);
  53   }
  54
  55   /**
  56    * Runs the HMMSearchThread: 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     HiddenMarkovModel hmm = getHmmProfile();
  65     if (hmm == null)
  66     {
  67       // shouldn't happen if we got this far
  68       Console.error("Error: no hmm for hmmsearch");
  69       return;
  70     }
  71
  72     SequenceI hmmSeq = hmm.getConsensusSequence();
  73     long msgId = System.currentTimeMillis();
  74     af.setProgressBar(MessageManager.getString("status.running_search"),
  75             msgId);
  76
  77     try
  78     {
  79       File hmmFile = FileUtils.createTempFile("hmm", ".hmm");
  80       File hitsAlignmentFile = FileUtils.createTempFile("hitAlignment",
  81               ".sto");
  82       File searchOutputFile = FileUtils.createTempFile("searchOutput",
  83               ".sto");
  84
  85       exportHmm(hmm, hmmFile.getAbsoluteFile());
  86
  87       boolean ran = runCommand(searchOutputFile, hitsAlignmentFile, hmmFile);
  88       if (!ran)
  89       {
  90         JvOptionPane.showInternalMessageDialog(af, MessageManager
  91                 .formatMessage("warn.command_failed", "hmmsearch"));
  92         return;
  93       }
  94
  95       importData(hmmSeq, hitsAlignmentFile, hmmFile, searchOutputFile);
  96       // TODO make realignment of search results a step at this level
  97       // and make it conditional on this.realign
  98     } catch (IOException | InterruptedException e)
  99     {
 100       e.printStackTrace();
 101     }
 102     finally
 103     {
 104       af.setProgressBar("", msgId);
 105     }
 106   }
 107
 108   /**
 109    * Executes an hmmsearch with the given hmm as input. The database to be
 110    * searched is a local file as specified by the 'Database' parameter, or the
 111    * current alignment (written to file) if none is specified.
 112    *
 113    * @param searchOutputFile
 114    * @param hitsAlignmentFile
 115    * @param hmmFile
 116    *
 117    * @return
 118    * @throws IOException
 119    */
 120   private boolean runCommand(File searchOutputFile, File hitsAlignmentFile,
 121           File hmmFile) throws IOException
 122   {
 123     String command = getCommandPath(HMMSEARCH);
 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, hmmFile);
 132
 133     return runCommand(args);
 134   }
 135
 136
 137   /**
 138    * Imports the data from the temporary file to which the output of hmmsearch
 139    * was directed. The results are optionally realigned using hmmalign.
 140    *
 141    * @param hmmSeq
 142    */
 143   private void importData(SequenceI hmmSeq, File inputAlignmentTemp,
 144           File hmmTemp, File searchOutputFile)
 145           throws IOException, InterruptedException
 146   {
 147     BufferedReader br = new BufferedReader(
 148             new FileReader(inputAlignmentTemp));
 149     try
 150     {
 151       if (br.readLine() == null)
 152       {
 153         JOptionPane.showMessageDialog(af,
 154                 MessageManager.getString("label.no_sequences_found"));
 155         return;
 156       }
 157       StockholmFile file = new StockholmFile(new FileParse(
 158               inputAlignmentTemp.getAbsolutePath(), DataSourceType.FILE));
 159       seqs = file.getSeqsAsArray();
 160
 161       readDomainTable(searchOutputFile, false);
 162
 163       if (searchAlignment)
 164       {
 165         recoverSequences(sequencesHash, seqs);
 166       }
 167
 168       // look for PP cons and ref seq in alignment only annotation
 169       AlignmentAnnotation modelpos = null, ppcons = null;
 170       for (AlignmentAnnotation aa : file.getAnnotations())
 171       {
 172         if (aa.sequenceRef == null)
 173         {
 174           if (aa.label.equals("Reference Positions")) // RF feature type in
 175                                                       // stockholm parser
 176           {
 177             modelpos = aa;
 178           }
 179           if (aa.label.equals("Posterior Probability"))
 180           {
 181             ppcons = aa;
 182           }
 183         }
 184       }
 185
 186
 187       int seqCount = Math.min(seqs.length, seqsToReturn);
 188       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 189       hmmSeq = hmmSeq.deriveSequence(); // otherwise all bad things happen
 190       hmmAndSeqs[0] = hmmSeq;
 191       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 192       if (modelpos != null)
 193       {
 194         // TODO need - get ungapped sequence method
 195         hmmSeq.setSequence(
 196                 hmmSeq.getDatasetSequence().getSequenceAsString());
 197         Annotation[] refpos = modelpos.annotations;
 198         // insert gaps to match with refseq positions
 199         int gc = 0, lcol = 0;
 200         for (int c = 0; c < refpos.length; c++)
 201         {
 202           if (refpos[c] != null && ("x".equals(refpos[c].displayCharacter)))
 203           {
 204             if (gc > 0)
 205             {
 206               hmmSeq.insertCharAt(lcol + 1, gc, '-');
 207             }
 208             gc = 0;
 209             lcol = c;
 210           }
 211           else
 212           {
 213             gc++;
 214           }
 215         }
 216       }
 217
 218       if (realign)
 219       {
 220         realignResults(hmmAndSeqs);
 221       }
 222       else
 223       {
 224         AlignmentI al = new Alignment(hmmAndSeqs);
 225         if (ppcons != null)
 226         {
 227           al.addAnnotation(ppcons);
 228         }
 229         if (modelpos != null)
 230         {
 231           al.addAnnotation(modelpos);
 232         }
 233         AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 234                 AlignFrame.DEFAULT_HEIGHT);
 235         String ttl = "hmmSearch of " + databaseName + " using "
 236                 + hmmSeq.getName();
 237         Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 238                 AlignFrame.DEFAULT_HEIGHT);
 239
 240         if (returnNoOfNewSeqs)
 241         {
 242           int nNew = checkForNewSequences();
 243           JvOptionPane.showMessageDialog(af.alignPanel, nNew + " "
 244                   + MessageManager.getString("label.new_returned"));
 245         }
 246
 247       }
 248
 249
 250       hmmTemp.delete();
 251       inputAlignmentTemp.delete();
 252       searchOutputFile.delete();
 253     } finally
 254     {
 255       if (br != null)
 256       {
 257         br.close();
 258       }
 259     }
 260   }
 261
 262   private int checkForNewSequences()
 263   {
 264     int nNew = seqs.length;
 265
 266     for (SequenceI resultSeq : seqs)
 267     {
 268       for (SequenceI aliSeq : alignment.getSequencesArray())
 269       {
 270         if (resultSeq.getName().equals(aliSeq.getName()))
 271         {
 272           nNew--;
 273           break;
 274         }
 275       }
 276     }
 277
 278     return nNew;
 279
 280   }
 281
 282   /**
 283    * Realigns the given sequences using hmmalign, to the HMM profile sequence
 284    * which is the first in the array, and opens the results in a new frame
 285    *
 286    * @param hmmAndSeqs
 287    */
 288   protected void realignResults(SequenceI[] hmmAndSeqs)
 289   {
 290     /*
 291      * and align the search results to the HMM profile
 292      */
 293     AlignmentI al = new Alignment(hmmAndSeqs);
 294     AlignFrame frame = new AlignFrame(al, 1, 1);
 295     List<ArgumentI> alignArgs = new ArrayList<>();
 296     String alignTo = hmmAndSeqs[0].getName();
 297     List<String> options = Collections.singletonList(alignTo);
 298     Option option = new Option(MessageManager.getString("label.use_hmm"),
 299             "", true, alignTo, alignTo, options, null);
 300     alignArgs.add(option);
 301     if (trim)
 302     {
 303       alignArgs.add(new BooleanOption(
 304               MessageManager.getString(TRIM_TERMINI_KEY),
 305               MessageManager.getString("label.trim_termini_desc"), true,
 306               true, true, null));
 307     }
 308     HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 309     hmmalign.run();
 310
 311     if (returnNoOfNewSeqs)
 312     {
 313       int nNew = checkForNewSequences();
 314       JvOptionPane.showMessageDialog(frame.alignPanel,
 315               nNew + " " + MessageManager.getString("label.new_returned"));
 316     }
 317   }
 318
 319 }