-
Notifications
You must be signed in to change notification settings - Fork 2
/
hierbaps.R
54 lines (40 loc) · 1.66 KB
/
hierbaps.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
## Author: Raymond Kiu [email protected]
## This R script computes phylogenetic clustering using hierbaps method
# Install packages
install.packages("rhierbaps")
install.packages("devtools")
#install.packages("ggtree")
install.packages("phytools")
# install ggtree - special source
source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")
# Load libraries
library(rhierbaps)
library(ggtree)
library(phytools)
# Set seed
set.seed(1234) # To enable reproduction of the same result
# Set working directory (path to your file)
setwd("U:/")
# Load FASTA file
snp.matrix <- load_fasta("coregene_snp.fa",keep.singletons = TRUE)
# Predict clustering (main run)
hb.results <- hierBAPS(snp.matrix, max.depth = 2, n.pops = 20, quiet = TRUE, n.cores = 4)
# n.cores only available in Linux R not in R studio
# Look into your data, first column should be your isolate names
head(hb.results$partition.df)
tail(hb.results$partition.df)
############### this section is additional and not tested ###########################
hb.results <- hierBAPS(snp.matrix, max.depth = 2, n.pops = 20, n.extra.rounds = Inf,
quiet = TRUE)
system.time(hierBAPS(snp.matrix, max.depth = 2, n.pops = 20, quiet = TRUE))
#> user system elapsed
#> 83.378 10.177 96.571
#check running time
# plotting:
newick.file.name <- file("coregene_snp.tree", package = "rhierbaps")
iqtree <- phytools::read.newick(coregene_snp.tree)
#####################################################################################
# To save files in csv:
write.csv(hb.results$partition.df, file = "hierbaps_partition.csv")
#save_lml_logs(hb.results, file.path(tempdir(), "hierbaps_logML.txt"))