# Hierarchical clustering: convert Mathematica cluster hierarchy to NEWICK format

Can anyone tell me a clever way to convert the hierarchical clustering object, created, say, by Agglomerate or DirectAgglomerate (of the built-in hierarchical clustering package) to the so-called Newick (or New Hampshire) format used in biology? It would convert

Cluster[Cluster[
Cluster["A", Cluster["H", "J", 1.52217, 1, 1], 28.8538, 1, 2],
Cluster[Cluster["C", "E", 10.1371, 1, 1], "D", 22.0063, 2, 1],
47.1129, 3, 3], Cluster[Cluster["B", Cluster["G", "I", 2.5374, 1, 1],
5.73533, 1, 2],"F", 13.6197, 3, 1], 64.5168, 6, 4]


to

(((A:28.85378,(H:1.52217,J:1.52217):27.33162):18.25912,
((C:10.13706,E:10.13706):11.86925,D:22.00630):25.10660):17.40389,
((B:5.73533,(G:2.53740,I:2.53740):3.19793):7.88433,F:13.61966):50.89713)


the latter composed of branch lengths rather than nodal abscissae.

-
There are some numbers in the first format that are not present in the second and vice versa. So, it's not just a bit of reshuffling and replacement but more. A more precise specification seems needed. – Sjoerd C. de Vries Jun 18 '14 at 6:17
In the "newick" format the values, which represent branch lengths, are pairwise differences of the nodal abscissae at the ends of each given branch. I.e., differences of the numbers in the Cluster objects. – user15994 Jun 18 '14 at 21:18

Until a cleaner approach is offered, consider the following brute-force method for a partial solution:

ClearAll[newickF];
Needs["HierarchicalClustering"]
newickF[clstr_] := Fold[Replace[#, #2, {0, Infinity}] &, clstr,
{Cluster[a_, b_, c_, _, _] :> {{a, c}, {b, c}},
{{{a_, r_}, {b_, r_}}, t_} :> {{{a, r}, {b, r}}, t - r},
{{{{a_, r_}, {b : {{_, s_}, {_, s_}}, t_}}, u_} :> {{{a, r}, {b, t}}, u - r},
{{{b : {{_, s_}, {_, s_}}, t_}, {a_, r_}},  u_} :> {{{b, t}, {a, r}}, u - r}},
{a_, b_?NumericQ} :> StringJoin[ToString[a], ":", ToString[b]]}]


OP's example:

clusters =  Cluster[Cluster[
Cluster["a", Cluster["h", "j", 1.52217, 1, 1], 28.8538, 1, 2],
Cluster[Cluster["c", "e", 10.1371, 1, 1], "d", 22.0063, 2, 1], 47.1129, 3, 3],
Cluster[Cluster["b", Cluster["g", "i", 2.5374, 1, 1], 5.73533, 1, 2],
"f", 13.6197, 3, 1], 64.5168, 6, 4];


In picture:

dplt1 = DendrogramPlot[clusters, LeafLabels -> (# &),
GridLines -> {None,
tt = Cases[clusters, Cluster[a_, b_, c_, d__] :> c, {0, Infinity}]},
GridLinesStyle -> Green, ImageSize -> 500, Axes -> {False, True},
AxesOrigin -> {.75, Automatic}, Ticks -> {Automatic, tt}]


newickF[clusters]
(* {"{{a:28.8538, {h:1.52217, j:1.52217}:27.3316}:18.2591,
{{c:10.1371, e:10.1371}:11.8692, d:22.0063}:25.1066}:17.4039",
"{{b:5.73533, {g:2.5374, i:2.5374}:3.19793}:7.88437, f:13.6197}:50.8971"} *)

rule = Line[a : {{x1_, y1_}, {x1_, y2_}, {x2_, z1_}, {x2_, z2_}}] :>
{Line[a], Directive[Thick, Opacity[.5], Red],
Line[{{x1, y1}, {x1, y2}}], Line[{{x2, z1}, {x2, z2}}], Opacity[1],
Text[Abs[y2 - y1], {x1 + .2, (y1 + y2)/2}, Automatic, {0, 1}],
Text[Abs[z1 - z2], {x2 + .2, (z1 + z2)/2}, Automatic, {0, 1}]}
dplt1 /. rule


Another example:

cl = Agglomerate[N@{1, 2, 10, 12, 3, 14, 15, 20, 26, 25, 27} -> CharacterRange["A", "K"],
DistanceFunction -> ManhattanDistance, Linkage -> "Centroid"]
(* Cluster[Cluster[Cluster["A", Cluster["B", "E", 1., 1, 1], 1.25, 1, 2],
Cluster[Cluster["C", "D", 2., 1, 1],
Cluster["F", "G", 1., 1, 1],  2.75, 2, 2], 9.24306, 3, 4],
Cluster[Cluster["J", Cluster["K", "I", 1., 1, 1], 1.25, 1, 2],
"H", 5.55556, 3, 1], 11.9209, 7, 4] *)

dplt2 = DendrogramPlot[cl, LeafLabels -> (# &),
GridLines -> {None,
tt = Cases[cl,  Cluster[a_, b_, c_, d__] :> c, {0, Infinity}]},
GridLinesStyle -> Green, ImageSize -> 500, Axes -> {False, True},
AxesOrigin -> {.75, Automatic}, Ticks -> {Automatic, tt}]


newickF[cl]
(* {"{{A:1.25, {B:1., E:1.}:0.25}:7.99306,
{{C:2., D:2.}:0.75, {F:1., G:1.}:1.75}:8.49306}:2.67786",
"{{J:1.25, {K:1., I:1.}:0.25}:4.30556, H:5.55556}:6.36536"} *)

dplt2 /. rule
`

-
Thank you this is very helpful, almost exactly what I need. But why in line 2 of the example immediately above is the value 8.49306 rather than 6.49306? (The latter is the branch length shown in the figure just above.) – user15994 Jun 18 '14 at 21:22
@user15994, it looks like i need to add another replacement rule - I will see if i can fix it.. – kglr Jun 19 '14 at 7:28