run_kdetrees.R

#!/usr/bin/env Rscript

# 7Feb16
# run_kdetrees.R

VERSION <- 'Version: 0.1 - 7Feb16' # runs kdetrees(all.trees.raw, distance = "dissimilarity", topo.only = TRUE &  topo.only = FALSE)
AUTHOR <- "Author: Pablo Vinuesa, CCG-UNAM"
REPOS <- "https://cran.itam.mx/"

# Note: development made on Tenerife@/home/vinuesa/Projects/marfil/PHYLOMARK/Enterobacter_MLSA_primers_Jan16/F2P_primers/dna_amps

#--------------------------------
#>>>>> FUNCTION DEFINITIONS <<<<<
#--------------------------------
print_help <- function(){
  cat("", VERSION, AUTHOR, "",
      ">>> USAGE: ~/R_code/scripts/run_kdetrees.R <src trees extension name> <file_with_multiple_input_newick_trees> <real; outlier detection tuning parameter [def: k=1.5; less is more stringent]",
      "",
      ">>> AIM: runs the non-parametric kdetrees test to find discordant phylogenetic trees using dist=dissimilarity and topo.only = TRUE & topo.only = FALSE",
      "     assumes trees are generated by the multispecies coalescent distribution and detects outlier trees that are unlikely produced by this distribution",
      "",
      ">>> OUTPUT: Returns diagnostic plots of the kernel density estimates (kde), along with a csv file with kde values and test results (ok|outlier), 
            for both topo.only = TRUE and topo.only = FALSE. Also writes summary stat text files",
      "",
      ">>> NOTES:",
      "# 1. ARGUMENTS:",
      " src tree extension names should be different from that of the file holding them all; e.g. <ph> vs <all_trees.tre>",
      "",
      ">>> ASSUMPTIONS ON DATA:",
      " runs kdeobj.diss.topo <- kdetrees(all.trees.raw, distance = \"dissimilarity\", topo.only = TRUE & topo.only = FALSE )",
      " which does not complain if trees are not rooted or have different number of terminals",
      "",
      ">>> REFERENCE: Weyenberg et al. 2014. Bioinformatics: 30(16):2280-2287; PMID:24764459",
      "",
      ">>> TODO:",
      "1. Integrate with compute_suppValStasts_and_RF-dist.R to provide a single dataframe for amplicon quality evaluation!",
      "2. To achieve 1, put this script into a function and run from within the previous script",
      "3. Need to explore the possibilities offered by package distory",
      "4. Add runmodes; runmode==1 could be an evaluation of k=0.75, 1, 1.25 and 1.5",
      "   An additional runmode could be used to get the good alns for phylogenomics (aln->concat->FastTree)",
      "   Alternatively or in addition to 1., this script sould be called from run_core_genome_FT_WAGGphylo.sh",
      "   by parsing the file kde_dfr_file_all_trees.tre.csv",
      "5. refactor code into subrutines for easier calling from other scripts",
      "6. Add getoptLong; see https://cran.r-project.org/web/packages/GetoptLong/index.html",
      "   or simply optparse: https://cran.r-project.org/web/packages/optparse/index.html",
      "",
      sep ="\n")
}
#-----------------------------------------------------------------------

