Wolfram Function Repository

Instant-use add-on functions for the Wolfram Language

Function Repository Resource:

 PhylogeneticTreePlot

Source Notebook

Plot a dendrogram for a set of genome nucleotide sequences

Contributed by: Daniel Lichtblau

ResourceFunction["PhylogeneticTreePlot"][seqs,slabels]

gives a dendrogram plot that groups the genetic sequences seqs and labels them with slabels.

Details and Options

ResourceFunction["PhylogeneticTreePlot"] uses an alignment-free method to compare pairs of sequences.
Input sequences should be strings comprised of the standard nucleotide letters {A,C,G,T}. Lowercase letters are also allowed. The character U is allowed and is replaced by T. All other characters, such as N, are removed.
ResourceFunction["PhylogeneticTreePlot"] uses Dendrogram and accepts all options for that function.
Most default option settings agree with those of Dendrogram.
Unlike Dendrogram, the DistanceFunction defaults to CosineDistance.
ResourceFunction["PhylogeneticTreePlot"] creates the dendrogram from vectors that are derived via dimensional reductions of the input genetic sequences.
By default, each sequence is first converted to an image using the Frequency Chaos Game Representation (FCGR).
FCGR images are reduced in dimension using a Fourier Cosine Transform (FCT). A further dimensional reduction is done using the Singular Value Decomposition (SVD) on vectors comprised of the flattened FCT matrices.
An explanation may be found in the articles noted in the Related Links and Source Metadata sections.
The dimensional reduction described above has certain parameters that in principle might be changed. ResourceFunction["PhylogeneticTreePlot"] uses a fixed set of values for these that has been seen to perform fairly well in practice.
In order to produce a reasonable grouping, ResourceFunction["PhylogeneticTreePlot"] requires moderately long genetic sequences containing at least a few thousand nucleotides.
The SVD step of the dimension reduction, when applied to vectors that came from a set of very similar genome sequences, will tend to give a result where the reduced vectors are dominated by a large first value that is approximately equal across the set. This tends to distort the distances between genomes. The option "IgnoreOutlier" (default: Automatic) is provided to address this. The automatic behavior is to remove the first components whenever the largest singular value is at least ten times the size of the next largest.
ResourceFunction["PhylogeneticTreePlot"] takes an option "VectorMethod". The Automatic default converts sequences to vectors using the methodology described above. A user-defined method can be specified with this option. One possibility is to use the resource function StringsToVectors.

Examples

Basic Examples (4) 

Obtain genes for a type of breast cancer from the entity store:

In[1]:=
InputForm[ brca2 = SemanticInterpretation["brca2 gene", AmbiguityFunction -> All]]
Out[361]=
In[362]:=
bcGenes = brca2[[1]]
Out[362]=

Obtain the nucleotide sequences and species names:

In[363]:=
refseqs = Through[bcGenes["ReferenceSequence"]]; species = Through[bcGenes["Species"]]
Out[363]=

Plot the phylogenetic tree derived from these genes:

In[364]:=
ResourceFunction["PhylogeneticTreePlot"][refseqs, species]
Out[364]=

One can force PhylogeneticTreePlot to ignore the largest singular value in the dimension-reduction:

In[365]:=
ResourceFunction["PhylogeneticTreePlot"][refseqs, species, "IgnoreOutlier" -> True]
Out[365]=

Neat Examples (10) 

Mammal benchmark (5) 

This example was introduced in the referenced article "A novel fast vector method for genetic sequence comparison" by Li, He, He and Yau and also used in "Alignment-free genomic sequence comparison using FCGR and signal processing".

Show an example that has appeared in literature on alignment-free methods for clustering gene sequences:

