Unifrac (8): understanding PCoA---calculate

I'm exploring Unifrac (first and second posts).

In the previous post, I showed code and some output related to our topic for this post. You'll need to generate those numbers to follow along. Of course, you could just paste 'em in:

[[-20 -16 -12 -8 -4 0 4 8 12 16] [ -2 1 2 2 -1 1 -1 -3 -2 5]]

The whole trick to PCoA is that it can take a "measure of association" between two samples, like a distance. So, let's use the values from before, calculate the distances between them, and then see how well PCoA recovers the original data. The second half of the code listing is at the bottom of the post. In the first segment of this second half, we simply define two functions to calculate the distances between all the points, using either the "manhattan" distance or the "euclidean" distance. We will only use the euclidean version for this analysis. If you want to see the results of the distance calculation, you'll have to run the code (manhattan) or modify it to print the euclidean result.

Next, we want to do the PCoA analysis. I did this two ways (to check myself). First, I looked under the hood at the PyCogent module cogent.cluster.metric_scaling, and I call the same functions used in metric_scaling.principal_coordinates_analysis. These include the following two functions as well as a numpy.linalg routine:

E = metric_scaling.make_E_matrix(eucl) F = metric_scaling.make_F_matrix(E) eigvals, eigvecs = eigh(F)

followed by a little massaging of the vectors to generate the output. As you will see if you run the code, the eigenvalues and then the last two rows of my result are:

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.3 36.34 ... [ 1.9 -1. -2. -1.9 1.2 -0.8 1.3 3.4 2.4 -4.5] [-18. -14. -10. -6. -2. 2. 6. 9.9 14. 18.1]

The other approach is to call PyCogent's function directly. The eigenvalues and the last two results of the call to

pca = metric_scaling.principal_coordinates_analysis point_matrix, eigvals = pca(eucl)

are

-0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 53.28 1320.32 ... [ 1.9 -1. -2. -1.9 1.2 -0.8 1.3 3.4 2.4 -4.5] [-18. -14. -10. -6. -2. 2. 6. 9.9 14. 18.1]

This looks pretty good, though I notice that my eigenvalues aren't matching. I'm missing a square root somewhere. I will postpone worrying about this for another day.

>>> 1320.32/53.28 24.780780780780781 >>> (36.34/7.3)**2 24.781302308125362

The original matrix is:

[[-20 -16 -12 -8 -4 0 4 8 12 16] [ -2 1 2 2 -1 1 -1 -3 -2 5]]

I think this is pretty remarkable. We've taken a matrix of distances between our sample points, done a little PCoA mumbo-jumbo, and regenerated something pretty close to the original.

However, I am a bit worried at what look to be bugs in the implementation:

• the original x-values are systematically transformed by + 2
• the original y-values are systematically transformed by * -1

This gives a shift and a reflection of the points, but the relationships revealed by the PCoA are identical to the original.

