<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/using-taxa-and-metacoder-to-explore-taxonomy-of-vancouvers-trees/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.yihanwu.ca/post/vancouver-trees-part-3/"> R on YIHAN WU</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>From the <a rel="nofollow" href="https://www.yihanwu.ca/post/vancouver-trees-part-2" target="_blank">last post</a>, Vancouver has several common genus in its collection, such as <em>Prunus</em> and <em>Acer</em>. So rather than analyzing the diversity of Vancouver tree species on a species level, we could, with the help of some R packages, visualize the variety of taxonomic groups.</p>
<p>First, we format the Vancouver tree dataset similar to before to get the proper scientific name. This is the same formatting done in the previous post.</p>
<pre><code>library(dplyr)
library(purrr)
library(tidyr)</code></pre>
<pre><code># read in  trees shp
# read in  trees df
trees_df&lt;-read.csv(&quot;StreetTrees_CityWide.csv&quot;, stringsAsFactors = F)
trees_df&lt;-trees_df[trees_df$LATITUDE > 40 &amp; trees_df$LONGITUDE &lt; -100,]
trees_df$SCIENTIFIC_NAME &lt;- paste(trees_df$GENUS_NAME, trees_df$SPECIES_NM)</code></pre>
<pre><code>trees_df$GENUS_NAME &lt;- as.character(trees_df$GENUS_NAME)
trees_df$SPECIES_NAME &lt;- as.character(trees_df$SPECIES_NAME)

species_correction &lt;- function(genus, species){
  if(species == &quot;SPECIES&quot;){
    return(genus)
  } else if (length(grep(&quot;\\sX$&quot;, species)) > 0){
    return(paste(genus,&quot; x &quot;, unlist(strsplit(species, &quot; &quot;))[[1]]))
  } else {
    return(paste(genus, species))
  }
}
trees_df$SCIENTIFIC_NAME &lt;- map2_chr(trees_df$GENUS_NAME, trees_df$SPECIES_NAME, species_correction)
head(trees_df)</code></pre>
<pre><code>##   TREE_ID CIVIC_NUMBER       STD_STREET NEIGHBOURHOOD_NAME
## 1     589         4694         W 8TH AV    WEST POINT GREY
## 2     592         1799 W KING EDWARD AV        SHAUGHNESSY
## 3     594         2401         W 1ST AV          KITSILANO
## 5     596         2428         W 1ST AV          KITSILANO
## 6     598         2428         W 1ST AV          KITSILANO
## 7     599         1725        BALSAM ST          KITSILANO
##          ON_STREET ON_STREET_BLOCK STREET_SIDE_NAME ASSIGNED
## 1        BLANCA ST            2400             EVEN        N
## 2 W KING EDWARD AV            1700              MED        N
## 3         W 1ST AV            2400              ODD        N
## 5         W 1ST AV            2400             EVEN        N
## 6         W 1ST AV            2400             EVEN        N
## 7         W 1ST AV            2400             EVEN        N
##   HEIGHT_RANGE_ID DIAMETER DATE_PLANTED PLANT_AREA ROOT_BARRIER CURB
## 1               5     32.0     20180128          B            N    Y
## 2               1      3.0     20100104         30            N    Y
## 3               4     14.0     20180128          B            N    Y
## 5               2     11.0     20180128          B            N    Y
## 6               2      3.0     20180128          B            N    Y
## 7               2     13.5     20180128          B            N    Y
##   CULTIVAR_NAME GENUS_NAME SPECIES_NAME              COMMON_NAME LATITUDE
## 1                    ALNUS       RUGOSA           SPECKLED ALDER 49.26532
## 2                    ABIES      GRANDIS                GRAND FIR 49.24952
## 3                    PINUS        NIGRA            AUSTRIAN PINE 49.27096
## 5                   PRUNUS    SARGENTII SARGENT FLOWERING CHERRY 49.27083
## 6                   PRUNUS    SARGENTII SARGENT FLOWERING CHERRY 49.27084
## 7                    PICEA        ABIES            NORWAY SPRUCE 49.27083
##   LONGITUDE  SCIENTIFIC_NAME
## 1 -123.2150     ALNUS RUGOSA
## 2 -123.1457    ABIES GRANDIS
## 3 -123.1599      PINUS NIGRA
## 5 -123.1602 PRUNUS SARGENTII
## 6 -123.1605 PRUNUS SARGENTII
## 7 -123.1600      PICEA ABIES</code></pre>
<p>To get taxonomic information, we can take advantage of the rOpenSci package <a rel="nofollow" href="https://github.com/ropensci/taxa" target="_blank"><code>taxa</code></a>. This package provides a function <code>lookup_tax_data</code> which will query NCBI taxonomy. We can then use the package <a rel="nofollow" href="https://grunwaldlab.github.io/metacoder_documentation/index.html" target="_blank"><code>metacoder</code></a> to visualize a “heat tree” for datasets of taxon. This package was originally developed for next generation sequencing and metabarcoding in identifying species from DNA. So the “heat tree” has been diverted from its original purpose but still works very well.</p>
<pre><code># install.packages(&quot;taxa&quot;, &quot;metacoder&quot;)
library(taxa)</code></pre>
<pre><code>## Warning: package 'taxa' was built under R version 3.5.1</code></pre>
<pre><code>library(metacoder)</code></pre>
<pre><code>## Warning: package 'metacoder' was built under R version 3.5.1</code></pre>
<pre><code>## This is metacoder verison 0.3.0 (stable). If you use metacoder for published research, please cite our paper:
## 
## Foster Z, Sharpton T and Grunwald N (2017). &quot;Metacoder: An R package for visualization and manipulation of community taxonomic diversity data.&quot; PLOS Computational Biology, 13(2), pp. 1-15. doi: 10.1371/journal.pcbi.1005404
## 
## Enter `citation(&quot;metacoder&quot;)` for a BibTeX entry for this citation.</code></pre>
<p>To save on space, we will make a table containing each species in our Vancouver dataset and the number of times it occurs.</p>
<pre><code>species_df &lt;- as.data.frame(table(trees_df$SCIENTIFIC_NAME))
names(species_df)&lt;- c(&quot;SCIENTIFIC_NAME&quot;, &quot;count&quot;)
head(species_df)</code></pre>
<pre><code>##      SCIENTIFIC_NAME count
## 1              ABIES    33
## 2     ABIES BALSAMEA     7
## 3     ABIES CONCOLOR     7
## 4      ABIES FRASERI     1
## 5      ABIES GRANDIS    65
## 6 ABIES NORDMANNIANA     1</code></pre>
<p>Querying the NCBI servers for information takes a long time depending on the number of request. If you sign up for a NCBI API key (<a rel="nofollow" href="https://ncbiinsights.ncbi.nlm.nih.gov/2017/11/02/new-api-keys-for-the-e-utilities/" target="_blank">instructions here</a>), and add it to your <code>.Renviron</code> file which is usually located in your home directory.</p>
<p>We are going to use the <code>taxize</code> function <code>classification</code> to obtain the taxonomy of each species. Then, we are inserting the taxonomy for each species back into the <code>trees_df</code> so each tree has an associated taxonomy in column.</p>
<pre><code>library(taxize)