# see http://www.inside-r.org/r-doc/base/file.copy
# for details on file manipulation from R
checkFileCreated <- function(F){
  if( file.exists(F) ){
    message("File ", F, " was created ...")
  }else{
    warning("File ", F,  " could not be written to disk!")
  }
} 
#-----------------------------------------------------------------------
run_kdetrees_diss_topo <- function(tf, k = 1.5)
{
  # compute kernel density estimate of input tree topology distribution  
  # using (distance = "dissimilarity", topo.only = TRUE)
  # trees do not require to be rooted
  #
  # ARGS: 
  #   tf = multinewick file name holding multiple newick strings/trees
  #   k  = the outlier sensintitivity constant: less is more sensitive (detects more outliers); historical default = 1.5
  # Returns:
  #   a kdeobj.diss.topo object
  fun_name <- "run_kdetrees_diss_topo"
  if (missing(tf))
  {
    stop(" Function ", fun_name, " requires a file name holding multiple newick strings/trees")
  }
  
  kdeobj.diss.topo <- kdetrees(tf, distance = "dissimilarity", k=k, topo.only = TRUE )
  return(kdeobj.diss.topo)
}
#-----------------------------------------------------------------------
run_kdetrees_diss_bl <- function(tf, k = 1.5)
{
  # compute kernel density estimate of input tree distribution with branch lengths
  # using (distance = "dissimilarity", topo.only = FALSE)
  # trees do not require to be rooted
  #
  # ARGS: 
  #   tf = multinewick file name holding multiple newick strings/trees
  #   k = the outlier sintitivity constant: less is more sensitive (detects more outliers); historical default = 1.5
  # Returns:
  #   a kdeobj.diss.topo object
  fun_name <- "run_kdetrees_diss_bl"
  if (missing(tf))
  {
    stop(" Function ", fun_name, " requires a file name holding multiple newick strings/trees")
  }
  
  kdeobj.diss.bl <- kdetrees(tf, distance = "dissimilarity", k=k, topo.only = FALSE )
  return(kdeobj.diss.bl)
}
#-----------------------------------------------------------------------
run_kdetrees_k_check <- function(tf, k_vec = c(0.75, 1, 1.25, 1.5))
{
  # get the kdeobj using distance = "dissimilarity", topo.only = TRUE; trees do not require to be rooted
  #kdeobj <- kdetrees(all.trees.raw, distance = "dissimilarity", topo.only = TRUE )
  fun_name <- "run_kdetrees_k_check"
  if (missing(tf))
  {
    stop(" Function ", fun_name, " requires a file name holding multiple newick strings/trees")
  }
  
  for (i in k_vec)
  {
    kde_obj <- paste("kdeobj.diss.topo.", i, sep="")
    kde_obj <- run_kdetrees_diss_topo(trees_file, i)
    
    
    return(kdeobj.diss.topo)
  }  
}
#-----------------------------------------------------------------------

#####################
##### MAIN CODE #####
#####################

#-------------------------
#>>>>> GET USER ARGS <<<<<
#-------------------------
argv <- commandArgs(TRUE)

if(length(argv) < 2)
{
  print_help()
  stop(" Usage: <src trees extension names> <file_with_multiple_input_newick_trees> <real; outlier detection tuning parameter [def: k=1.5; less is more stringent]")
}

tree_ext <- as.character(argv[1])
input_trees_file <- as.character(argv[2])
if(length(argv) == 2)
{
  k.in <- 1.5
}else{ k.in <- as.numeric(argv[3]) }

message("#>>> Running with arguments: ", tree_ext, " ", input_trees_file, " k=", k.in, " ...")

#--------------------------
#>>>>> LOAD LIBRARIES <<<<<
#--------------------------
# see ?kdetrees for more info on the function
# see help(package=kdetrees)
library("stringr")
library("ape")
library("kdetrees")
library("vioplot")

# initialize vars
no.tips.vec <- c()
no.tips.dfr <- c()
col.vec.topo <- c()
flag.vec.topo <-c()
flag.vec.bl <-c()
combined.flag.vec.bl.topo <- c()
kde_bl_topo_test <-c()

# 1. get the list of tree files
# >>> pass the tree_ext arg to list.files() funct as a regex
rgx <- paste("\\.", tree_ext, "$", sep= "")
files <- list.files(pattern=rgx)

# check there are tree files with tree_ext extension in the working directory
if(length(files) == 0) stop("There are no tree files with ", tree_ext, " extension in the working directory! Will stop now ...")

# save old par() to reset after manipulating it for the plots
opar <- par(no.readonly = TRUE)

# Need to find out how to create a multiphylo object by directly reading in multiple trees
# like from files; may need tree.names = files
#for (i in files){
#  all.trees.raw <- read.tree(file=files[i], tree.names = files, keep.multi = TRUE )
#}

# read all trees from a concatenated file holding them all
# system("cat *.ph > all_trees.raw.tre")
all.trees.raw <- read.tree(file=input_trees_file)

# get the kdeobj using distance = "dissimilarity", topo.only = TRUE; trees do not require to be rooted
#kdeobj <- kdetrees(all.trees.raw, distance = "dissimilarity", topo.only = TRUE )

kdeobj.diss.bl <-run_kdetrees_diss_bl(all.trees.raw, k.in)
kdeobj.diss.topo <-run_kdetrees_diss_topo(all.trees.raw, k.in)
#kdeobj.diss.topo <- kdetrees(all.trees.raw, distance = "dissimilarity", k=k.in, topo.only = TRUE )
#kdeobj.diss.bl <- kdetrees(all.trees.raw, distance = "dissimilarity", k=k.in, topo.only = FALSE )