[UPDATE: <whacks self on head> Duh...I bet that in view of the last line above, if we calculate the distance we'll get the same result as we did for the original matrix. So, of course, there's no reason that PCoA could possibly tell the difference. Presumably, it's just centered its result a little better around (0,0). ]

#============================================ def delta(n1,n2): return abs(n2-n1) def manhattan_dist(x1,y1,x2,y2): return delta(x1,x2) + delta(y1,y2) def euclidean_dist(x1,y1,x2,y2): sumsq = delta(x1,x2)**2 + delta(y1,y2)**2 return np.sqrt(sumsq) eucl = list() manh = list() for i in range(len(x)): for j in range(len(x)): x1,y1 = A[0][i], A[1][i] x2,y2 = A[0][j], A[1][j] manh.append(manhattan_dist(x1,y1,x2,y2)) eucl.append(euclidean_dist(x1,y1,x2,y2)) manh = np.array(manh) manh.shape = (N,N) eucl = np.array(eucl) eucl.shape = (N,N) print 'manhattan\n', manh, '\n' #============================================ # PyCogent stuff starts here: # I call their functions E = metric_scaling.make_E_matrix(eucl) F = metric_scaling.make_F_matrix(E) eigvals, eigvecs = eigh(F) t_eigvecs = eigvecs.transpose() print 'eigvals' for e in eigvals: print round(e,2), print print 't_eigvecs' for row in t_eigvecs: print np.round(row,decimals=1) eigval_root = np.sqrt(abs(eigvals)) print 'eigval_root\n' for e in eigval_root: print round(e,2), print final = t_eigvecs * eigval_root[:,np.newaxis] print 'my calculation' for row in final: print np.round(row,decimals=1) print # They call their functions print 'PyCogent' pca = metric_scaling.principal_coordinates_analysis point_matrix, eigvals = pca(eucl) print 'eigvals' for e in eigvals: print round(e,2), print print 'point_matrix' for row in point_matrix: print np.round(row,decimals=1) print print 'original matrix (A)\n', A

Unifrac (7): understanding PCoA---setup

I'm exploring Unifrac (first and second posts).

One of the nice aspects of UniFrac is the "multivariate analysis" it can do including PCoA (Principal Coordinates Analysis). I'm trying to understand how this works, and I have two posts to go through it. This first post is similar to stuff we've done before (e.g. here), but I want to try to present a coherent picture. Sorry for the repetition. The code is presented in two parts overall, and then today's section has two parts also. We begin by getting a bunch of integers in x-coordinates over the range from -20 to +20, and in y-coordinates randomly between -5 and +5. We turn those into a 2D array and plot them using matplotlib as red circles.

A [[-20 -16 -12 -8 -4 0 4 8 12 16] [ -2 1 2 2 -1 1 -1 -3 -2 5]]

Next, the array is rotated through 45 degrees and replotted (blue circles).

In the second part (of the first half), we use standard Principal Component Analysis (PCA) on these data points to recover their long dimension (now rotated) back to the x-axis. We calculate the covariance matrix for these two sets of points, and then do np.linalg.eigh(m) to get the eigenvalues and eigenvectors.

cov [[ 74.08888889 70.35555556] [ 70.35555556 78.53333333]] eigenvalues [ 5.9204692 146.70175303] eigenvectors [[-0.71818168 0.69585564] [ 0.69585564 0.71818168]]

The eigenvalues come out with the largest one last. We grab the first eigenvector (the last in the array we received), and do

np.dot(ev1,B)

which gives back the x-values we started with.

[-20. -16. -12. -8. -4. 0. 4. 8. 12. 16.]

We can do the same with the second eigenvector to get the y-values. We can combine these two operations by re-ordering the eigenvectors:

print 'C' T2 = np.array([ev1,ev2]) C = np.dot(T2,B) for row in C: print np.round(row) print

C [-20. -16. -12. -8. -4. 0. 4. 8. 12. 16.] [-2. 1. 2. 2. -1. 1. -1. -3. -2. 5.]

In the figure, the C array has been plotted as black cross-bars, which match the original red circles pretty well. Next time, we'll take the same data points and look at them from the perspective of PCoA. We'll also see how PyCogent fares with our data, and how reasonable a picture it can estimate of what we started with.

Here's the code:

import numpy as np from pylab import * from cogent.cluster import metric_scaling plotting = len(sys.argv) > 1 np.random.seed(1357) x = range(-20,20,4) y = np.random.random_integers(-5,5,size=len(x)) A = np.array([x,y]) N = len(x) A.shape=(2,N) print 'A\n', A n = 1.0/np.sqrt(2) T = np.array([[n,-n],[n,n]]) B = np.dot(T,A) if plotting: scatter(A[0],A[1],s=250,color='r') scatter(B[0],B[1],s=250,color='b') #============================================ m = np.cov(B) print 'cov\n', m e = np.linalg.eigh(m) print 'eigenvalues\n', e[0] print 'eigenvectors\n', e[1] ev1 = e[1][1] print np.round(np.dot(ev1,B)) ev2 = e[1][0] print np.round(np.dot(ev2,B)) print print 'C' T2 = np.array([ev1,ev2]) C = np.dot(T2,B) for row in C: print np.round(row) print if plotting: scatter(C[0],C[1],s=100,c='k',marker='+') R = (-25,25) xlim(R) ylim(R) savefig('plot.png')

Unifrac (6): understanding how PyCogent does PCoA

I'm exploring Unifrac (first and second posts).

Continuing with the example from last time, I am trying to understand how the PCoA analysis works in UniFrac. I have made some progress, with a lot of help from Justin Kuczynski. The theoretical basis for what we're going to do is given in a book called "Principles of Multivariate Analysis" by A. Krzanowski (Amazon, Google Books, p 107), (among other sources).

I don't understand it (yet).

But the basic approach is that we're going to do something similar to what is done in Principal Component Analysis (PCA, see here, for example), however, the mathematical treatment is different. In PCA we have data, values for a number of measurements (x,y,z, up to n dimensions) for each of a number of samples or entities. PCA will shows projections of the data onto two-dimensions (planes), where the planes are chosen so as to maximize the spread of the overall variance in that plane. We're reducing multi-dimensional data to fewer dimensions, while preserving as much of the structure as we can.

And in PCA, again, we make a covariation matrix, and then get the eigenvalues and eigenvectors for that matrix.

In the UniFrac flavor of Principal Coordinates Analysis (PCoA) we have a distance measurement on a number of samples (based of course, on the phylogenetic distance as measured by the UniFrac metric). What Krzanowski shows is how to do something similar to PCA but starting with the distance data, or indeed, any "measure of association."

Without worrying too much further about the individual steps, here they are. First, we recreate the distance matrix from our example:

import numpy as np from cogent.cluster import metric_scaling from numpy.linalg import eigh AB = 0.69276243 AC = 0.6549957 BC = 0.60625465 M = np.array( [0,AB,AC,AB,0,BC,AC,BC,0] ) M.shape = (3,3) print 'M\n', M, '\n'

M [[ 0. 0.69276243 0.6549957 ] [ 0.69276243 0. 0.60625465] [ 0.6549957 0.60625465 0. ]]

Now, we have two simple arithmetic operations on the data:

# step 1 E = metric_scaling.make_E_matrix(M) print 'E\n', E print 'M**2/-2\n', M**2/-2, '\n'

E [[-0. -0.23995989 -0.21450968] [-0.23995989 -0. -0.18377235] [-0.21450968 -0.18377235 -0. ]] M**2/-2 [[-0. -0.23995989 -0.21450968] [-0.23995989 -0. -0.18377235] [-0.21450968 -0.18377235 -0. ]]

# step 2 F = metric_scaling.make_F_matrix(E) print 'F\n', F row_means = np.mean(E,axis=0) col_means = np.mean(E,axis=1) overall_mean = np.mean(E) F2 = E - row_means - col_means + overall_mean print 'F2\n', F2, '\n\n'

F [[ 0.16114818 -0.08905749 -0.07209069] [-0.08905749 0.14065662 -0.05159913] [-0.07209069 -0.05159913 0.12368982]] F2 [[ 0.16114818 -0.08905749 -0.07209069] [-0.08905749 0.14065662 -0.05159913] [-0.07209069 -0.05159913 0.12368982]]

# step 3 eigvals, eigvecs = eigh(F) print 'eigvals\n', eigvals print 'eigvecs\n', eigvecs print eigval_root = np.sqrt(abs(eigvals)) print 'eigval_root\n', eigval_root print t_eigvecs = eigvecs.transpose() print 't_eigvecs\n', t_eigvecs print print t_eigvecs * eigval_root[:,np.newaxis]

As you can see, the eigenvalues do not come out in sorted order. I have copied the code that UniFrac uses to massage the vectors from this point, but I'm not sure I understand it yet---and not sure there isn't a simpler way :).