In[366]:=
genomesMammal = {"V00662", "D38116", "D38113", "D38114", "X99256", "Y18001", "AY863426", "D38115", "NC_002083", "U20753", "U96639", "AJ002189", "AF010406", "AF533441", "V00654", "AY488491", "EU442884", "EF551003", "EF551002", "X97336", "Y07726", "DQ402478", "AF303110", "AF303111", "EF212882", "AJ001588", "X88898", "NC_002764", "AJ238588", "AJ001562", "X72204", "NC_005268", "NC_007441", "NC_008830", "NC_001788", "NC_001321", "NC_005270", "NC_001640", "NC_005275", "NC_006931", "NC_010640"}; speciesMammal = {"Human", "Pygmy chimpanzee", "Common chimpanzee", "Gorilla", "Gibbon", "Baboon", "Vervet monkey", "Bornean orangutan", "Sumatran orangutan", "Cat", "Dog", "Pig", "Sheep", "Goat", "Cow", "Buffalo", "Wolf", "Tiger", "Leopard", "Indian rhinoceros", "White rhinoceros", "Black bear", "Brown bear", "Polar bear", "Giant panda", "Rabbit", "Hedgehog", "Macaca thibet", "Squirrel", "Dormouse", "Blue whale", "Bowhead whale", "Chiru", "Common warthog", "Donkey", "Fin whale", "Gray whale", "Horse", "Indus river dolphin", "North pacific right whale", "Taiwan serow"}; stringColor[str_String] := Which[ StringContainsQ[str, {"Human", "chimpanzee", "Gorilla", "orangutan"}], Red, StringContainsQ[str, {"Gibbon", "Baboon", "monkey", "thibet"}], Orange, StringContainsQ[str, {"Tiger", "Leopard", "Cat"}], Darker[Green], StringContainsQ[str, {"Dog", "Wolf"}], Lighter[Green], StringContainsQ[ str, {"Cow", "Buffalo", "Goat", "Sheep", "Chiru", "serow"}], Darker[Brown], StringContainsQ[str, "rhino"], Lighter[Brown], StringContainsQ[str, {"bear", "panda"}], Black, StringContainsQ[str, {"mouse", "Squirrel", "Rabbit"}], Gray, StringContainsQ[str, {"Pig", "warthog"}], Darker[Blue], StringContainsQ[str, {"Donkey", "Horse"}], Blue, StringContainsQ[str, {"whale", "dolphin"}], Purple, StringContainsQ[str, "Hedgehog"], Darker[Gray]] speciesMammalColored = Map[Style[#, stringColor[#]] &, speciesMammal]
Out[369]=

Use the resource function ImportFASTA to download the genetic sequences from GenBank for the previously shown accession identifiers:

In[370]:=
getFASTA = ResourceFunction["ImportFASTA"]; AbsoluteTiming[ sequencesMammal = Map[getFASTA, genomesMammal][[All, 2, 1]];]
Out[371]=

Plot the phylogenetic tree:

In[372]:=
ResourceFunction[ "PhylogeneticTreePlot"][sequencesMammal, speciesMammalColored, PlotLabel -> "Mammalian mitochondrial DNA"]
Out[372]=

Use a non-default method for clustering:

In[373]:=
ResourceFunction[ "PhylogeneticTreePlot"][sequencesMammal, speciesMammalColored, PlotLabel -> "Mammalian mitochondrial DNA", ClusterDissimilarityFunction -> "Average"]
Out[373]=

Use the default n-gram based method of the resource function StringsToVectors to convert the sequences:

In[374]:=
ResourceFunction[ "PhylogeneticTreePlot"][sequencesMammal, speciesMammalColored, "VectorMethod" -> ResourceFunction["StringsToVectors"], PlotLabel -> "Mammalian mitochondrial DNA"]
Out[374]=

SARS-CoV-2 (5) 

We now apply this to samples of the novel coronavirus SARS-CoV-2, which is the virus responsible for the COVID-19 pandemic. For purposes of getting better comparisons, we remove trailing A nucleotides (their presence or absence seems to be a vagary of the sequencing protocols). We also restrict to a tighter size range that still contains most of the fully sequenced samples. For purposes of visualization we will color the labels by location and date of sample collection.

Obtain data from the Wolfram Data Repository:

In[375]:=
fullGenes = ResourceData["Genetic Sequences for the SARS-CoV-2 Coronavirus"][ Select[StringContainsQ[#GenBankTitle, "complete genome"] &], {"Sequence", "GeographicLocation", "CollectionDate"}]; geneLists = Normal[fullGenes];

Remove trailing A nucleotides and retain only a narrow size range:

In[376]:=
dropTrailingA[seq_] := StringReplace[seq, StartOfString ~~ Shortest[a__] ~~ ("A" ..) ~~ EndOfString :> a]; normalGL = Map[Values, geneLists]; midsize = Select[normalGL, 29825 <= StringLength[dropTrailingA[#[[1]]]] <= 29885 &]; {fullseqs0, places, dates} = Transpose[midsize]; fullseqs = Map[dropTrailingA, fullseqs0];

Create labels colored both by collection locations and dates:

In[377]:=
dates1 = Map[First, dates]; ndates1 = Map[(#[[1]] - 2019)*365 + #[[2]]*30.4 + If[Length[#] == 3, #[[3]], 15] &, dates1]; ndates = Rescale[ndates1]; colorDateLabels = Thread[Style[dates1, Map[Darker[Hue[#]] &, ndates + .05]]]; placeNames = Map[If[Head[#] === Entity, #[[2]], Head[#]] &, places]; placeSet = Union[placeNames]; nplaces = Length[placeSet]; crange = Range[nplaces]/(nplaces + 1); biggest = Ordering[ Map[Length, Map[ColorData[#, "ColorList"] &, ColorData["Indexed"]]], -1]; placeColors = Take[ColorData[ColorData["Indexed"][[biggest[[1]]]], "ColorList"], UpTo[nplaces]]; colorPlaceRule = Thread[placeSet -> placeColors]; colorPlaceLabels = Thread[Style[placeNames, placeNames /. colorPlaceRule]]; colorPlaceDateLabels = Thread[{colorPlaceLabels, colorDateLabels}];

Show the phylogenetic tree:

In[378]:=
ResourceFunction[ "PhylogeneticTreePlot"][fullseqs, colorPlaceDateLabels, AspectRatio -> 5, ImageSize -> 750]
Out[378]=

This example illustrates the importance of clipping off the first singular value in the dimension reduction, in cases where it is much larger than the next one. In effect one component will dominate the rest, in a way that will tend to make the distances less meaningful. We show this now by forcing the large value to be retained. The result is a tree that is unbalanced in terms of branch lengths:

In[379]:=
ResourceFunction[ "PhylogeneticTreePlot"][fullseqs, colorPlaceDateLabels, "IgnoreOutlier" -> False, AspectRatio -> 5, ImageSize -> 750]
Out[379]=

Possible Issues (3) 

The implementation was changed in April 2020 to use the System` context Dendrogram instead of HierarchicalClustering`DendrogramPlot. This necessitated small changes to the argument structure and also means there are some changes in the functional behavior. For example, Dendrogram did not at that time properly handle CosineDistance as a DistanceFunction setting. We give the prior implementation here as a convenience for those who might have preferred it:

In[380]:=
Options[OldPhylogeneticTreePlot] = Join[Options[ HierarchicalClustering`DendrogramPlot], {AspectRatio -> 1.2, ImageSize -> 600, "IgnoreOutlier" -> Automatic}] /. (HierarchicalClustering`Orientation -> aa_) :> (HierarchicalClustering`Orientation -> Left); OldPhylogeneticTreePlot[geneseqs_, labels_, opts : OptionsPattern[]] := OldPhylogeneticTreePlot[geneseqs, labels, "", opts] OldPhylogeneticTreePlot[geneseqs_, labels_, plabel_, opts : OptionsPattern[]] := With[{freq = 30, keep = 40, dim = 7}, Module[ {ftts, vecs, ratio, size, io, optsd}, {ratio, size, io} = ({AspectRatio, ImageSize, "IgnoreOutlier"} /. {opts}) /. Thread[{AspectRatio, ImageSize, "IgnoreOutlier"} -> {OptionValue[AspectRatio], OptionValue[ImageSize], OptionValue["IgnoreOutlier"]}]; Block[{$ContextPath}, Needs["HierarchicalClustering`"]]; optsd = FilterRules[{opts}, Options[HierarchicalClustering`DendrogramPlot]]; ftts = Map[Developer`ToPackedArray[ processNucleotideString[#, dim, freq]] &, geneseqs]; vecs = FTTsToVectors[ftts, keep, io]; HierarchicalClustering`DendrogramPlot[vecs, HierarchicalClustering`LeafLabels -> labels, PlotLabel -> plabel, AspectRatio -> ratio, ImageSize -> size, Apply[Sequence, optsd], Apply[Sequence, Options[ ResourceFunction["PhylogeneticTreePlot"]]]] ]]

Set up the mammal benchmark:

In[381]:=
genomesMammal = {"V00662", "D38116", "D38113", "D38114", "X99256", "Y18001", "AY863426", "D38115", "NC_002083", "U20753", "U96639", "AJ002189", "AF010406", "AF533441", "V00654", "AY488491", "EU442884", "EF551003", "EF551002", "X97336", "Y07726", "DQ402478", "AF303110", "AF303111", "EF212882", "AJ001588", "X88898", "NC_002764", "AJ238588", "AJ001562", "X72204", "NC_005268", "NC_007441", "NC_008830", "NC_001788", "NC_001321", "NC_005270", "NC_001640", "NC_005275", "NC_006931", "NC_010640"}; speciesMammalColored = {\!\(\* StyleBox["\"\<Human\>\"", StripOnInput->False, FontColor->RGBColor[1, 0, 0], "NodeID" -> 110]\), \!\(\* StyleBox["\"\<Pigmy chimpanzee\>\"", StripOnInput->False, FontColor->RGBColor[1, 0, 0], "NodeID" -> 111]\), \!\(\* StyleBox["\"\<Common chimpanzee\>\"", StripOnInput->False, FontColor->RGBColor[1, 0, 0], "NodeID" -> 112]\), \!\(\* StyleBox["\"\<Gorilla\>\"", StripOnInput->False, FontColor->RGBColor[1, 0, 0], "NodeID" -> 113]\), \!\(\* StyleBox["\"\<Gibbon\>\"", StripOnInput->False, FontColor->RGBColor[1, Rational[1, 3], Rational[1, 3]], "NodeID" -> 114]\), \!\(\* StyleBox["\"\<Baboon\>\"", StripOnInput->False, FontColor->RGBColor[1, Rational[1, 3], Rational[1, 3]], "NodeID" -> 115]\), \!\(\* StyleBox["\"\<Vervet monkey\>\"", StripOnInput->False, FontColor->RGBColor[1, Rational[1, 3], Rational[1, 3]], "NodeID" -> 116]\), \!\(\* StyleBox["\"\<Bornean orangutan\>\"", StripOnInput->False, FontColor->RGBColor[1, 0, 0], "NodeID" -> 117]\), \!\(\* StyleBox["\"\<Sumatran orangutan\>\"", StripOnInput->False, FontColor->RGBColor[1, 0, 0], "NodeID" -> 118]\), \!\(\* StyleBox["\"\<Cat\>\"", StripOnInput->False, FontColor->RGBColor[0., 0., 0.65], "NodeID" -> 119]\), \!\(\* StyleBox["\"\<Dog\>\"", StripOnInput->False, FontColor->RGBColor[0, 0, 1], "NodeID" -> 120]\), \!\(\* StyleBox["\"\<Pig\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 121]\), \!\(\* StyleBox["\"\<Sheep\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 122]\), \!\(\* StyleBox["\"\<Goat\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 123]\), \!\(\* StyleBox["\"\<Cow\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 124]\), \!\(\* StyleBox["\"\<Buffalo\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 125]\), \!\(\* StyleBox["\"\<Wolf\>\"", StripOnInput->False, FontColor->RGBColor[0, 0, 1], "NodeID" -> 126]\), \!\(\* StyleBox["\"\<Tiger\>\"", StripOnInput->False, FontColor->RGBColor[0., 0., 0.65], "NodeID" -> 127]\), \!\(\* StyleBox["\"\<Leopard\>\"", StripOnInput->False, FontColor->RGBColor[0., 0., 0.65], "NodeID" -> 128]\), \!\(\* StyleBox["\"\<Indian rhinoceros\>\"", StripOnInput->False, FontColor->RGBColor[0.36, 0.24, 0.12], "NodeID" -> 129]\), \!\(\* StyleBox["\"\<White rhinoceros\>\"", StripOnInput->False, FontColor->RGBColor[0.36, 0.24, 0.12], "NodeID" -> 130]\), \!\(\* StyleBox["\"\<Black bear\>\"", StripOnInput->False, FontColor->RGBColor[0., 0., 0.30000000000000004`], "NodeID" -> 131]\), \!\(\* StyleBox["\"\<Brown bear\>\"", StripOnInput->False, FontColor->RGBColor[0., 0., 0.30000000000000004`], "NodeID" -> 132]\), \!\(\* StyleBox["\"\<Polar bear\>\"", StripOnInput->False, FontColor->RGBColor[0., 0., 0.30000000000000004`], "NodeID" -> 133]\), \!\(\* StyleBox["\"\<Giant panda\>\"", StripOnInput->False, FontColor->RGBColor[0., 0., 0.30000000000000004`], "NodeID" -> 134]\), \!\(\* StyleBox["\"\<Rabbit\>\"", StripOnInput->False, FontColor->RGBColor[0.2, 0.2, 0.2], "NodeID" -> 135]\), \!\(\* StyleBox["\"\<Hedgehog\>\"", StripOnInput->False, FontColor->RGBColor[0.4, 0.4, 0.4], "NodeID" -> 136]\), \!\(\* StyleBox["\"\<Macaca thibet\>\"", StripOnInput->False, FontColor->RGBColor[1, Rational[1, 3], Rational[1, 3]], "NodeID" -> 137]\), \!\(\* StyleBox["\"\<Squirrel\>\"", StripOnInput->False, FontColor->GrayLevel[0], "NodeID" -> 138]\), \!\(\* StyleBox["\"\<Dormouse\>\"", StripOnInput->False, FontColor->GrayLevel[0], "NodeID" -> 139]\), \!\(\* StyleBox["\"\<Blue whale\>\"", StripOnInput->False, FontColor->RGBColor[0, Rational[2, 3], 0], "NodeID" -> 140]\), \!\(\* StyleBox["\"\<Bowhead whale\>\"", StripOnInput->False, FontColor->RGBColor[0, Rational[2, 3], 0], "NodeID" -> 141]\), \!\(\* StyleBox["\"\<Chiru\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 142]\), \!\(\* StyleBox["\"\<Common warthog\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 143]\), \!\(\* StyleBox["\"\<Donkey\>\"", StripOnInput->False, FontColor->RGBColor[0.36, 0.24, 0.12], "NodeID" -> 144]\), \!\(\* StyleBox["\"\<Fin whale\>\"", StripOnInput->False, FontColor->RGBColor[0, Rational[2, 3], 0], "NodeID" -> 145]\), \!\(\* StyleBox["\"\<Gray whale\>\"", StripOnInput->False, FontColor->RGBColor[0, Rational[2, 3], 0], "NodeID" -> 146]\), \!\(\* StyleBox["\"\<Horse\>\"", StripOnInput->False, FontColor->RGBColor[0.36, 0.24, 0.12], "NodeID" -> 147]\), \!\(\* StyleBox["\"\<Indus river dolphin\>\"", StripOnInput->False, FontColor->RGBColor[0, Rational[2, 3], 0], "NodeID" -> 148]\), \!\(\* StyleBox["\"\<North pacific right whale\>\"", StripOnInput->False, FontColor->RGBColor[0, Rational[2, 3], 0], "NodeID" -> 149]\), \!\(\* StyleBox["\"\<Taiwan serow\>\"", StripOnInput->False, FontColor->RGBColor[0.7333333333333333, 0.6, 0.4666666666666667], "NodeID" -> 150]\)}; getFASTA = ResourceFunction["ImportFASTA"]; AbsoluteTiming[ sequencesMammal = Map[getFASTA, genomesMammal][[All, 2, 1]];]
Out[381]=

Now show the old method, using the cosine distance:

In[382]:=
OldPhylogeneticTreePlot[sequencesMammal, speciesMammalColored, "Mammalian mitochondrial DNA", DistanceFunction -> CosineDistance]
Out[382]=

Version History

  • 4.0.1 – 22 March 2024
  • 4.0.0 – 09 April 2020
  • 3.0.0 – 18 February 2020
  • 2.0.0 – 31 December 2019
  • 1.0.0 – 20 August 2019

Source Metadata

  • Citation:
    • Lichtblau, D., "Linking Fourier and PCA methods for image look-up." In: Proceedings of the 2016 18th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC). Piscataway, NJ: IEEE, 105–110, 2016. DOI: 10.1109/SYNASC.2016.028
    • Lichtblau, D., "Alignment-free genomic sequence comparison using FCGR and signal processing." BMC Bioinformatics, vol. 20, no. 742, 2019. DOI: 10.1186/s12859-019-3330-3
    • Li, Y., He, L., He, R. L., Yau, S. S.-T., "A novel fast vector method for genetic sequence comparison." Scientific Reports, vol. 7, no. 12226, 2017. DOI: 10.1038/s41598-017-12493-2
    • Clark, K., Karsch-Mizrachi, I., Lipman, D. J., Ostell, J., Sayers, E. W., "GenBank." Nucleic Acids Research, vol. 44, no. D67-D72, 2016. DOI: 10.1093/nar/gkv1276

Related Resources

License Information