# print overview stats
kde_stats_file <- paste("kde_stats_", input_trees_file, ".out", sep = "")
sink(file = kde_stats_file, type=c("output"))
kdeobj.diss.bl
kdeobj.diss.topo
sink()

checkFileCreated(kde_stats_file)

# These lines, using dist=geodesic do not work; require rooting
#kdeobj.geod.bl <- kdetrees(all.trees.raw, distance = "geodesic", outgroup= "Cronobacter_sakazakii_ATCC_BAA-894", topo.only = FALSE )
#kdeobj.geod.topo <- kdetrees(all.trees.raw, distance = "geodesic", outgroup="Cronobacter_sakazakii_ATCC_BAA-894", topo.only = TRUE )
# print overview stats
#kdeobj.geod.bl
#kdeobj.geod.topo

# make parallel boxplots 
svg(file="parallel_bxplots_kdeDensity_dist_dissim_topo_TRUE-FALSE.svg")
layout(matrix( c(1,2), 1, 2, byrow = TRUE) )
boxplot(kdeobj.diss.topo$density, main="dist=dissim., topo.only=T")
boxplot(kdeobj.diss.bl$density, main="dist=dissim., topo.only=F")
dev.off()
par(opar)

checkFileCreated("parallel_bxplots_kdeDensity_dist_dissim_topo_TRUE-FALSE.svg")

# print the bad files to screen
topo.outlier.tree.idx <- kdeobj.diss.topo$i
message("there are ", length(topo.outlier.tree.idx), " outlier trees")

bl.outlier.tree.idx <- kdeobj.diss.bl$i
message("there are ", length(bl.outlier.tree.idx), " outlier trees")

# create a color vector, to plot the good tree points in blue and outliers in black
src.tree.idx = 1:length(all.trees.raw)
col.vec.topo <- src.tree.idx %in% topo.outlier.tree.idx
col.vec.topo <- ifelse(col.vec.topo, col.vec.topo <- c("black"), col.vec.topo <- c("blue"))

col.vec.bl <- src.tree.idx %in% bl.outlier.tree.idx
col.vec.bl <- ifelse(col.vec.bl, col.vec.bl <- c("black"), col.vec.bl <- c("blue"))

# create plot of kde density points and a boxplot summarithing their distribution
# Note the use of fig= for fine control of placement
svg(file="dotplot_and_bxplot_kdeDensity_dist_dissim_topo_TRUE.svg")
#layout(matrix( c(1,2), 1, 2, byrow = TRUE) )
main_txt <- paste("dist=dissim., topo.only=T, k=", k.in, sep = "")
par(fig=c(0, 0.9, 0, 0.5))
plot(kdeobj.diss.topo$density, col = col.vec.topo, main = main_txt )

# notice the use of new=T; otherwise the boxplot would wipe out the 1st plot!

par(fig=c(0.75, 1, 0, 0.5), new=T) 
boxplot(kdeobj.diss.topo$density, col = col.vec.topo, axes=FALSE, main="k=1.5")

main_txt <- paste("dist=dissim., topo.only=F, k=", k.in, sep = "")
par(fig=c(0, 0.9, 0.5, 1), new = T)
plot(kdeobj.diss.bl$density, col = col.vec.bl, main = main_txt )

par(fig=c(0.75, 1, 0.5, 1), new=T) 
boxplot(kdeobj.diss.bl$density, col = col.vec.bl, axes=FALSE, main="k=1.5")

dev.off()
par(opar)

checkFileCreated("dotplot_and_bxplot_kdeDensity_dist_dissim_topo_TRUE.svg")

# print the outlier files to 
kde_outlier_files <- paste("kde_outlier_files_", input_trees_file, ".out", sep = "")
sink( file = kde_outlier_files, type=c("output") )
message("# These are the outlier trees using dist=dissimilarity and topo.only=T, with k = ", k.in)
files[topo.outlier.tree.idx]
message("===============================================================================================")
message("# These are the outlier trees using dist=dissimilarity and topo.only=F, with k = ", k.in)
files[bl.outlier.tree.idx]
sink()

