The longest common subsequence and sequence alignment problems have been studied extensively and they can be regarded as the relationship measurement between sequences. However, most of them treat sequences evenly or consider only two sequences. Recently, with the rise of whole-genome duplication research, the doubly conserved synteny relationship among three sequences should be considered. It is a brand new model to find a merging way for understanding the interleaving relationship of sequences. Here, we define the merged LCS problem for measuring the interleaving relationship among three sequences. An O (n3 ) algorithm is first proposed for solving the problem, where n is the sequence length. We further discuss the variant version of this problem with the block information. For the blocked merged LCS problem, we propose an algorithm with time complexity O (n2 m), where m is the number of blocks. An improved O (n2 + nm2 ) algorithm is further proposed for the same blocked problem.
int BMLCS_DP(vector< int> T, vector< int> A, vector< int> B, int Tlen, int Alen, int Blen, vector< int> EOBA, vector< int> EOBB) {
vector< vector< int>> vat(Alen + 1, vector< int>(Tlen + 1));
vector< vector< int>> vbt(Blen + 1, vector< int>(Tlen + 1));
vector< vector< int>> vabt1(Alen + 1, vector< int>(Blen + 1));
vector< vector< int>> vabt2(Alen + 1, vector< int>(Blen + 1));
//use LCS algorithm to calcluate similiar between sequence A and sequence T
vat[0][0] = 0;
for (int i = 1; i <= Alen; ++i) {
vat[i][0] = 0;
}
for (int j = 1; j <= Tlen; ++j) {
vat[0][j] = 0;
}
for (int i = 1; i <= Alen; ++i) {
for (int j = 1; j <= Tlen; ++j) {
if (A[i] != T[j]) {
vat[i][j] = max(vat[i - 1][j], vat[i][j - 1]);
}
else {
vat[i][j] = vat[i - 1][j - 1] + 1;
}
}
}
//use LCS algorithm to calcluate similiar between sequence B and sequence T
vbt[0][0] = 0;
for (int i = 1; i <= Blen; ++i) {
vbt[i][0] = 0;
}
for (int j = 1; j <= Tlen; ++j) {
vbt[0][j] = 0;
}
for (int i = 1; i <= Blen; ++i) {
for (int j = 1; j <= Tlen; ++j) {
if (B[i] != T[j]) {
vbt[i][j] = max(vbt[i - 1][j], vbt[i][j - 1]);
}
else {
vbt[i][j] = vbt[i - 1][j - 1] + 1;
}
}
}
for (int i = 0; i <= Alen; ++i) {
for (int j = 0; j <= Blen; ++j) {
vabt1[i][j] = 0;
}
}
for (int i = 0; i <= Alen; ++i) {
for (int j = 0; j <= Blen; ++j) {
vabt2[i][j] = 0;
}
}
//use BMLCS-DP algorithm
int max = 0;
int cali = 0, calj = 0;
int tmpb = 0;
int layer = 1;
while (1) {
vabt1[0][0] = 0;
vabt2[0][0] = 0;
if (layer % 2 == 1) {
cali = 0, calj = 0;
for (int i = 1; i <= Alen; ++i) {
vabt1[i][0] = vat[i][layer];
if (i < EOBA[cali]) { // ai!=EOB amd bj = EOB
for (int j = 0; j < EOBB.size(); ++j) {
max = 0;
tmpb = EOBB[j];
if (A[i] == T[layer]) {
if ((vabt2[i - 1][tmpb] + 1) > max)
max = vabt2[i - 1][tmpb] + 1;
}
else {
if (vabt1[i - 1][tmpb] > max)
max = vabt1[i - 1][tmpb];
if (vabt2[i][tmpb] > max)
max = vabt2[i][tmpb];
}
vabt1[i][tmpb] = max;
}
}
else {// i = EOBA[cali]
calj = 0;
for (int j = 1; j <= Blen; ++j) {
vabt1[0][j] = vbt[j][layer];
max = 0;
if (j < EOBB[calj]) {// ai=EOB amd bj!= EOB
if (B[j] == T[layer]) {
if ((vabt2[i][j - 1] + 1) > max)
max = vabt2[i][j - 1] + 1;
}
else {
if (vabt1[i][j - 1] > max)
max = vabt1[i][j - 1];
if (vabt2[i][j] > max)
max = max = vabt2[i][j];
}
vabt1[i][j] = max;
}
else { // ai=EOB amd bj= EOB
if (A[i] == T[layer]) {
if ((vabt2[i - 1][j] + 1) > max)
max = vabt2[i - 1][j] + 1;
}
else {
if (vabt1[i - 1][j] > max)
max = vabt1[i - 1][j];
}
if (B[j] == T[layer]) {
if ((vabt2[i][j - 1] + 1) > max)
max = vabt2[i][j - 1] + 1;
}
else {
if (vabt1[i][j - 1] > max)
max = vabt1[i][j - 1];
}
if (A[i] != T[layer] || B[j] != T[layer]) {
if (vabt2[i][j] > max)
max = vabt2[i][j];
}
vabt1[i][j] = max;
++calj;
}
}
++cali;
}
}
if (layer == Tlen)
return vabt1[Alen][Blen];
++layer;
}
else {
cali = 0, calj = 0;
for (int i = 1; i <= Alen; ++i) {
vabt2[i][0] = vat[i][layer];
if (i < EOBA[cali]) { // ai!=EOB amd bj = EOB
for (int j = 0; j < EOBB.size(); ++j) {
max = 0;
tmpb = EOBB[j];
if (A[i] == T[layer]) {
if ((vabt1[i - 1][tmpb] + 1) > max)
max = vabt1[i - 1][tmpb] + 1;
}
else {
if (vabt2[i - 1][tmpb] > max)
max = vabt2[i - 1][tmpb];
if (vabt1[i][tmpb] > max)
max = vabt1[i][tmpb];
}
vabt2[i][tmpb] = max;
}
}
else {// i = EOBA[cali]
calj = 0;
for (int j = 1; j <= Blen; ++j) {
vabt2[0][j] = vbt[j][layer];
max = 0;
if (j < EOBB[calj]) {// ai=EOB amd bj!= EOB
if (B[j] == T[layer]) {
if ((vabt1[i][j - 1] + 1) > max)
max = vabt1[i][j - 1] + 1;
}
else {
if (vabt2[i][j - 1] > max)
max = vabt2[i][j - 1];
if (vabt1[i][j] > max)
max = max = vabt1[i][j];
}
vabt2[i][j] = max;
}
else { // ai=EOB amd bj= EOB
if (A[i] == T[layer]) {
if ((vabt1[i - 1][j] + 1) > max)
max = vabt1[i - 1][j] + 1;
}
else {
if (vabt2[i - 1][j] > max)
max = vabt2[i - 1][j];
}
if (B[j] == T[layer]) {
if ((vabt1[i][j - 1] + 1) > max)
max = vabt1[i][j - 1] + 1;
}
else {
if (vabt2[i][j - 1] > max)
max = vabt2[i][j - 1];
}
if (A[i] != T[layer] || B[j] != T[layer]) {
if (vabt1[i][j] > max)
max = vabt1[i][j];
}
vabt2[i][j] = max;
++calj;
}
}
++cali;
}
}
if (layer == Tlen)
return vabt2[Alen][Blen];
++layer;
}
}
}