eigvals [ 5.20417043e-18 1.80259579e-01 2.45235038e-01] eigvecs [[-0.57735027 0.18984739 0.79411878] [-0.57735027 0.59280334 -0.56147205] [-0.57735027 -0.78265073 -0.23264672]] eigval_root [ 2.28126509e-09 4.24569875e-01 4.95212115e-01] t_eigvecs [[-0.57735027 -0.57735027 -0.57735027] [ 0.18984739 0.59280334 -0.78265073] [ 0.79411878 -0.56147205 -0.23264672]] [[ -1.31708902e-09 -1.31708902e-09 -1.31708902e-09] [ 8.06034836e-02 2.51686440e-01 -3.32289924e-01] [ 3.93257240e-01 -2.78047763e-01 -1.15209477e-01]]

The final result bears some resemblance to what fast_unifrac gave us (see last time).

================================================================ Type Label vec_num-0 vec_num-1 vec_num-2 ---------------------------------------------------------------- Eigenvectors A 0.39 0.08 -0.00 Eigenvectors B -0.28 0.25 -0.00 Eigenvectors C -0.12 -0.33 -0.00 Eigenvalues eigenvalues 0.25 0.18 0.00 Eigenvalues var explained (%) 57.64 42.36 0.00 ----------------------------------------------------------------