species_classifications &lt;- classification(species_df$SCIENTIFIC_NAME, db = &quot;ncbi&quot;, rows = 1)</code></pre>
<p>This creates a list of data frames containing the taxonomy. We can then extract the taxonomy information by iterating over the classification list and drop all species which do not have taxonomy.</p>
<pre><code>class_to_df &lt;- function(x, ranking){
  if(!is.na(x)){ 
  return(x$name[x$rank == ranking])
    }
}

# species_df &lt;- species_df[species_df$SCIENTIFIC_NAME != &quot;PINUS THUNBERGII&quot;, ]
species_df$superkingdom &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;superkingdom&quot;))
species_df$kingdom &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;kingdom&quot;))
species_df$phylum &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;phylum&quot;))
species_df$subphylum &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;subphylum&quot;))
species_df$subclass &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;subclass&quot;))
species_df$order &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;order&quot;))
species_df$family &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;family&quot;))
species_df$genus &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;genus&quot;))
species_df$species &lt;- as.character(sapply(species_classifications, class_to_df, ranking=&quot;species&quot;))

species_df &lt;- species_df[species_df$superkingdom != &quot;NULL&quot;, ]</code></pre>
<p>We can then join the taxonomy back to the original data set so each tree has its own taxonomy. We also drop the other columns not related to the taxonomy.</p>
<pre><code>trees_df_class &lt;- left_join(trees_df, species_df)
kept_columns &lt;- c('SCIENTIFIC_NAME', 'phylum', 'subphylum', 'order', 'family', 'genus', 'species')
trees_df_class &lt;- trees_df_class[, kept_columns]
## write.csv(trees_df_class, &quot;trees_taxa.csv&quot;, row.names=F)</code></pre>
<p>Next, we use the <code>parse_tax_data</code> function to … … parse the taxonomy columns in the data frame into a <code>tax_map</code> object.</p>
<p>By making a <code>ranks</code> vector, we tell <code>parse_tax_data</code> which columns contain the taxonomy and the ranking order.</p>
<pre><code>ranks &lt;- c('phylum', 'subphylum', 'order', 'family', 'genus', 'species')
obj &lt;- parse_tax_data(trees_df_class, class_cols = ranks, named_by_rank = T)</code></pre>
<pre><code>## Warning: The following 7305 of 124349 input indexes have `NA` in their classifications:
##    29, 33, 37, 38, 64 ... 124328, 124329, 124342, 124347</code></pre>
<p>The resulting object returned from <code>parse_tax_data</code> is actually a special <code>R6</code> object rather than the typical <code>data.frame</code> or S4 object. But since <code>metacoder</code> is based on using <code>taxa</code>, it can directly interface with the <code>R6</code> object produced and we can produce a “heat tree” without any further manipulation.</p>
<pre><code>set.seed(123)