checkFileCreated(kde_outlier_files)

# construct a dataframe
flag.vec.topo <- src.tree.idx %in% topo.outlier.tree.idx
flag.vec.topo <- ifelse(flag.vec.topo, flag.vec.topo <- c("outlier"), flag.vec.topo <- c("ok")) 

flag.vec.bl <- src.tree.idx %in% bl.outlier.tree.idx
flag.vec.bl <- ifelse(flag.vec.bl, flag.vec.bl <- c("outlier"), flag.vec.bl <- c("ok")) 

#flag.vec.topo.comb <- src.tree.idx %in% topo.outlier.tree.idx
#flag.vec.bl.comb <- src.tree.idx %in% bl.outlier.tree.idx
#combined.flag.vec.bl.topo <- unique(union(flag.vec.topo.comb, flag.vec.bl.comb))
#combined.flag.vec.bl.topo <- ifelse(combined.flag.vec.bl.topo, combined.flag.vec.bl.topo <- c("outlier"), combined.flag.vec.bl.topo <- c("ok")) 

kde.densities.vec.topo <- kdeobj.diss.topo$density
kde.densities.vec.bl <- kdeobj.diss.bl$density

#kde.trees.dfr <- data.frame(file=files, kde_topo_dens=kde.densities.vec.topo, kde_topo_test=flag.vec.topo, kde_bl_dens=kde.densities.vec.bl, kde_bl_test=flag.vec.bl, kde_bl_topo_test=combined.flag.vec.bl.topo)
kde.trees.dfr <- data.frame(file=files, kde_topo_dens=kde.densities.vec.topo, kde_topo_test=flag.vec.topo, kde_bl_dens=kde.densities.vec.bl, kde_bl_test=flag.vec.bl)

# Genarate a new variable kde_bl_topo_test, that combines the outilers found by both kde_topo_test & kde_bl_test 
kde.trees.dfr <- within(kde.trees.dfr, {
  kde_bl_topo_test[kde_topo_test == "ok" & kde_bl_test == "ok"] <- "ok"
  kde_bl_topo_test[kde_topo_test == "outlier"] <- "outlier"
  kde_bl_topo_test[kde_bl_test == "outlier"] <- "outlier"
})

# write kde.trees.dfr to file:
kde_dfr_file <- paste("kde_dfr_file_", input_trees_file, ".tab", sep = "")
write.table(kde.trees.dfr, file=kde_dfr_file, row.names = FALSE, sep = "\t", quote = FALSE)
checkFileCreated(kde_dfr_file)

# make histograms to summarize distributions of tree KDEs
# 1. get colors
# http://stackoverflow.com/questions/21858394/partially-color-histogram-in-r
# Here's the method I mentioned in comments:
# Make some test data (you should do this in your question!)
# test = runif(10000,-2,0)
# get R to compute the histogram but not plot it:
# h = hist(test, breaks=100,plot=FALSE)
# Your histogram is divided into three parts:
# ccat = cut(h$breaks, c(-Inf, -0.6, -0.4, Inf))
# plot with this palette, implicit conversion of factor to number indexes the palette:
# plot(h, col=c("white","green","red")[ccat])

topo.cutoff <- max(kde.trees.dfr$kde_topo_dens[kde.trees.dfr$kde_topo_test == "outlier"])
h.topo <- hist(kde.trees.dfr$kde_topo_dens, breaks = 50, plot=FALSE) # probability = TRUE, <== does not like it!?
ccat.topo <- cut(h.topo$breaks, c(-Inf, topo.cutoff, Inf) )

bl.cutoff <- max(kde.trees.dfr$kde_bl_dens[kde.trees.dfr$kde_bl_test == "outlier"])
h.bl <- hist(kde.trees.dfr$kde_bl_dens, breaks = 50, plot=FALSE) # probability = TRUE, <== does not like it!?
ccat.bl <- cut(h.bl$breaks, c(-Inf, bl.cutoff, Inf) )

hist_plot_file <- paste("kde_hist_plot_file_", input_trees_file, ".svg", sep = "")
svg(file=hist_plot_file)
layout(matrix( c(1,2), 2, 1, byrow = TRUE) )
plot(h.topo, main="kde topo", col=c("black", "blue")[ccat.topo])
rug(jitter(kde.trees.dfr$kde_topo_dens))
#lines(density(kde.trees.dfr$kde_topo_dens), lwd=2)
plot(h.bl, main="kde bl", col=c("black", "blue")[ccat.bl])
rug(jitter(kde.trees.dfr$kde_bl_dens))
#lines(density(kde.trees.dfr$kde_bl_dens), lwd=2)
dev.off()