Can you see that in the final output from the script, we have the x-values for each sample in order in the third row of the vector, and the y-values in the second row? The 3rd dimension is in the first row.

While there are still plenty of mysteries, I believe we have at least a handle on what is going on here. That's enough for now.

Unifrac (5): results using PyCogent

Last time, I used a set of simulated sequences and their phylogenetic tree as input to the web interface of Unifrac (the original version). This time, I propose to repeat the analysis on my laptop using PyCogent. (See the PyCogent category in the sidebar for more from me on this, or start with this post). I'm just going to do the analysis here. Next time I want to look the distinction between PCA and PCoA.

Provided that we've done our imports, and loaded both the tree and the environment, the actual analysis runs in a single line. The results are shown first, followed by the code. The commented lines in the code were used in the first implementation based on the usage example, but aren't needed here.

There are a couple of things to notice. As for the rest of PyCogent, objects bearing results act like dicts, and we recover the elements of interest with result['pcoa'] (for example). Second, the distance matrix (a table of UniFrac metrics for pairwise comparisons between each of the three samples) is slightly different than what we got from the web interface. The comparison of A with C gives exactly the same result, but the other two have changed. Therefore, it is not surprising that the PCoA results are different. It is however, a bit of a shock to see that the small change in the ED (Environmental Distance) matrix has rotated the plot through about 170°.

