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       readDomainTable(searchOutputFile, false);
 161
 162       if (searchAlignment)
 163       {
 164         recoverSequences(sequencesHash, seqs);
 165       }
 166
 167       // look for PP cons and ref seq in alignment only annotation
 168       AlignmentAnnotation modelpos = null, ppcons = null;
 169       for (AlignmentAnnotation aa : file.getAnnotations())
 170       {
 171         if (aa.sequenceRef == null)
 172         {
 173           if (aa.label.equals("Reference Positions")) // RF feature type in
 174                                                       // stockholm parser
 175           {
 176             modelpos = aa;
 177           }
 178           if (aa.label.equals("Posterior Probability"))
 179           {
 180             ppcons = aa;
 181           }
 182         }
 183       }
 184
 185
 186       int seqCount = Math.min(seqs.length, seqsToReturn);
 187       SequenceI[] hmmAndSeqs = new SequenceI[seqCount + 1];
 188       hmmSeq = hmmSeq.deriveSequence(); // otherwise all bad things happen
 189       hmmAndSeqs[0] = hmmSeq;
 190       System.arraycopy(seqs, 0, hmmAndSeqs, 1, seqCount);
 191       if (modelpos != null)
 192       {
 193         // TODO need - get ungapped sequence method
 194         hmmSeq.setSequence(
 195                 hmmSeq.getDatasetSequence().getSequenceAsString());
 196         Annotation[] refpos = modelpos.annotations;
 197         // insert gaps to match with refseq positions
 198         int gc = 0, lcol = 0;
 199         for (int c = 0; c < refpos.length; c++)
 200         {
 201           if (refpos[c] != null && ("x".equals(refpos[c].displayCharacter)))
 202           {
 203             if (gc > 0)
 204             {
 205               hmmSeq.insertCharAt(lcol + 1, gc, '-');
 206             }
 207             gc = 0;
 208             lcol = c;
 209           }
 210           else
 211           {
 212             gc++;
 213           }
 214         }
 215       }
 216
 217       if (realign)
 218       {
 219         realignResults(hmmAndSeqs);
 220       }
 221       else
 222       {
 223         AlignmentI al = new Alignment(hmmAndSeqs);
 224         if (ppcons != null)
 225         {
 226           al.addAnnotation(ppcons);
 227         }
 228         if (modelpos != null)
 229         {
 230           al.addAnnotation(modelpos);
 231         }
 232         AlignFrame alignFrame = new AlignFrame(al, AlignFrame.DEFAULT_WIDTH,
 233                 AlignFrame.DEFAULT_HEIGHT);
 234         String ttl = "hmmSearch of " + databaseName + " using "
 235                 + hmmSeq.getName();
 236         Desktop.addInternalFrame(alignFrame, ttl, AlignFrame.DEFAULT_WIDTH,
 237                 AlignFrame.DEFAULT_HEIGHT);
 238
 239         if (returnNoOfNewSeqs)
 240         {
 241           int nNew = checkForNewSequences();
 242           JvOptionPane.showMessageDialog(af.alignPanel, nNew + " "
 243                   + MessageManager.getString("label.new_returned"));
 244         }
 245
 246       }
 247
 248
 249       hmmTemp.delete();
 250       inputAlignmentTemp.delete();
 251       searchOutputFile.delete();
 252     } finally
 253     {
 254       if (br != null)
 255       {
 256         br.close();
 257       }
 258     }
 259   }
 260
 261   private int checkForNewSequences()
 262   {
 263     int nNew = seqs.length;
 264
 265     for (SequenceI resultSeq : seqs)
 266     {
 267       for (SequenceI aliSeq : alignment.getSequencesArray())
 268       {
 269         if (resultSeq.getName().equals(aliSeq.getName()))
 270         {
 271           nNew--;
 272           break;
 273         }
 274       }
 275     }
 276
 277     return nNew;
 278
 279   }
 280
 281   /**
 282    * Realigns the given sequences using hmmalign, to the HMM profile sequence
 283    * which is the first in the array, and opens the results in a new frame
 284    *
 285    * @param hmmAndSeqs
 286    */
 287   protected void realignResults(SequenceI[] hmmAndSeqs)
 288   {
 289     /*
 290      * and align the search results to the HMM profile
 291      */
 292     AlignmentI al = new Alignment(hmmAndSeqs);
 293     AlignFrame frame = new AlignFrame(al, 1, 1);
 294     List<ArgumentI> alignArgs = new ArrayList<>();
 295     String alignTo = hmmAndSeqs[0].getName();
 296     List<String> options = Collections.singletonList(alignTo);
 297     Option option = new Option(MessageManager.getString("label.use_hmm"),
 298             "", true, alignTo, alignTo, options, null);
 299     alignArgs.add(option);
 300     if (trim)
 301     {
 302       alignArgs.add(new BooleanOption(
 303               MessageManager.getString(TRIM_TERMINI_KEY),
 304               MessageManager.getString("label.trim_termini_desc"), true,
 305               true, true, null));
 306     }
 307     HmmerCommand hmmalign = new HMMAlign(frame, alignArgs);
 308     hmmalign.run();
 309
 310     if (returnNoOfNewSeqs)
 311     {
 312       int nNew = checkForNewSequences();
 313       JvOptionPane.showMessageDialog(frame.alignPanel,
 314               nNew + " " + MessageManager.getString("label.new_returned"));
 315     }
 316   }
 317
 318 }