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