R Markdown

#This script shows Step 10 of the Unit 4 project
#Download from Canvas and test with example.nwk
#Then change ALL comments and code to read in your tree.nwk file
library(ape) #remember to install package first
## Warning: package 'ape' was built under R version 4.0.3
#Create Newick file of the tree from file downloaded from NCBI Virus
myTree <- read.tree('example.nwk')

#Check that root is the reference genome, "NC_045512.2||China|2019-12"
rerootTree <- root(myTree, myTree$tip.label[9])

#Add formatting to the tree for the different groups
#Create vector to store colors you want for each branch
edgecol <- c()
edgecol[1:20] = "gray" #set all lines to gray
edgecol[7:10]="green" #change color for edge lines connected to B.1.1.7 variants

#Create vector to store colors you want for each label
#Reference is black
#Utah is red
#Spain in purple
tipcol = c("red","purple","purple","red","purple","purple","red","purple","black","red","red")

#Rename tips in rerooted tree
rerootTree$tip.label
##  [1] "MT334526.1_||USA|2020-03-07"   "MT292571.1_||Spain|2020-03-09"
##  [3] "MT292574.1_||Spain|2020-03-02" "MT334533.1_||USA|2020-03-13"  
##  [5] "MT292570.1_||Spain|2020-03-10" "MT292569.1_||Spain|2020-03-09"
##  [7] "MT334529.1_||USA|2020-03-10"   "MT292572.1_||Spain|2020-03-10"
##  [9] "NC_045512.2||China|2019-12"    "MT334523.1_||USA|2020-03-05"  
## [11] "MT334531.1_||USA|2020-03-12"
rerootTree$tip.label <- c('Utah_Mar7','Spain_Mar9','Spain_Mar2',
                      'Utah_Mar13','Spain_Mar10','Spain_Mar9_2',
                      'Utah_Mar10','Spain_Mar10_2','Reference',
                      'Utah_Mar5','Utah_Mar12')

Including Plots

#Use plot.phylo function to create radial style tree
plot.phylo(rerootTree, "radial",
           main="Phylogenetic tree of COVID-19 strains",
           edge.color = edgecol,
           tip.color = tipcol,
           show.tip.label = TRUE, 
           font = .5,
           cex = .5)
legend("bottomleft",
       legend = c("Utah", "Spain", "Reference"), #Name of groups
       fill = c("red","purple","black"),       # Color of the squares
       border = "black", # Color of the border of the squares
       cex =.6, #sets legend size
       xpd= TRUE) #places outside plot area
legend("bottomright",
       legend = c("Matches B.1.1.7", "Not a match"), #Name of groups
       col = c("green","gray"),       # Color of the squares
       border = "black", # Color of the border of the squares
       cex =.6,#sets legend size
       lwd=2, #sets line weight
       xpd=TRUE) #places outside plot area