sequence.cc

Libgist is an implementation of the Generalized Search Tree, a template index structure that makes i

// sequence.cc					-*- c++ -*-
// Copyright (c) 1997, 1998, Regents of the University of California
// $Id: sequence.cc,v 1.3 2000/03/15 00:24:40 mashah Exp $

#include <stdio.h>

// VCPORT_B
#ifdef WIN32
#include <iostream>

using namespace std;
#else
#include <iostream.h>

#endif
// VCPORT_E

#include "sequence.h"
#include <string.h> // for strlen

#define max(a, b)  ((a) < (b) ? (b) : (a))


ScoringMatrix::Matrix blosum62Matrix = {
    {4, -1, -2, -2, 0, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -3, -2, 0},
    {-1, 5, 0, -2, -3, 1, 0, -2, 0, -3, -2, 2, -1, -3, -2, -1, -1, -3, -2, -3},
    {-2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3, 0, -2, -3, -2, 1, 0, -4, -2, -3},
    {-2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4, -3, -3},
    {0, -3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1},
    {-1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2, -1, -2},
    {-1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3, 1, -2, -3, -1, 0, -1, -3, -2, -2},
    {0, -2, 0, -1, -3, -2, -2, 6, -2, -4, -4, -2, -3, -3, -2, 0, -2, -2, -3, -3},
    {-2, 0, 1, -1, -3, 0, 0, -2, 8, -3, -3, -1, -2, -1, -2, -1, -2, -2, 2, -3},
    {-1, -3, -3, -3, -1, -3, -3, -4, -3, 4, 2, -3, 1, 0, -3, -2, -1, -3, -1, 3},
    {-1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, -2, 2, 0, -3, -2, -1, -2, -1, 1},
    {-1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, -1, -3, -1, 0, -1, -3, -2, -2},
    {-1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, 0, -2, -1, -1, -1, -1, 1},
    {-2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, -4, -2, -2, 1, 3, -1},
    {-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, -1, -1, -4, -3, -2},
    {1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, 1, -3, -2, -2},
    {0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, -2, -2, 0},
    {-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, 2, -3},
    {-2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7, -1},
    {0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4}
};
ScoringMatrix blosum62(blosum62Matrix);

ScoringMatrix::Matrix blosum45Matrix = {
    {5, -2, -1, -2, -1, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -2, -2, 0},
    {-2, 7, 0, -1, -3, 1, 0, -2, 0, -3, -2, 3, -1, -2, -2, -1, -1, -2, -1, -2},
    {-1, 0, 6, 2, -2, 0, 0, 0, 1, -2, -3, 0, -2, -2, -2, 1, 0, -4, -2, -3},
    {-2, -1, 2, 7, -3, 0, 2, -1, 0, -4, -3, 0, -3, -4, -1, 0, -1, -4, -2, -3},
    {-1, -3, -2, -3, 12, -3, -3, -3, -3, -3, -2, -3, -2, -2, -4, -1, -1, -5, -3, -1},
    {-1, 1, 0, 0, -3, 6, 2, -2, 1, -2, -2, 1, 0, -4, -1, 0, -1, -2, -1, -3},
    {-1, 0, 0, 2, -3, 2, 6, -2, 0, -3, -2, 1, -2, -3, 0, 0, -1, -3, -2, -3},
    {0, -2, 0, -1, -3, -2, -2, 7, -2, -4, -3, -2, -2, -3, -2, 0, -2, -2, -3, -3},
    {-2, 0, 1, 0, -3, 1, 0, -2, 10, -3, -2, -1, 0, -2, -2, -1, -2, -3, 2, -3},
    {-1, -3, -2, -4, -3, -2, -3, -4, -3, 5, 2, -3, 2, 0, -2, -2, -1, -2, 0, 3},
    {-1, -2, -3, -3, -2, -2, -2, -3, -2, 2, 5, -3, 2, 1, -3, -3, -1, -2, 0, 1},
    {-1, 3, 0, 0, -3, 1, 1, -2, -1, -3, -3, 5, -1, -3, -1, -1, -1, -2, -1, -2},
    {-1, -1, -2, -3, -2, 0, -2, -2, 0, 2, 2, -1, 6, 0, -2, -2, -1, -2, 0, 1},
    {-2, -2, -2, -4, -2, -4, -3, -3, -2, 0, 1, -3, 0, 8, -3, -2, -1, 1, 3, 0},
    {-1, -2, -2, -1, -4, -1, 0, -2, -2, -2, -3, -1, -2, -3, 9, -1, -1, -3, -3, -3},
    {1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -3, -1, -2, -2, -1, 4, 2, -4, -2, -1},
    {0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -1, -1, 2, 5, -3, -1, 0},
    {-2, -2, -4, -4, -5, -2, -3, -2, -3, -2, -2, -2, -2, 1, -3, -4, -3, 15, 3, -3},
    {-2, -1, -2, -2, -3, -1, -2, -3, 2, 0, 0, -1, 0, 3, -3, -2, -1, 3, 8, -1},
    {0, -2, -3, -3, -1, -3, -3, -3, -3, 3, 1, -2, 1, 0, -3, -1, 0, -3, -1, 5}
};
ScoringMatrix blosum45(blosum45Matrix);