checkFileCreated(hist_plot_file)

# make violin plots to summarize distributions of tree KDEs
violin_plot_file <- paste("violin_plot_file_", input_trees_file, ".svg", sep = "")
svg(file=violin_plot_file)
layout(matrix( c(1,2), 2, 1, byrow = TRUE) )
vioplot( kde.trees.dfr$kde_bl_dens, col="gold", names = c("kde for tree distributions with branch lengths") ) 
vioplot( kde.trees.dfr$kde_topo_dens, col="gold", names = c("kde for topology distributions") ) 
dev.off()

checkFileCreated(violin_plot_file)

# Try doing a consensus;
# Error in FUN(X[[i]], ...) : one tree has a different number of tips
# consensus(good.trees, p = 0.5)


# exit without saving workspace
q(save = "no")




##### THIS SNIPPET WORKS, USING ONLY ROOTED TREES
# however, many trees are skipped; need to revise;
# for the moment, use the code/script above



# remove previous ones
system("rm rooted_*ph all_trees_rooted.tre")

# 1 get the no of tips of each src tree 
#   and construct a dfr with the file_name and no.tips
files <- list.files(pattern="\\.ph$")
for ( i in 1:length(files) ){ 
   tr <- read.tree(files[i])
   no.tips <- length(tr$tip.label)
   no.tips.vec <- c(no.tips.vec, no.tips)
   
   #message(files[i], "\t", no.tips)
   # construct a dataframe with "file", "no_tips"
   no.tips.dfr <- rbind( no.tips.dfr, data.frame(files[i], no.tips))
}
names(no.tips.dfr) <- c("file", "no_tips")

# hist the (no.tips.vec and get the max.no.tips
str(no.tips.dfr)
hist(no.tips.vec)

max.no.tips <- max(no.tips.dfr$no_tips)

# subset the dfr no.tips.dfr, to get the files with no_tips == max.no.tips
good.files.dfr <- no.tips.dfr[no.tips.dfr$no_tips == max.no.tips, ]
#good.files <- good.files.dfr$file
good.files <- as.character(good.files.dfr$file)


# read those trees
for (i in 1:length(good.files) ){ 
      tr <- read.tree(good.files[i])
      no.tips <- length(tr$tip.label)
      
     # if ( i == 1 & no.tips == max.no.tips) { t.labs <- tr$tip.label }
     #kdetrees requires that all src trees to be compared have the same no. of leaves 
     # ant that the tree labels are equal
     #if (no.tips == max(no.tips.vec) & tr$tip.label == t.labs ){
     #out.check <- grep(outgr, tr) 
     #if(out.check == outgr){
          rtr <-  root(tr, "Cronobacter_sakazakii_ATCC_BAA-894", resolve.root = TRUE)
          rtr.file <- paste("rooted_", good.files[i], sep ="")
          if( is.rooted(rtr)){
             message("writing file", rtr.file, " ...")
              write.tree(rtr, file=rtr.file)
          }
          else{ warning("could not root tree", rtr.file)}
          #else { next }
       #}
       #else { next }
   #}
}



# files <- list.files(pattern="rooted_*\\.ph$")
# for (i in 1:length(files) ){ 
#   tr <- read.tree(files[i])
#   no.tips <- length(tr$tip.label)
#   no.tips.vec <- c(no.tips.vec, no.tips)
# }
# hist(no.tips.vec)

system("cat rooted_*ph > all_trees_rooted.tre")

#Error in dist.multiPhylo(trees, ...) : Not all trees have the same tips.
# they all have 28 tips, but no all are equal!?



good.rooted.trees <- read.tree("all_trees_rooted.tre")
#kdetrees.complete("all_trees_rooted.tre")
kdetrees(good.rooted.trees, distance = "dissimilarity", topo.only = TRUE )
kdeobj <- kdetrees(good.rooted.trees, distance = "dissimilarity", topo.only = TRUE )
topo.outlier.tree.idx <- print(kdeobj$i)
plot(kdeobj$density)
boxplot(kdeobj$density)