[UPDATE: I forgot to mention that I had altered PyCogent slightly to always do PCoA. I think what you will need to do with this code is to pass an argument to fast_unifrac of modes = ['distance_matrix','pcoa'], but I haven't tested it yet.]

(array([[ 0. , 0.69276243, 0.6549957 ], [ 0.69276243, 0. , 0.60625465], [ 0.6549957 , 0.60625465, 0. ]]), ['A', 'B', 'C']) ================================================================ Type Label vec_num-0 vec_num-1 vec_num-2 ---------------------------------------------------------------- Eigenvectors A 0.39 0.08 -0.00 Eigenvectors B -0.28 0.25 -0.00 Eigenvectors C -0.12 -0.33 -0.00 Eigenvalues eigenvalues 0.25 0.18 0.00 Eigenvalues var explained (%) 57.64 42.36 0.00 ----------------------------------------------------------------

from cogent import LoadTree #from cogent.parse.tree import DndParser #from cogent.maths.unifrac.fast_tree import UniFracTreeNode from cogent.maths.unifrac.fast_unifrac import fast_unifrac # convert our environment file repr back into a dict FH = open('environ.txt','r') data = FH.read().strip() FH.close() L = data.split('\n') env_dict = dict() for e in L: seq_name, sample = e.strip().split() env_dict[seq_name] = {sample:1} tr = LoadTree('samples.tree') #tree_str = tr.getNewick(with_distances=True) #tr = DndParser(tree_str, UniFracTreeNode) result = fast_unifrac(tr, env_dict) print result['distance_matrix'] print result['pcoa'] ''' R code: x = c(0.39,-0.28,-0.12) y = c(0.08,0.25,-0.33) plot(x,y,pch=c('A','B','C'),cex=2) '''

Unifrac (4): Web results

I'm exploring UniFrac in a series of posts (first here). Last time, we produced a set of simulated sequences from three samples, or environments. I ended up fiddling with that a bit more. I changed this line in my code to drop two of the A sequences out (A11 and A12 in the original example). I did this because the initial set did not show any significant difference between samples.

A = {0:5,1:5,2:0,3:0,4:0,5:0,6:1,7:1}

(Because the sequences are numbered in order, in the new set the labels A11 and A12 have been reassigned to the previous A13 and A14). The sequences are organized to make a phylogenetic tree which was written to disk as 'samples.tree'. I also wrote a simple script to process the titles of the sequences and save them in a file named 'environ.txt'. The first line is:

A1 A

Now let's use these two files as input to UniFrac (web interface). We compute the UniFrac metric or distance for each pair of samples. As you can see, we're going to do 1000 randomizations (this requires simple registration with the site):

and evaluate the signficance of the results by repeated randomization of the labels. The raw p-values

need to be corrected by multiplying by the total number of comparisons (Bonferroni correction).

Notice that only A and B are significantly different, although it was a near-run thing with A and C. We download the actual UniFrac values for comparison with what PyCogent gives us. The text version of the file looks like this:

. A B C A 0 0.666882546144 0.654995703131 B 0.666882546144 0 0.576182797703 C 0.654995703131 0.576182797703 0

We also ask for a PCA (Principal Coordinates Analysis). Here is the plot of the two principal coordinates:

And here are the eigenvalues and eigenvectors in a data file:

pc vector number 1 2 3 A -0.396005816461 -0.0214251320143 -3.04168679166e-09 B 0.22015648746 -0.276525789753 -3.04168679166e-09 C 0.175849329001 0.297950921767 -3.04168679166e-09 eigvals 0.236212472152 0.165700300462 2.77555756156e-17 percent variation explained 58.7720740039 41.2279259961 6.90587050397e-15

Of course, we can always replot using R:

> x = c(-0.396,0.220,0.176) > y = c(-0.021,-0.277,0.298) > plot(x,y,pch=c('A','B','C'),cex=2)

You can make this as fancy as you like.

Unifrac (3): Simulating sequences

As part of my exploration of Unifrac (first and second posts), I'm going to need some sequences.

I posted about this before, then pulled that post. In fact, I pushed out this very post yesterday, and am now editing heavily. I'm finding a conflict between what's needed to make a good example for showing how the calculations work, and what's needed for a good example in terms of the final result. I've decided to keep this here, but be aware that we may not use these sequences going forward.

This script fetches a bunch of prototype sequences from Genbank. Then it distributes them to three "samples" according to a hard-coded distribution. In the process, the sequences are mutagenized a bit. Finally, the sequences are aligned, rooted (by hand) and then, as shown above, a plot is made in R of the phylogenetic tree.

[UPDATE: I modified the code slightly to give a bit more heterogeneity to the sequences. The code is updated, and here is a plot of the new phylogenetic tree.

Also, I fixed an awkward point about the old code: I have made sure that the tree we use is a real, rooted tree by incorporating Thermotoga as an outgroup.]



Here's the code (R code at the end):

import random, os, sys random.seed(137) from cogent import LoadSeqs, DNA, LoadTree from cogent.db.ncbi import EFetch from cogent.app.muscle import align_unaligned_seqs from cogent.app.fasttree import build_tree_from_alignment def fetch_ncbi_data(ofile,s): # get the seqs from Genbank input = [e.split() for e in s.strip().split('\n')] id_list = [t[0] for t in input] names = [t[1] for t in input] ef = EFetch(id=','.join(id_list), rettype='fasta') data = ef.read().strip() # title lines are too long, replace by genus_species rL = list() for i,e in enumerate(data.split('\n\n')): old_title, seq = e.strip().split('\n',1) new_title = '>' + names[i] seq = seq[:500] rL.append('\n'.join([new_title,seq])) FH = open(ofile,'w') FH.write('\n\n'.join(rL)) FH.close() def mutagenize(seq, mrate=1): L = list(seq) D = { 'A':'CGT', 'C':'AGT', 'G':'ACT', 'T':'ACG' } N = int(mrate / 100.0 * len(seq)) X = len(seq) for i in range(N): j = random.choice(range(X)) nt = L[j] if not nt in 'ACGT': continue L[j] = random.choice(D[nt]) return ''.join(L) def distribute_seqs(ifile,ofile): # set up our samples FH = open(ifile,'r') data = FH.read().strip().split('\n\n') FH.close() seqs = list() for e in data: title,seq = e.split('\n',1) seqs.append(''.join(seq.split())) outgroup = '>Thermotoga\n' + seqs.pop() A = {0:5,1:5,2:0,3:1,4:0,5:1,6:1,7:1} # A has lots of Firmicutes B = {0:0,1:1,2:5,3:5,4:1,5:0,6:1,7:1} # B has Bacteroidetes C = {0:1,1:0,2:1,3:0,4:5,5:5,6:1,7:1} # C has enterics dL = [A,B,C] L = list() for distr, sample in zip(dL,list('ABC')): counter = 1 for k in distr: seq = seqs[k] n = distr[k] for i in range(n): if n == 1: mrate = 5 else: mrate = random.choice((1,2,3)) copy = mutagenize(seq[:],mrate) name = sample + str(counter) L.append(DNA.makeSequence(copy,name)) counter += 1 FH = open(ofile,'w') L = [seq.toFasta() for seq in L] L.insert(0,outgroup) FH.write('\n\n'.join(L)) FH.close() def align_seqs(ifile,ofile): seqs = LoadSeqs(ifile, moltype=DNA, aligned=False) aln = align_unaligned_seqs(seqs, DNA) aln.writeToFile(ofile) return aln def get_tree(ifile): aln = LoadSeqs(ifile, moltype=DNA, aligned=True) tr = build_tree_from_alignment(aln,moltype=DNA) return tr #=============================================== s = ''' AY005045.1 Streptococcus_mitis_bv2 D83363.1 Staphylococcus_epidermidis_14990 L14639.1 Capnocytophaga_gingivalis AB053940.1 Tannerella_forsythensis_HA3 EU009197.1 Shigella_sonnei_FBD023 AB435616.1 Serratia_marcescens_JCM24201 AB302401.1 Pseudomonas_cinnamophila AF411020.1 Achromobacter_xylosoxidans_AU1011 AJ401017.1 Thermotoga_maritima_SL7 ''' fn1 = 'rRNA_gb.fasta' fn2 = 'samples.fasta' fn3 = 'samples.aln.fasta' fn4 = 'samples.tree' if not os.path.exists(fn1) or False: fetch_ncbi_data(fn1,s) if not os.path.exists(fn2) or False: distribute_seqs(fn1,fn2) if not os.path.exists(fn3) or True: aln = align_seqs(fn2,fn3) tr = get_tree(fn3) # re-root manually print tr.asciiArt() n = tr.getNodeMatchingName('Thermotoga') for a in n.ancestors(): print a.Name tr2 = tr.rootedAt(n.ancestors()[0].Name) tree_str = tr2.getNewick(with_distances=True) FH = open(fn4,'w') FH.write(tree_str + '\n') FH.close() tr = LoadTree(fn4) print tr.asciiArt() ''' R code: library(ape) setwd('Desktop') tr = read.tree('temp/samples.tree') colors = rep('black',length(tr$tip.label)) o = grep('A',tr$tip.label) colors[o] = 'red' o = grep('B',tr$tip.label) colors[o] = 'blue' o = grep('C',tr$tip.label) colors[o] = 'darkgreen' plot(tr,tip.color=colors) axisPhylo() '''

Unifrac analysis (2): P-test

I'm planning to discuss UniFrac in a series of posts (first here). But before moving further into it, I would like to define another statistic called the Phylogenetic-test or P-test. This test is available in UniFrac, but it was introduced earlier by Andrew Martin (Martin 2002 PMID 12147459).

Consider the same phylogenetic tree as in the first post. Here, the environment (sample) labels for a few of the sequences have been exchanged to provide a bit more heterogeneity in the tree.

The criterion of parsimony is employed to characterize the tree with respect to the labels, by computing the minimum number of "changes" (switches from one community to another), to explain the observed distribution. For example, we can imagine (in some theoretical sense) to have started with a blue label, and switched to red (and then to blue again) as indicated on the left; or alternatively, to have started with red and changed to blue, on the right. Here, the minimum number of changes needed to "explain" the tree is four.

The significance of the result is evaluated by label randomization (as for UniFrac itself), and the result is referred to as a P-test. If the randomization tests yield a statistic as extreme as that for the actual sample less than 5% of the time, then the distribution is considered to be significant. For a sample where one environment is distinct from another, fewer label switches are required. With completely separate distributions, the number of changes needed is one.

Unifrac analysis: introduction

I'm going to take a few posts to explore the "Unifrac tool for comparing microbial community diversity in a phylogenetic context." UniFrac was developed by Rob Knight and Cathy Lozupone and their colleagues at the University of Colorado, and has become an important statistical apprach in microbial ecology.

The methodological papers are:

• Lozupone
2005 PMID 16332807

• Lozupone
2006 PMID 16893466

• Hamady
2010 PMID 19710709

(The last one is behind a firewall and likely to remain there, unfortunately).

Both the original version of UniFrac, and a new implementation called FastUniFrac (Hamady paper above) have web interfaces. FastUniFrac is also part of PyCogent.

The specific goal is to understand the population structure of microbial communities. In "classical" ecology, we might count how many species are present in a given environment ("lions and tigers and bears"---to quote the Tin Woodsman), and the numbers of each. In microbial ecology, the data we have are DNA sequences, typically, 16S rRNA gene sequences. Since these are homologous sequences, we can obtain additional information including the phylogenetic distances between different individuals. UniFrac makes use of phylogenetic distance in quantifying community structure.

There is the complication of deciding when two sequences belong to the same "species" or OTU (operational taxonomic unit), which we will ignore for the moment.

Knight et al introduced the UniFrac "metric." Consider two idealized phylogenetic trees. In the first one, above left, the sequences from two environments (red and blue) are cleanly separated into monophyletic clusters. We count the amount of "unique branch length," the fraction of total branch length that leads exclusively to one or the other type of sequence, shown in either red or blue color in the figure.

In contrast, the second tree (right) has the identical structure but the environment labels have been switched for one pair of sequences. Now, a substantial fraction of branch length on the tree is not unique to one or the other label (shared branches are colored black in the figure).

The statistical significance of the UniFrac metric is determined by a resampling technique, randomizing the labels on sequences (in the same proportion as in the original sample). In a large number of replicates (say, 1000), the method counts the number of times the statistic is more extreme than for the observed tree. If this number is < 0.05, we consider the result significant.

One point that I'm pretty sure but not certain about:
• randomized samples contain the same total of blue or red sequences as in the original.

[UPDATE: an earlier version of this post had an error. The correct statement is that:
• the lengths of the terminal branches are considered when doing the calculation.]