tax_graph &lt;- obj %>% 
  filter_taxa(nchar(taxon_names)!=0) %>%
  filter_taxa(taxon_ranks == &quot;genus&quot;, supertaxa = TRUE) %>%
  heat_tree(node_label = taxon_names, 
            node_size = n_obs, 
            node_label_size = n_obs, 
            node_color = n_obs, 
            node_size_range = c(0.01, .05),  
            node_color_range = RColorBrewer::brewer.pal(5, &quot;BuGn&quot;), 
            node_label_size_range = c(0.02, 0.05),
            node_color_axis_label = &quot;Number of Trees&quot;)

tax_graph</code></pre>
<p><img src="https://i1.wp.com/www.yihanwu.ca/post/2018-11-15-vancouver-trees-part-3_files/figure-html/unnamed-chunk-11-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>The <code>tax_map</code> is filtered to only the genus level in order to produce a legible visualization compared to a species level heat tree.</p>
<p>Node size is related to number of observations. The most common genus in Vancouver are Prunus and Acer.</p>
<pre><code>sessionInfo()</code></pre>
<pre><code>## R version 3.5.0 (2018-04-23)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17134)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_Canada.1252  LC_CTYPE=English_Canada.1252   
## [3] LC_MONETARY=English_Canada.1252 LC_NUMERIC=C                   
## [5] LC_TIME=English_Canada.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] taxize_0.9.3    metacoder_0.3.0 taxa_0.3.1      tidyr_0.8.1    
## [5] purrr_0.2.5     dplyr_0.7.5    
## 
## loaded via a namespace (and not attached):
##  [1] zoo_1.8-2          tidyselect_0.2.4   xfun_0.3          
##  [4] reshape2_1.4.3     lattice_0.20-35    colorspace_1.3-2  
##  [7] htmltools_0.3.6    yaml_2.1.19        rlang_0.2.2       
## [10] pillar_1.2.3       glue_1.2.0         RColorBrewer_1.1-2
## [13] bindrcpp_0.2.2     foreach_1.4.4      bindr_0.1.1       
## [16] plyr_1.8.4         stringr_1.3.1      ggfittext_0.6.0   
## [19] munsell_0.5.0      blogdown_0.8       gtable_0.2.0      
## [22] codetools_0.2-15   evaluate_0.10.1    labeling_0.3      
## [25] knitr_1.20         parallel_3.5.0     curl_3.2          
## [28] Rcpp_0.12.17       backports_1.1.2    scales_0.5.0      
## [31] jsonlite_1.5       ggplot2_3.0.0      digest_0.6.15     
## [34] stringi_1.1.7      bookdown_0.7.17    grid_3.5.0        
## [37] rprojroot_1.3-2    tools_3.5.0        magrittr_1.5      
## [40] lazyeval_0.2.1     tibble_1.4.2       bold_0.5.0        
## [43] crul_0.5.2         crayon_1.3.4       ape_5.1           
## [46] pkgconfig_2.0.1    data.table_1.11.4  xml2_1.2.0        
## [49] assertthat_0.2.0   rmarkdown_1.10     reshape_0.8.7     
## [52] httr_1.3.1         iterators_1.0.9    R6_2.2.2          
## [55] igraph_1.2.1       nlme_3.1-137       compiler_3.5.0</code></pre>
<p>Notes:</p>
<p>The <code>taxa</code> package has a function <code>lookup_tax_data</code> which queries NCBI for and creates a <code>tax_map</code> object. However, I have problems using <code>tree_df</code>, possibly due to the large amount of associated data which has to go into the <code>tax_map</code> object. Additionally, the <code>taxon_ranks</code> are not created directly when using <code>lookup_tax_data</code>. The example code for using <code>lookup_tax_data</code> is below. Also weirdly, <code>lookup_tax_data</code> breaks when querying for <em>Pinus thunbergii</em> the last three times I tried.</p>
<pre><code>trees_df &lt;- trees_df[trees_df$SCIENTIFIC_NAME != &quot;PINUS THUNBERGII&quot;, c(&quot;SCIENTIFIC_NAME&quot;, &quot;GENUS_NAME&quot;)]
tax_env &lt;- lookup_tax_data(trees_df, type = &quot;taxon_name&quot;, column = &quot;SCIENTIFIC_NAME&quot;, database = &quot;ncbi&quot;)</code></pre>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://www.yihanwu.ca/post/vancouver-trees-part-3/"> R on YIHAN WU</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/KMLieZ7XQEE" height="1"></a></div>