R version 4.1.0 (2021-05-18) -- "Camp Pontanezen"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-conda-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.


   Welcome to the Galaxy R Interactive Environment!

The convenience functions gx_put(), gx_get(), and gx_save() are available to you to interact with your current Galaxy history.

  gx_get(42) - Fetch dataset 42 from your Galaxy history. Returns the path where it is available
  gx_put(filename) - Push a dataset to Galaxy
  gx_save() - Save .RHistory, .RData to your Galaxy environment

Some packages are pre-installed for you; RCurl, XML, shiny, ggvis, dplyr, tidyr, ggplot2, reshape2, RODBC, Bioconductor core packages, edgeR, Gviz, Rgraphviz, limma, DESeq2, cummeRbund, affay, Rsamtools 

 You can install new packages with conda/mamba. For example you can do `mamba install bioconductor-tximport`. With the default configuration, all packages from conda-forge and bioconda can be installed.`
> conda install r-ape==5.5
Error: unexpected symbol in "conda install"
> conda install r-ape==5.5
Error: unexpected symbol in "conda install"
> conda install r-ape
Error: unexpected symbol in "conda install"
> library("ape")
Error in library("ape") : there is no package called ‘ape’
> 
> library("ape")
> 
> treefile <- gx_get(###)
+     
+ # Parse the tree
+ tree <- read.tree(treefile)
+ # Plot it
+ plot(tree)
Error: unexpected symbol in:
"# Plot it
plot"
> 
> library("ape")
> 
> treefile <- gx_get(###)
+     
+ # Parse the tree
+ tree <- read.tree(treefile)
+ # Plot it
+ plot(tree)
Error: unexpected symbol in:
"# Plot it
plot"
> 
> #Parse the tree
> tree <- read.tree(treefile)
Error in scan(file = file, what = "", sep = "\n", quiet = TRUE, skip = skip,  : 
  object 'treefile' not found
> #Plot it
> plot(tree)
Error in plot(tree) : object 'tree' not found
> 
> 
> tree <- read.tree(treefile)
Error in scan(file = file, what = "", sep = "\n", quiet = TRUE, skip = skip,  : 
  object 'treefile' not found
> 
> plot(tree)
Error in plot(tree) : object 'tree' not found
> 
> # Parse the tree
> tree <- read.tree(treefile)
Error in scan(file = file, what = "", sep = "\n", quiet = TRUE, skip = skip,  : 
  object 'treefile' not found
> # Plot it
> plot(tree)
Error in plot(tree) : object 'tree' not found
> 
> treefile <- gx_get(119)
DEBUG:root:Host IP determined to be b'172.17.0.1\n'
DEBUG:root:TestURL url=https://usegalaxy.eu obj=False
/opt/miniconda/lib/python3.9/site-packages/bioblend/galaxy/histories/__init__.py:120: FutureWarning: The history_id parameter is deprecated, use the show_history() method to view details of a history for which you know the ID.
  warnings.warn(
DEBUG:urllib3.connectionpool:Starting new HTTPS connection (1): usegalaxy.eu:443
DEBUG:urllib3.connectionpool:https://usegalaxy.eu:443 "GET /api/histories HTTP/1.1" 200 None
DEBUG:root:TestURL url=https://usegalaxy.eu state=success
DEBUG:root:Downloading gx=GalaxyInstance object for Galaxy at https://usegalaxy.eu history=f64524b0a58b291b dataset=119
DEBUG:urllib3.connectionpool:Starting new HTTPS connection (1): usegalaxy.eu:443
DEBUG:urllib3.connectionpool:https://usegalaxy.eu:443 "GET /api/histories/f64524b0a58b291b/contents HTTP/1.1" 200 None
DEBUG:urllib3.connectionpool:Starting new HTTPS connection (1): usegalaxy.eu:443
DEBUG:urllib3.connectionpool:https://usegalaxy.eu:443 "GET /api/datasets/4838ba20a6d867656ef2359c3da9770b?hda_ldda=hda HTTP/1.1" 200 None
DEBUG:urllib3.connectionpool:Starting new HTTPS connection (1): usegalaxy.eu:443
DEBUG:urllib3.connectionpool:https://usegalaxy.eu:443 "GET /api/histories/f64524b0a58b291b/contents/4838ba20a6d867656ef2359c3da9770b/display?to_ext=nhx HTTP/1.1" 200 1109
/import/119
> # Parse the tree
> tree <- read.tree(treefile)
> # Plot it
> plot(tree)
> 
> # Root the tree
> tree_rooted <- root(tree, "ERR313115")
> 
> # Remove the outgroup to make distances within MTB clearer
> tree_rooted <- drop.tip(tree_rooted, "ERR313115")
> tree_rooted$root.edge <- 0.005
> plot(tree_rooted, root.edge = T, cex=0.6)
> 
> # Assign lineages to samples, as identified by TB-profiler
> 
> mtbc_lineages <- c(
+     "ERR181435" = "L7",
+     "ERR313115" = "canettii",
+     "ERR551620" = "L5",
+     "ERR1203059" = "L5",
+     "ERR2659153" = "orygis",
+     "ERR2704678" = "L3",
+     "ERR2704679" = "L1",
+     "ERR2704687" = "L6",
+     "ERR5987300" = "L2",
+     "ERR5987352" = "L4",
+     "ERR6362078" = "L2",
+     "ERR6362138" = "L2",
+     "ERR6362139" = "L4",
+     "ERR6362156" = "L2",
+     "ERR6362253" = "L2",
+     "ERR6362333" = "L2",
+     "ERR6362484" = "L4",
+     "ERR6362653" = "L2",
+     "SRR998584" = "L5",
+     "SRR13046689" = "bovis"
+ )
> 
> # Replace the tree labels
> tree_lineages <- tree_rooted
> tree_lineages$tip.label <- as.character(mtbc_lineages[tree_rooted$tip.label])
> 
> # Define some colors for the lineages
> 
> color_code_lineages = c(
+     L1 = "#ff00ff",
+     L2 = "#0000ff",
+     L3 = "#a000cc",
+     L4 = "#ff0000",
+     L5 = "#663200",
+     L6 = "#00cc33",
+     L7 = "#ede72e",
+     bovis="black",
+     orygis="black")
> 
> pal_lineages <- as.character(color_code_lineages[tree_lineages$tip.label])
> 
> # Plot the old and new tree version next to each other
> par(mfrow = c(1, 2))
> plot(tree_rooted,cex = 0.7, root.edge = TRUE)
> plot(tree_lineages,cex = 0.8, tip.color = pal_lineages, root.edge = TRUE)
> 
> # Same as before, but with DR profiles instead of lineages
> 
> mtbc_dr <- c(
+     "ERR181435" = "Sensitive",
+     "ERR313115" = "Sensitive",
+     "ERR551620" = "MDR",
+     "ERR1203059" = "Sensitive",
+     "ERR2659153" = "Sensitive",
+     "ERR2704678" = "Sensitive",
+     "ERR2704679" = "Sensitive",
+     "ERR2704687" = "Sensitive",
+     "ERR5987300" = "PreXDR",
+     "ERR5987352" = "PreMDR",
+     "ERR6362078" = "MDR",
+     "ERR6362138" = "MDR",
+     "ERR6362139" = "PreMDR",
+     "ERR6362156" = "PreXDR",
+     "ERR6362253" = "MDR",
+     "ERR6362333" = "PreXDR",
+     "ERR6362484" = "PreMDR",
+     "ERR6362653" = "MDR",
+     "SRR998584" = "Sensitive",
+     "SRR13046689" = "Other"
+ )
> 
> tree_dr <- tree_rooted
> tree_dr$tip.label <- as.character(mtbc_dr[tree_rooted$tip.label])
> 
> color_code_dr = c(
+     Sensitive = "#ff00ff",
+     PreXDR = "#0000ff",
+     PreMDR = "#a000cc",
+     MDR = "#ff0000",
+     Other = "#663200"
+ )
> 
> pal_dr <- as.character(color_code_dr[tree_dr$tip.label])
> 
> par(mfrow = c(1, 2))
> plot(tree_rooted,cex = 0.7, root.edge = TRUE)
> plot(tree_dr,cex = 0.8, tip.color = pal_dr, root.edge = TRUE)
> 
> 
> 
> # Rescale branch lengths (here called edge lenghts)
> genome_size = 4411532
> alignment_length = 18077
> invariant_sites = genome_size - alignment_length
> 
> tree_rescaled <- tree_rooted
> tree_rescaled$edge.length <- ((tree_rescaled$edge.length * alignment_length) / genome_size )
> tree_rescaled$root.edge <- ((tree_rescaled$root.edge * alignment_length) / genome_size )
> 
> par(mfrow = c(1, 2))
> plot(tree_rooted,cex = 0.7, root.edge = T, main = "original")
> axisPhylo()
> plot(tree_rescaled,cex = 0.7,root.edge = T, main = "rescaled")
> axisPhylo()
> dev.off()
null device 
          1 
> # Let's also remove the two outlier strains, they would cause troubles and are anyway useless
> tree_rescaled <- drop.tip(tree_rescaled, c("ERR1203059", "ERR5987300"))
> 
> # Estimate dates: translate phylogenetic distance into years by assuming a mutation rate and the number of generations per year
> mutation_rate = 2.01e-10
> generations_per_year = 200
> 
> ## Ape has a function, estimate.dates(), to date a tree by assuming a specific mutation rate
> node.date <- estimate.dates(
+     tree_rescaled,
+     node.dates = rep(0, length(tree_rescaled$tip.label)), # set sampling dates to 0
+     mu = (mutation_rate * generations_per_year) # mutation rate per year
+ )
> 
> tree_rescaled$edge.length <- node.date[tree_rescaled$edge[, 2]] - node.date[tree_rescaled$edge[, 1]]
> tree_rescaled$tip.label <- as.character(mtbc_lineages[tree_rescaled$tip.label])
> 
> plot(tree_rescaled, cex = 0.6, main = "Dated phylogeny (years)")
> axisPhylo()