src/jalview/analysis/scoremodels/PIDModel.java

   1 package jalview.analysis.scoremodels;
   2
   3 import jalview.api.analysis.PairwiseScoreModelI;
   4 import jalview.api.analysis.SimilarityParamsI;
   5 import jalview.api.analysis.SimilarityScoreModelI;
   6 import jalview.datamodel.AlignmentView;
   7 import jalview.math.Matrix;
   8 import jalview.math.MatrixI;
   9 import jalview.util.Comparison;
  10
  11 /**
  12  * A class to provide sequence pairwise similarity based on residue identity
  13  */
  14 public class PIDModel implements SimilarityScoreModelI,
  15         PairwiseScoreModelI
  16 {
  17
  18   @Override
  19   public String getName()
  20   {
  21     return "% Identity (PID)";
  22   }
  23
  24   @Override
  25   public boolean isDNA()
  26   {
  27     return true;
  28   }
  29
  30   @Override
  31   public boolean isProtein()
  32   {
  33     return true;
  34   }
  35
  36   /**
  37    * Answers 1 if c and d are the same residue (ignoring case), and not gap
  38    * characters. Answers 0 for non-matching or gap characters.
  39    */
  40   @Override
  41   public float getPairwiseScore(char c, char d)
  42   {
  43     c = toUpper(c);
  44     d = toUpper(d);
  45     if (c == d && !Comparison.isGap(c))
  46     {
  47       return 1f;
  48     }
  49     return 0f;
  50   }
  51
  52   /**
  53    * @param c
  54    */
  55   protected static char toUpper(char c)
  56   {
  57     if ('a' <= c && c <= 'z')
  58     {
  59       c += 'A' - 'a';
  60     }
  61     return c;
  62   }
  63
  64   @Override
  65   public MatrixI findSimilarities(AlignmentView seqData)
  66   {
  67     // TODO reuse code in ScoreMatrix instead somehow
  68     String[] seqs = seqData.getSequenceStrings(' ');
  69     return findSimilarities(seqs, SimilarityParams.Jalview);
  70   }
  71
  72   /**
  73    * Compute percentage identity scores, using the gap treatment and
  74    * normalisation specified by the options parameter
  75    *
  76    * @param seqs
  77    * @param options
  78    * @return
  79    */
  80   protected MatrixI findSimilarities(String[] seqs,
  81           SimilarityParamsI options)
  82   {
  83     double[][] values = new double[seqs.length][];
  84     for (int row = 0; row < seqs.length; row++)
  85     {
  86       values[row] = new double[seqs.length];
  87       for (int col = 0; col < seqs.length; col++)
  88       {
  89         double total = computePID(seqs[row], seqs[col], options);
  90         values[row][col] = total;
  91       }
  92     }
  93     return new Matrix(values);
  94   }
  95
  96   /**
  97    * Computes a percentage identity for two sequences, using the algorithm
  98    * choices specified by the options parameter
  99    *
 100    * @param seq1
 101    * @param seq2
 102    * @param options
 103    * @return
 104    */
 105   protected double computePID(String seq1, String seq2,
 106           SimilarityParamsI options)
 107   {
 108     int len1 = seq1.length();
 109     int len2 = seq2.length();
 110     int width = Math.max(len1, len2);
 111     int total = 0;
 112     int divideBy = 0;
 113
 114     for (int i = 0; i < width; i++)
 115     {
 116       if (i >= len1 || i >= len2)
 117       {
 118         /*
 119          * off the end of one sequence; stop if we are only matching
 120          * on the shorter sequence length, else treat as trailing gap
 121          */
 122         if (options.denominateByShortestLength())
 123         {
 124           break;
 125         }
 126         if (options.denominatorIncludesGaps())
 127         {
 128           divideBy++;
 129         }
 130         if (options.matchGaps())
 131         {
 132           total++;
 133         }
 134         continue;
 135       }
 136       char c1 = seq1.charAt(i);
 137       char c2 = seq2.charAt(i);
 138       boolean gap1 = Comparison.isGap(c1);
 139       boolean gap2 = Comparison.isGap(c2);
 140
 141       if (gap1 && gap2)
 142       {
 143         /*
 144          * gap-gap: include if options say so, if so
 145          * have to score as identity; else ignore
 146          */
 147         if (options.includeGappedColumns())
 148         {
 149           divideBy++;
 150           total++;
 151         }
 152         continue;
 153       }
 154
 155       if (gap1 || gap2)
 156       {
 157         /*
 158          * gap-residue: include if options say so,
 159          * count as match if options say so
 160          */
 161         if (options.denominatorIncludesGaps())
 162         {
 163           divideBy++;
 164         }
 165         if (options.matchGaps())
 166         {
 167           total++;
 168         }
 169         continue;
 170       }
 171
 172       /*
 173        * remaining case is gap-residue
 174        */
 175       if (toUpper(c1) == toUpper(c2))
 176       {
 177         total++;
 178       }
 179       divideBy++;
 180     }
 181
 182     return divideBy == 0 ? 0D : 100D * total / divideBy;
 183   }
 184 }