+
+ /**
+ * calculate the mean score of the alignment
+ * mean score is equal to the score of an alignmenet of two sequences with randomly shuffled AA sequence composited of the same AA as the two original sequences
+ *
+ */
+ public void meanScore()
+ {
+ //int length = (indelfreeAstr1.length() > indelfreeAstr2.length()) ? indelfreeAstr1.length() : indelfreeAstr2.length();
+ int length = indelfreeAstr1.length(); //both have the same length
+ //create HashMap for counting residues in each sequence
+ HashMap<Character, Integer> seq1ResCount = new HashMap<Character, Integer>();
+ HashMap<Character, Integer> seq2ResCount = new HashMap<Character, Integer>();
+
+ // for both sequences (String indelfreeAstr1 or 2) create a key for the residue and add 1 each time its encountered
+ for (char residue: indelfreeAstr1.toCharArray())
+ {
+ seq1ResCount.putIfAbsent(residue, 0);
+ seq1ResCount.replace(residue, seq1ResCount.get(residue) + 1);
+ }
+ for (char residue: indelfreeAstr2.toCharArray())
+ {
+ seq2ResCount.putIfAbsent(residue, 0);
+ seq2ResCount.replace(residue, seq2ResCount.get(residue) + 1);
+ }
+
+ // meanscore = for each residue pair get the number of appearance and add (countA * countB * pairwiseScore(AB))
+ // divide the meanscore by the sequence length afterwards
+ float _meanscore = 0;
+ for (char resA : seq1ResCount.keySet())
+ {
+ for (char resB : seq2ResCount.keySet())
+ {
+ int countA = seq1ResCount.get(resA);
+ int countB = seq2ResCount.get(resB);
+
+ float scoreAB = scoreMatrix.getPairwiseScore(resA, resB);
+
+ _meanscore += countA * countB * scoreAB;
+ }
+ }
+ _meanscore /= length;
+ this.meanScore = _meanscore;
+ }
+
+ public float getMeanScore()
+ {
+ return this.meanScore;
+ }
+
+ /**
+ * calculate the hypothetic max score using the self-alignment of the sequences
+ */
+ public void hypotheticMaxScore()
+ {
+ int _hmsA = 0;
+ int _hmsB = 0;
+ for (char residue: indelfreeAstr1.toCharArray())
+ {
+ _hmsA += scoreMatrix.getPairwiseScore(residue, residue);
+ }
+ for (char residue: indelfreeAstr2.toCharArray())
+ {
+ _hmsB += scoreMatrix.getPairwiseScore(residue, residue);
+ }
+ this.hypotheticMaxScore = (_hmsA < _hmsB) ? _hmsA : _hmsB; // take the lower self alignment
+
+ }
+
+ public int getHypotheticMaxScore()
+ {
+ return this.hypotheticMaxScore;
+ }
+
+ /**
+ * create strings based of astr1 and astr2 but without gaps
+ */
+ public void getIndelfreeAstr()
+ {
+ int n = astr1.length(); // both have the same length
+ for (int i = 0; i < n; i++)
+ {
+ if (Character.isLetter(astr1.charAt(i)) && Character.isLetter(astr2.charAt(i))) // if both sequences dont have a gap -> add to indelfreeAstr
+ {
+ this.indelfreeAstr1 += astr1.charAt(i);
+ this.indelfreeAstr2 += astr2.charAt(i);
+ }
+ }
+ }
+
+ /**
+ * calculates the overall score of the alignment
+ * preprescore = sum of all scores - all penalties
+ * if preprescore < 1 ~ alignmentScore = Float.NaN >
+ * alignmentScore = ((preprescore - meanScore) / (hypotheticMaxScore - meanScore)) * coverage
+ */
+ public void scoreAlignment() throws RuntimeException
+ {
+
+ getIndelfreeAstr();
+ meanScore();
+ hypotheticMaxScore();
+ // cannot calculate score because denominator would be zero
+ if (this.hypotheticMaxScore == this.meanScore)
+ {
+ throw new IllegalArgumentException(String.format("hypotheticMaxScore (%8.2f) == meanScore (%8.2f) - division by 0", hypotheticMaxScore, meanScore));
+ }
+ //int n = (astr1.length() > astr2.length()) ? astr1.length() : astr2.length();
+ int n = indelfreeAstr1.length();
+
+ float score = 0;
+ boolean aGapOpen = false;
+ boolean bGapOpen = false;
+ for (int i = 0; i < n; i++)
+ {
+ char char1 = indelfreeAstr1.charAt(i);
+ char char2 = indelfreeAstr2.charAt(i);
+ boolean aIsLetter = Character.isLetter(char1);
+ boolean bIsLetter = Character.isLetter(char2);
+ if (aIsLetter && bIsLetter) // if pair -> get score
+ {
+ score += scoreMatrix.getPairwiseScore(char1, char2);
+ } else if (!aIsLetter && !bIsLetter) { // both are gap -> skip
+ } else if ((!aIsLetter && aGapOpen) || (!bIsLetter && bGapOpen)) { // one side gapopen -> score - gap_extend
+ score -= GAP_EXTEND_COST;
+ } else { // no gap open -> score - gap_open
+ score -= GAP_OPEN_COST;
+ }
+ // adjust GapOpen status in both sequences
+ aGapOpen = (!aIsLetter) ? true : false;
+ bGapOpen = (!bIsLetter) ? true : false;
+ }
+
+ float preprescore = score; // if this score < 1 --> alignment score = Float.NaN
+ score = (score - this.meanScore) / (this.hypotheticMaxScore - this.meanScore);
+ int[] _max = MiscMath.findMax(new int[]{astr1.replace("-","").length(), astr2.replace("-","").length()}); // {index of max, max}
+ float coverage = (float) n / (float) _max[1]; // indelfreeAstr length / longest sequence length
+ float prescore = score; // only debug
+ score *= coverage;
+
+ System.out.println(String.format("prepre-score: %f, pre-score: %f, longlength: %d\nscore: %f, mean: %f, max: %d", preprescore, prescore, _max[1], score, this.meanScore, this.hypotheticMaxScore));
+ this.alignmentScore = (preprescore < 1) ? Float.NaN : score;
+ }