src/jalview/schemes/HMMMatchScoreColourScheme.java

   1 package jalview.schemes;
   2
   3 import jalview.api.AlignViewportI;
   4 import jalview.datamodel.AnnotatedCollectionI;
   5 import jalview.datamodel.HiddenMarkovModel;
   6 import jalview.datamodel.SequenceI;
   7 import jalview.util.ColorUtils;
   8 import jalview.util.Comparison;
   9
  10 import java.awt.Color;
  11 import java.io.BufferedReader;
  12 import java.io.FileNotFoundException;
  13 import java.io.FileReader;
  14 import java.io.IOException;
  15 import java.util.ArrayList;
  16 import java.util.HashMap;
  17 import java.util.List;
  18 import java.util.Map;
  19
  20
  21
  22
  23 public class HMMMatchScoreColourScheme extends ResidueColourScheme
  24 {
  25
  26   private Map<Character, Map<Integer, Map<String, Double>>> probabilities;
  27
  28   private List<Integer> ranges;
  29
  30   private static double binSize;
  31
  32   public class MatchProbReader
  33   {
  34     private BufferedReader reader;
  35
  36     MatchProbReader() throws FileNotFoundException
  37     {
  38       reader = new BufferedReader(
  39               new FileReader("resources/ProbabilityOfMatch"));
  40     }
  41
  42     /*
  43     public Map<Character, Map<String, Double>>  getProbabilities() throws IOException
  44     {
  45       Map<Character, Map<String, Double>> probabilities = new HashMap<>();
  46
  47       String[] alphabet = reader.readLine().split("\\,");
  48       for (int i = 1; i < alphabet.length - 1; i++)
  49       {
  50         probabilities.put(alphabet[i].replaceAll("\\ ", "").charAt(0),
  51                 new HashMap<>());
  52       }
  53
  54       String line = reader.readLine();
  55       while(line != null)
  56       {
  57         String[] contents = line.split("\\,");
  58
  59         for(int i = 1; i < contents.length; i++)
  60         {
  61           probabilities.get(alphabet[i].replaceAll("\\ ", "").charAt(0))
  62                   .put(contents[0], Double
  63                           .valueOf(contents[i].replaceAll("\\ ", "")));
  64         }
  65         line = reader.readLine();
  66
  67       }
  68       reader.close();
  69       return probabilities;
  70     }
  71     */
  72
  73     public Map<Character, Map<Integer, Map<String, Double>>> getProbabilities()
  74             throws IOException
  75     {
  76
  77       Map<Character, Map<Integer, Map<String, Double>>> probabilities = new HashMap<>();
  78
  79       ranges = new ArrayList<>();
  80       ranges.add(0);
  81
  82       binSize = Double.valueOf((reader.readLine().replaceAll("\\ ", "")));
  83       String line = reader.readLine();
  84       char c = line.charAt(0);
  85
  86       while (line != null)
  87       {
  88         line = reader.readLine();
  89         while (line != null && line.split("\\,").length != 1)
  90         {
  91           String[] llrs = line.split("\\,");
  92           String[] counts = reader.readLine().split("\\,");
  93           int range = Integer.valueOf(llrs[0]);
  94
  95           if (!ranges.contains(range))
  96           {
  97             ranges.add(range);
  98           }
  99           if (!probabilities.containsKey(c))
 100           {
 101             probabilities.put(c, new HashMap<>());
 102           }
 103           probabilities.get(c).put(range, new HashMap<>());
 104
 105           for (int i = 1; i < llrs.length; i++)
 106           {
 107             probabilities.get(c).get(range).put(
 108                     llrs[i].replaceAll("\\ ", ""),
 109                     Double.valueOf(counts[i].replaceAll("\\ ", "")));
 110           }
 111
 112           line = reader.readLine();
 113         }
 114         if (line != null)
 115         {
 116           c = line.charAt(0);
 117         }
 118       }
 119
 120       return probabilities;
 121     }
 122
 123   }
 124
 125   private static final Color INSERTION_COLOUR = Color.white;
 126
 127   /*
 128    * the aligned HMM consensus sequence to use as reference for colouring
 129    */
 130   private SequenceI hmmSeq;
 131
 132   private HiddenMarkovModel hmm;
 133
 134   private Map<Character, Float> frequencies;
 135
 136   /**
 137    * Constructor given a list of Hidden Markov Model consensus sequences. The
 138    * first sequence provides the HMM profile from which we can read the emission
 139    * probabilities that determine the colour.
 140    *
 141    * @param hmmSeqs
 142    * @throws IOException
 143    */
 144   public HMMMatchScoreColourScheme(List<SequenceI> hmmSeqs)
 145   {
 146     hmmSeq = hmmSeqs.isEmpty() ? null : hmmSeqs.get(0);
 147     hmm = hmmSeq == null ? null : hmmSeq.getHMM();
 148
 149     try
 150     {
 151       MatchProbReader probabilityReader = new MatchProbReader();
 152       probabilities = probabilityReader.getProbabilities();
 153     } catch (IOException e)
 154     {
 155       System.out.println(e.getStackTrace());
 156     }
 157   }
 158
 159   /**
 160    * Default constructor (required by ColourSchemes.loadColourSchemes)
 161    */
 162   public HMMMatchScoreColourScheme()
 163   {
 164   }
 165
 166   @Override
 167   public Color findColour(char symbol, int column, SequenceI seq,
 168           String consensusResidue, float pid)
 169   {
 170     if (seq == null)
 171     {
 172       return null;
 173     }
 174     return findColour(symbol, column, seq.gapMap().length);
 175   }
 176
 177   // TODO change documentation
 178   /**
 179    * Returns the colour at a particular symbol at a column in the alignment:
 180    * <ul>
 181    * <li>white for a gap</li>
 182    * <li>red for an insertion</li>
 183    * <li>orange for negative information content</li>
 184    * <li>white to blue for increasing information content</li>
 185    * </ul>
 186    *
 187    * @param symbol
 188    * @param column
 189    * @return
 190    */
 191   private Color findColour(char symbol, int column, int length)
 192   {
 193     if (getHmm() == null || Comparison.isGap(symbol))
 194     {
 195       return Color.white;
 196     }
 197     if (Comparison.isGap(hmmSeq.getCharAt(column)))
 198     {
 199       return INSERTION_COLOUR;
 200     }
 201     if (Character.isLowerCase(symbol))
 202     {
 203       symbol = Character.toUpperCase(symbol);
 204     }
 205
 206     double prob = 0;
 207     if (hmm.getBackgroundFrequencies().containsKey(symbol))
 208     {
 209       int lengthBin = getLengthBin(length);
 210
 211       double llr = Math
 212               .log(getHmm().getMatchEmissionProbability(column, symbol)
 213                       / hmm.getBackgroundFrequencies().get(symbol));
 214
 215       if (!probabilities.containsKey(symbol)
 216               || !probabilities.get(symbol).get(lengthBin)
 217               .containsKey(format(llr)))
 218       {
 219         return new Color(140, 140, 140);
 220       }
 221
 222
 223       prob = probabilities.get(symbol).get(lengthBin).get(format(llr));
 224     }
 225     else
 226     {
 227       return new Color(140, 140, 140);
 228     }
 229
 230     Color colour = Color.ORANGE;
 231     if (prob >= 0.5)
 232     {
 233
 234       colour = ColorUtils.getGraduatedColour((float) prob, 0.5f,
 235               Color.WHITE, 1f,
 236               Color.blue);
 237     }
 238     else
 239     {
 240       colour = ColorUtils.getGraduatedColour((float) prob, 0f, Color.red,
 241               0.5f, Color.WHITE);
 242     }
 243
 244     return colour;
 245   }
 246
 247   public static String format(Double d)
 248   {
 249     String formatArg = String.valueOf(binSize);
 250
 251     // if bin size, need format "%.n" where n is number of decimal places
 252     if (binSize < 1)
 253     {
 254       formatArg = "." + formatArg.split("\\.")[1].length();
 255     }
 256
 257     Double rounded = Math.round(d / binSize) * binSize;
 258     String formatted = String.format("%" + formatArg + "f", rounded);
 259
 260     // format sometimes returns a number rounded to 0 as -0
 261     // this ensures output will always be 0
 262     if (Double.valueOf(formatted) == 0)
 263     {
 264       formatted = "0";
 265     }
 266     return formatted;
 267
 268   }
 269
 270   /**
 271    * Answers the maximum possible value of information score (log ratio), for use
 272    * in scaling a graduated colour range
 273    *
 274    * @return
 275    */
 276   protected float getMaxInformationScore()
 277   {
 278     return 0f;
 279   }
 280
 281   /**
 282    * Answers a new colour scheme instance based on the HMM of the first sequence
 283    * in ac that has an HMM
 284    */
 285   @Override
 286   public ColourSchemeI getInstance(AlignViewportI viewport,
 287           AnnotatedCollectionI ac)
 288   {
 289     return newInstance(ac);
 290   }
 291
 292   /**
 293    * Constructor given a sequence collection
 294    *
 295    * @param ac
 296    */
 297   public HMMMatchScoreColourScheme(AnnotatedCollectionI ac)
 298   {
 299     this(ac.getHmmSequences());
 300   }
 301
 302   /**
 303    * Answers a new instance of the colour scheme for the given HMM
 304    *
 305    * @param ac
 306    * @return
 307    */
 308   protected HMMMatchScoreColourScheme newInstance(AnnotatedCollectionI ac)
 309   {
 310     return new HMMMatchScoreColourScheme(ac);
 311   }
 312
 313   @Override
 314   public boolean isSimple()
 315   {
 316     return false;
 317   }
 318
 319   /**
 320    * Answers true if the sequence collection has an HMM consensus sequence and
 321    * that the first HMM sequence contains background frequencies, else false
 322    */
 323   @Override
 324   public boolean isApplicableTo(AnnotatedCollectionI ac)
 325   {
 326     return !ac.getHmmSequences().isEmpty() && ac.getHmmSequences().get(0)
 327             .getHMM().getBackgroundFrequencies() != null;
 328   }
 329
 330   protected Map<Character, Float> getFrequencies()
 331   {
 332     return frequencies;
 333   }
 334
 335   protected void setFrequencies(Map<Character, Float> frequencies)
 336   {
 337     this.frequencies = frequencies;
 338   }
 339
 340   protected HiddenMarkovModel getHmm()
 341   {
 342     return hmm;
 343   }
 344
 345   protected SequenceI getHmmSequence()
 346   {
 347     return hmmSeq;
 348   }
 349
 350   @Override
 351   public String getSchemeName()
 352   {
 353     return JalviewColourScheme.HMMMatchScore.toString();
 354   }
 355
 356   private int getLengthBin(int l)
 357   {
 358     for (int i = 1; i < ranges.size(); i++)
 359     {
 360       if (l >= ranges.get(i - 1) && l < ranges.get(i))
 361       {
 362         return ranges.get(i);
 363       }
 364     }
 365     return -1;
 366   }
 367 }
 368
 369