Mahalanobis distance as measure of dissimilarity between strings (sequences)

Question

I'm doing some research about methods for distance-based comparison of composition of biological sequences (genes, proteins).

Suppose I have two strings (named X and Y) of different lengths, but from a finite alphabet (A, C, T, G):

X = 'ACGT'
Y = 'ACGTA'

The difference between two strings can be quantified by calculating distance between their transition matrices. To do so, we can calculate how many times each letter from the alphabet is present in each string. We obtain two vectors representing letter counts for the sequences:

x = [1,1,1,1] 
y = [2,1,1,1]

Then I can calculate Euclidean distance:

d(x,y) = [(1-2)^2 + (1-1)^2 + (1-1)^2 + (1-1)^2]^0.5 = 1^0.5 = 1

I can't figure out how to calculate the mahalanobis distance. I would be grateful if someone could employ my example and show me how to calculate the mahalanobis distance.

Answer 1

Manual calculation of Mahalanobis Distance is simple but unfortunately a bit lengthy:

>>> # here's the formula i'll use to calculate M/D:
>>> md = (x - y) * LA.inv(R) * (x - y).T

In other words, Mahalanobis distance is the difference (of the 2 data vecctors) multiplied by the inverse of the covariance matrix multiplied by the transpose of the difference (of the same 2 vectors, x & y)

>>> # your 2 data points whose Mahalanobis distance you wish to calculate
>>> x = NP.mat("1 1 1 1")
>>> y = NP.mat("2 1 1 1")

>>> # not enough data supplied in the OP to properly calculate covariance matrix,
>>> # so we'll make some up--a 10 rows of data points of same dimension as x & y
>>> #partition your data into classes (e.g., if you have two classes,
>>> # put all class I data points in one array & all class II points in another)

>>> # for instance pretend 'a' below is the matrix of of your data points
>>> (like x & y) all assigned to the same class
>>> a = NP.random.randint(0, 5, 40).reshape(10, 4)
>>> a
  array([[1, 2, 2, 1],
         [3, 0, 4, 4],
         [2, 3, 1, 1],
         [1, 0, 3, 0],
         [4, 4, 3, 2],
         [4, 0, 0, 4],
         [4, 4, 0, 1],
         [4, 1, 2, 1],
         [4, 0, 3, 4],
         [2, 2, 4, 1]])

>>> # "mean center" this data prior to calculating covariance matrix
>>> mx = NP.mean(a, axis=0)
>>> a1 = a - mx

>>> # sanity check:
>>> NP.mean(a1, axis=0)
  array([ 0., -0., -0.,  0.])

>>> # calculate coveriance matrix of the mean-centered data matrix, a1
>>> R = NP.corrcoef(a1, rowvar=0)
>>> R
  array([[ 1.   ,  0.084, -0.281,  0.561],
         [ 0.084,  1.   , -0.284, -0.461],
         [-0.281, -0.284,  1.   ,  0.059],
         [ 0.561, -0.461,  0.059,  1.   ]])

>>> # quick sanity check(s): 
>>> # (i) is cov matrix n x n? and a; and
>>> # (ii) main diagonal consists of all '1's 
>>> # (because a number and itself of course have perfect covariance)

>>> # repeat those 2 steps (mean center + calculate covariance matrix)
>>> # for the other data matrices (comprised of data points 
>>> # in the remaining classes).

>>> # next calculate 'pooled covariance matrix' by taking weighted average 
>>> of these covariance marices (weighted according to number of rows in 
>>> # the original data matrices used to calculate the covariance matrices

>>> # convert element-wise NumPy arrays to linear algebra matrices
>>> R = NP.matrix(R)    

>>> # calculate the inverse of the weighted average covariance matrix
>>> RI = LA.inv(R)

>>> # now just plug the values into the Mahalanobis code i recited near the top
>>> # we'll do it step-wise so we can see intermediate results:
>>> # another sanity check: we are calculating a distance obviously so the final
>>> # should be a 1 x 1 matrix (scalar)

>>> xy_diff = x - y
>>> a = xy_diff * RI
>>> a
 matrix([[-2.034,  0.737, -0.452,  1.508]])

>>> b = xy_diff.T
>>> a * b
  matrix([[2.043]])     # the Mahalanobis distance for the 2 vectors, x & y

Other (faster) ways to calculate Mahalanobis distance:

The excellent matrix computation mega-library for Python, SciPy, actually has a module "spatial" which inclues a good Mahalanobis function. I can recommend this highly (both the library and the function); I have used this function many times and on several ocassions i cross-verified the results with those from other libraries.

Or you can use R, which has a bult-in function of the same name to calculate M/D, mahalanobis. A concise and useful help page for this function can be accessed by typing in the R interpreter:

?mahalanobis

Finally, i am quite sure that other formulations of Mahalanobis Distance can be found in various R libraries, particularly in some of the libraries in the Bioconductor Project which contains a huge set of R libraries, or "Packages", for the quantitative study of life sciences) then you can calculate Mahalanobis distance using a built-in function of the same name ("mahalanobis.") The reason i mention this is that these domain-specific formulations are likely to have helper functions to save time on the tedious predicate steps e.g., mean-centering and calculating the weighted average covariance matrix.

Answer 2

The Mahalanobis distance is only defined between statistical distributions of (the same number of) numerical values. What you have is a list of bases. You therefore first need to cast your problem into the framework of the Mahalanobis distance.

This can be done as follows. However, I'm not sure whether this procedure has any useful interpretation.

One simple possibility would be to consider that each base (ATGC) has a certain probability of appearing, in each segment. You can estimate the "central" value p of the probability distribution for each base by simply dividing the number of occurrences of each base by the total number of bases in your string. The variance of the probability distribution of each base can be estimated through the binomial law, which applies here: variance = p*(p-1)/N, where N is the number of bases in your string.

You need to calculate probability distributions for the bases of each of your two strings x and y. In order to calculate the covariance matrix S required by the Mahalanobis distance, you can first assume that the distributions of bases in both strings are independent. This gives you a diagonal S matrix, whose diagonal elements are simply covariance(B, B) = var(B_x) + var(B_y), where base B is A, T, G or C, and where var(B_x) is the variance of base B in string x (and the same for y).

This will give you a Mahalanobis distance between two statistical distributions of variables (the probability of finding base B for each ATGC base). However, the interpretation of this distance (what it really represents) is not obvious; again, I am not sure that it has any useful meaning. One reason is that the probability distributions that you define (through the binomial distribution) represent your knowledge of the exact probability of finding a given base in each string: this probability distribution does not represent the real probabilities of finding such or such base in a string.

Anyway, you will get a number that looks like a distance, with the above procedure. This number is zero if the two strings contain the same number of each base, which is at least consistent with the idea of distance.