ScoringMatrix::Matrix blosum80Matrix = {
    {5, -2, -2, -2, -1, -1, -1, 0, -2, -2, -2, -1, -1, -3, -1, 1, 0, -3, -2, 0},
    {-2, 6, -1, -2, -4, 1, -1, -3, 0, -3, -3, 2, -2, -4, -2, -1, -1, -4, -3, -3},
    {-2, -1, 6, 1, -3, 0, -1, -1, 0, -4, -4, 0, -3, -4, -3, 0, 0, -4, -3, -4},
    {-2, -2, 1, 6, -4, -1, 1, -2, -2, -4, -5, -1, -4, -4, -2, -1, -1, -6, -4, -4},
    {-1, -4, -3, -4, 9, -4, -5, -4, -4, -2, -2, -4, -2, -3, -4, -2, -1, -3, -3, -1},
    {-1, 1, 0, -1, -4, 6, 2, -2, 1, -3, -3, 1, 0, -4, -2, 0, -1, -3, -2, -3},
    {-1, -1, -1, 1, -5, 2, 6, -3, 0, -4, -4, 1, -2, -4, -2, 0, -1, -4, -3, -3},
    {0, -3, -1, -2, -4, -2, -3, 6, -3, -5, -4, -2, -4, -4, -3, -1, -2, -4, -4, -4},
    {-2, 0, 0, -2, -4, 1, 0, -3, 8, -4, -3, -1, -2, -2, -3, -1, -2, -3, 2, -4},
    {-2, -3, -4, -4, -2, -3, -4, -5, -4, 5, 1, -3, 1, -1, -4, -3, -1, -3, -2, 3},
    {-2, -3, -4, -5, -2, -3, -4, -4, -3, 1, 4, -3, 2, 0, -3, -3, -2, -2, -2, 1},
    {-1, 2, 0, -1, -4, 1, 1, -2, -1, -3, -3, 5, -2, -4, -1, -1, -1, -4, -3, -3},
    {-1, -2, -3, -4, -2, 0, -2, -4, -2, 1, 2, -2, 6, 0, -3, -2, -1, -2, -2, 1},
    {-3, -4, -4, -4, -3, -4, -4, -4, -2, -1, 0, -4, 0, 6, -4, -3, -2, 0, 3, -1},
    {-1, -2, -3, -2, -4, -2, -2, -3, -3, -4, -3, -1, -3, -4, 8, -1, -2, -5, -4, -3},
    {1, -1, 0, -1, -2, 0, 0, -1, -1, -3, -3, -1, -2, -3, -1, 5, 1, -4, -2, -2},
    {0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -2, -1, -1, -2, -2, 1, 5, -4, -2, 0},
    {-3, -4, -4, -6, -3, -3, -4, -4, -3, -3, -2, -4, -2, 0, -5, -4, -4, 11, 2, -3},
    {-2, -3, -3, -4, -3, -2, -3, -4, 2, -2, -2, -3, -2, 3, -4, -2, -2, 2, 7, -2},
    {0, -3, -4, -4, -1, -3, -3, -4, -4, 3, 1, -3, 1, -1, -3, -2, 0, -3, -2, 4}
};
ScoringMatrix blosum80(blosum80Matrix);

bool ScoringMatrix::_initialized = false;
int ScoringMatrix::_index[256];

ScoringMatrix::ScoringMatrix(Matrix m)
{
    int i, j;
    for (i = 0; i < 20; i++) {
        for (j = 0; j < 20; j++) {
	    matrix[i][j] = m[i][j];
	}
    }
}


int
ScoringMatrix::score(char a, char b) const
{
    return(matrix[_index[a]][_index[b]]);
}

void
ScoringMatrix::init()
{
    if (_initialized) {
        return;
    }

    int i;
    for (i = 0; i < 255; i++) {
        _index[i] = -1;
    }

    // assign indices to amino acids
    _index['A'] = 0;
    _index['R'] = 1;
    _index['N'] = 2;
    _index['D'] = 3;
    _index['C'] = 4;
    _index['Q'] = 5;
    _index['E'] = 6;
    _index['G'] = 7;
    _index['H'] = 8;
    _index['I'] = 9;
    _index['L'] = 10;
    _index['K'] = 11;
    _index['M'] = 12;
    _index['F'] = 13;
    _index['P'] = 14;
    _index['S'] = 15;
    _index['T'] = 16;
    _index['W'] = 17;
    _index['Y'] = 18;
    _index['V'] = 19;

    _initialized = true;
}


///////////////////////////////////////////////////////////////////////////////
// Sequence
///////////////////////////////////////////////////////////////////////////////

Sequence::Sequence(const char* s) : seq(s)
{
}

void
Sequence::localAlign(
    const Sequence& s,
    const ScoringMatrix& m,
    int gapPenalty,
    double& score,
    double& prob)
{
    ScoringMatrix::init();

    static int count = 0;
    int* alignMatrix = new int[(strlen(seq) + 1) * (strlen(s.seq) + 1)];
    int slen = strlen(seq) + 1;

    int i, j;
    int maxScore = 0;
    for (i = 0; i <= strlen(seq); i++) {
        for (j = 0; j <= strlen(s.seq); j++) {
	    if (j == 0 || i == 0) {
	        alignMatrix[i * slen + j] = 0;
	    } else {
		int a = alignMatrix[(i-1) * slen + j-1] + m.score(seq[i-1], s.seq[j-1]);
		int b = alignMatrix[(i-1) * slen + j] - gapPenalty;
		int c = alignMatrix[i * slen + j-1] - gapPenalty;
	        alignMatrix[i * slen + j] = max(0, max(a, max(b, c)));
	    }
	    if (alignMatrix[i * slen + j] > maxScore) {
	        maxScore = alignMatrix[i * slen + j];
	    }

	    // debug
	    //cout << alignMatrix[i * slen + j] << " ";
	}

	// debug
	//cout << endl;
    }

    delete [] alignMatrix;
    score = maxScore;

    // debug
    count++;
    cout << count << endl;
}