1 package fr.orsay.lri.varna.models.treealign;
2 import java.util.ArrayList;
9 * @author Raphael Champeimont
11 public class ExampleDistance3 implements TreeAlignLabelDistanceSymmetric<RNANodeValue2> {
12 public double f(RNANodeValue2 v1, RNANodeValue2 v2) {
16 } else if (!v2.isSingleNode()) { // v2 is a list of bases
17 // We insert all bases, with a cost of 1 for each base.
18 return v2.getNodes().size();
19 } else { // v2 is a single node
22 } else if (!v1.isSingleNode()) { // v1 is a list of bases
24 return v1.getNodes().size();
25 } else if (!v2.isSingleNode()) { // v2 is a list of bases
26 // We compute the sequence distance
27 return alignSequenceNodes(v1, v2).getDistance();
28 } else { // v2 is a single node
29 return 2 + v1.getNodes().size();
31 } else { // v1 is a single node
32 // all the same as when v1 == null
35 } else if (!v2.isSingleNode()) { // v2 is a list of bases
36 return 2 + v2.getNodes().size();
37 } else { // v2 is a single node
38 String l1 = v1.getNode().getLeftNucleotide();
39 String r1 = v1.getNode().getRightNucleotide();
40 String l2 = v2.getNode().getLeftNucleotide();
41 String r2 = v2.getNode().getRightNucleotide();
42 // We have cost(subst((x,y) to (x',y'))) = 1
43 // when x != x' and y != y'.
44 // It means it is less than substituting 2 non-paired bases
45 return (!l1.equals(l2) ? 0.5 : 0)
46 + (!r1.equals(r2) ? 0.5 : 0);
52 public class SequenceAlignResult {
53 private double distance;
54 private int[][] alignment;
56 public double getDistance() {
59 public void setDistance(double distance) {
60 this.distance = distance;
63 /** The result array is a matrix of height 2
64 * and width at most length(sequence A) + length(sequence B).
65 * with result[0] is the alignment for A
66 * and result[1] the alignment for B.
67 * The alignment consists int the indexes of the original
68 * bases positions, with -1 when there is no match.
70 public int[][] getAlignment() {
73 public void setAlignment(int[][] alignment) {
74 this.alignment = alignment;
80 * Align two sequences contained in nodes.
81 * Both nodes have to be non-single nodes, otherwise an
82 * RNANodeValue2WrongTypeException exception will be thrown.
84 public SequenceAlignResult alignSequenceNodes(RNANodeValue2 v1, RNANodeValue2 v2) {
85 char[] A = v1.computeSequence();
86 char[] B = v2.computeSequence();
87 return alignSequences(A, B);
91 * Align sequences using the Needleman-Wunsch algorithm.
92 * Time: O(A.length * B.length)
93 * Space: O(A.length * B.length)
94 * Space used by the returned object: O(A.length + B.length)
96 public SequenceAlignResult alignSequences(char[] A, char[] B) {
97 SequenceAlignResult result = new SequenceAlignResult();
99 final int la = A.length;
100 final int lb = B.length;
101 double[][] F = new double[la+1][lb+1];
102 int[][] decisions = new int[la+1][lb+1];
103 final double d = 1; // insertion/deletion cost
104 final double substCost = 1; // substitution cost
105 for (int i=0; i<=la; i++)
107 for (int j=0; j<=lb; j++)
109 for (int i=1; i<=la; i++)
110 for (int j=1; j<=lb; j++)
114 double match = F[i-1][j-1] + (A[i-1] == B[j-1] ? 0 : substCost);
115 double delete = F[i-1][j] + d;
116 if (match < delete) {
123 double insert = F[i][j-1] + d;
129 decisions[i][j] = decision;
132 result.setDistance(F[la][lb]);
134 int[][] alignment = computeAlignment(F, decisions, A, B);
135 result.setAlignment(alignment);
140 private int[][] computeAlignment(double[][] F, int[][] decisions, char[] A, char[] B) {
141 // At worst the alignment will be of length (A.length + B.length)
142 ArrayList<Integer> AlignmentA = new ArrayList<Integer>(A.length + B.length);
143 ArrayList<Integer> AlignmentB = new ArrayList<Integer>(A.length + B.length);
146 while (i > 0 && j > 0)
148 int decision = decisions[i][j];
167 throw (new Error("Bug in ExampleDistance3: decision = " + decision));
183 // Convert the ArrayLists to the right format:
184 // We need to reverse the list and change the numbering of
185 int l = AlignmentA.size();
186 int[][] result = new int[2][l];
187 for (i=0; i<l; i++) {
188 result[0][i] = AlignmentA.get(l-1-i);
189 result[1][i] = AlignmentB.get(l-1-i);