--- /dev/null
+package jalview.analysis.scoremodels;
+
+import jalview.api.analysis.PairwiseScoreModelI;
+import jalview.api.analysis.SimilarityParamsI;
+import jalview.datamodel.AlignmentView;
+import jalview.math.Matrix;
+import jalview.math.MatrixI;
+import jalview.util.Comparison;
+
+/**
+ * A class to provide sequence pairwise similarity based on residue identity.
+ * Instances of this class are immutable and thread-safe.
+ */
+public class PIDModel extends SimilarityScoreModel implements
+ PairwiseScoreModelI
+{
+ private static final String NAME = "PID";
+
+ /**
+ * Constructor
+ */
+ public PIDModel()
+ {
+ }
+
+ @Override
+ public String getName()
+ {
+ return NAME;
+ }
+
+ /**
+ * Answers null for description. If a display name is needed, use getName() or
+ * an internationalized string built from the name.
+ */
+ @Override
+ public String getDescription()
+ {
+ return null;
+ }
+
+ @Override
+ public boolean isDNA()
+ {
+ return true;
+ }
+
+ @Override
+ public boolean isProtein()
+ {
+ return true;
+ }
+
+ /**
+ * Answers 1 if c and d are the same residue (ignoring case), and not gap
+ * characters. Answers 0 for non-matching or gap characters.
+ */
+ @Override
+ public float getPairwiseScore(char c, char d)
+ {
+ c = toUpper(c);
+ d = toUpper(d);
+ if (c == d && !Comparison.isGap(c))
+ {
+ return 1f;
+ }
+ return 0f;
+ }
+
+ /**
+ * @param c
+ */
+ protected static char toUpper(char c)
+ {
+ if ('a' <= c && c <= 'z')
+ {
+ c += 'A' - 'a';
+ }
+ return c;
+ }
+
+ /**
+ * Computes similarity scores based on pairwise percentage identity of
+ * sequences. For consistency with Jalview 2.10.1's SeqSpace mode PCA
+ * calculation, the percentage scores are rescaled to the width of the
+ * sequences (as if counts of identical residues). This method is thread-safe.
+ */
+ @Override
+ public MatrixI findSimilarities(AlignmentView seqData,
+ SimilarityParamsI options)
+ {
+ String[] seqs = seqData.getSequenceStrings(Comparison.GAP_DASH);
+
+ MatrixI result = findSimilarities(seqs, options);
+
+ result.multiply(seqData.getWidth() / 100d);
+
+ return result;
+ }
+
+ /**
+ * A distance score is computed in the usual way (by reversing the range of
+ * the similarity score results), and then rescaled to percentage values
+ * (reversing the rescaling to count values done in findSimilarities). This
+ * method is thread-safe.
+ */
+ @Override
+ public MatrixI findDistances(AlignmentView seqData,
+ SimilarityParamsI options)
+ {
+ MatrixI result = super.findDistances(seqData, options);
+
+ if (seqData.getWidth() != 0)
+ {
+ result.multiply(100d / seqData.getWidth());
+ }
+
+ return result;
+ }
+
+ /**
+ * Compute percentage identity scores, using the gap treatment and
+ * normalisation specified by the options parameter
+ *
+ * @param seqs
+ * @param options
+ * @return
+ */
+ protected MatrixI findSimilarities(String[] seqs,
+ SimilarityParamsI options)
+ {
+ // TODO reuse code in ScoreMatrix instead somehow
+ double[][] values = new double[seqs.length][];
+ for (int row = 0; row < seqs.length; row++)
+ {
+ values[row] = new double[seqs.length];
+ for (int col = 0; col < seqs.length; col++)
+ {
+ double total = computePID(seqs[row], seqs[col], options);
+ values[row][col] = total;
+ }
+ }
+ return new Matrix(values);
+ }
+
+ /**
+ * Computes a percentage identity for two sequences, using the algorithm
+ * choices specified by the options parameter
+ *
+ * @param seq1
+ * @param seq2
+ * @param options
+ * @return
+ */
+ public static double computePID(String seq1, String seq2,
+ SimilarityParamsI options)
+ {
+ int len1 = seq1.length();
+ int len2 = seq2.length();
+ int width = Math.max(len1, len2);
+ int total = 0;
+ int divideBy = 0;
+
+ for (int i = 0; i < width; i++)
+ {
+ if (i >= len1 || i >= len2)
+ {
+ /*
+ * off the end of one sequence; stop if we are only matching
+ * on the shorter sequence length, else treat as trailing gap
+ */
+ if (options.denominateByShortestLength())
+ {
+ break;
+ }
+ if (options.includeGaps())
+ {
+ divideBy++;
+ }
+ if (options.matchGaps())
+ {
+ total++;
+ }
+ continue;
+ }
+ char c1 = seq1.charAt(i);
+ char c2 = seq2.charAt(i);
+ boolean gap1 = Comparison.isGap(c1);
+ boolean gap2 = Comparison.isGap(c2);
+
+ if (gap1 && gap2)
+ {
+ /*
+ * gap-gap: include if options say so, if so
+ * have to score as identity; else ignore
+ */
+ if (options.includeGappedColumns())
+ {
+ divideBy++;
+ total++;
+ }
+ continue;
+ }
+
+ if (gap1 || gap2)
+ {
+ /*
+ * gap-residue: include if options say so,
+ * count as match if options say so
+ */
+ if (options.includeGaps())
+ {
+ divideBy++;
+ }
+ if (options.matchGaps())
+ {
+ total++;
+ }
+ continue;
+ }
+
+ /*
+ * remaining case is gap-residue
+ */
+ if (toUpper(c1) == toUpper(c2))
+ {
+ total++;
+ }
+ divideBy++;
+ }
+
+ return divideBy == 0 ? 0D : 100D * total / divideBy;
+ }
+}
git://source.jalview.org / jalview.git / blobdiff