EBSP PLOT

asli.sa...@gmail.com

unread,

Dec 21, 2018, 10:59:49 AM12/21/18

to beast-users

When I want to make EBSP plot in BEAST 2.4.8. with EBSP log file, following scripts are not working in R STUDİO.

>source("../scripts/plotEBSP.R")

>plotEBSP("EBSP.log", useHPD=FALSE, log="y") in

Do you know any other r scripts for making EBSP plot?

Tim McInerney

unread,

Jan 2, 2019, 9:40:12 PM1/2/19

to beast-users

This is the script I use :)

require(ggplot2)

require(ggthemes)

require(scales)

require(gridExtra)

require(RColorBrewer)

trans.csv = function(infile) {

f = read.table(infile, skip = 1, sep = '\t', header = TRUE)

write.csv(f, infile, row.names = FALSE)

}

trans.log = function(infile, t=1, fix.names = TRUE, remove.na = FALSE, remove.root = TRUE) {

outfile = sub("^([^.]*).*", "\\1", infile) # Outfile name by stripping the extensions and path from the infile

raw.data = read.csv(infile, header = TRUE) # Read in the .log.csv file from the skyline output

if (fix.names==TRUE) {

names(raw.data) = c('Time', 'Mean', 'Median', 'Lower', 'Upper') # Change names of the columns

}

Index = c(1:length(raw.data$Time)) # Create an Index whose length is the number of rows in the data frame

raw.data$Index = Index # Insert the Index as a new column in the raw data file

trans.data = data.frame(Index) # Create a data frame for the transformed data and insert the Index as a new column in the raw data file

for (x in raw.data$Index){

trans.data$Time[x] = raw.data$Time[x]

trans.data$Mean[x] = raw.data$Mean[x] / t

trans.data$Median[x] = raw.data$Median[x] / t

trans.data$Lower[x] = raw.data$Lower[x] / t

trans.data$Upper[x] = raw.data$Upper[x] / t

}

if (remove.root==TRUE) {

trans.data = trans.data[-nrow(trans.data),] # Use this only if you wish to remove the last row (the 'root row') from the data frame

}

if (remove.na==TRUE) {

trans.data = trans.data[complete.cases(trans.data),] # Removes NAs

}

trans.data$Index = NULL # Remove the index column

write.csv(trans.data, paste0(outfile,'_convertedSkyline.csv'), row.names = FALSE)

return(paste0(outfile,'_convertedSkyline.csv'))

}

skyline.plot = function(infile, ext='.png', out.of.africa=FALSE, holocene=FALSE, pleistocene=FALSE, pliocene=FALSE, pop='', species='H.sapiens', skyline.colour="#4DAF4A", save=FALSE) {

outfile = sub("^([^.]*).*", "\\1", infile) # Outfile name by stripping the extensions and path from the infile

Pop.Hist = read.csv(infile, header = TRUE)

species = species

skyline = ggplot(Pop.Hist, aes(x = Time)) +

geom_ribbon(aes(ymin=Lower, ymax=Upper), fill = skyline.colour) +

geom_line(aes(y=Median), size = 1) +

theme_bw() +

{if (holocene==TRUE)geom_vline(xintercept=c(11700), linetype=1, alpha = 0.5, colour = "blue", size = 1) # Holocene marker

} +

{if (pleistocene==TRUE)geom_vline(xintercept=c(2.588E6), linetype=1, alpha = 0.5, colour = "blue", size = 1) # Pleistocene marker

} +

{if (pliocene==TRUE)geom_vline(xintercept=c(5.333E6), linetype=1, alpha = 0.5, colour = "blue", size = 1) # Pleistocene marker

} +

{if (out.of.africa==TRUE)

annotate('text', x = c(37000,53000,81000,100000), y = max(Pop.Hist$Upper) * 1.20, label = c('A','B','C','D'), size = rel(5))

} +

{if (out.of.africa==TRUE)

annotate('rect', xmin = c(29000,47000,73000,94000), xmax = c(45000,59000,89000,106000),ymin = 0, ymax = Inf, alpha = 0.25)

} +

scale_y_log10(labels = comma) +

scale_x_continuous(labels = comma) +

#scale_y_log10(labels = comma, breaks = c(1E2,1E3,1E4,1E5,1E6,1E7,1E8,1E9,1E10)) +

#scale_x_continuous(labels = comma,breaks = seq(0,1E7,1E4)) +

labs(x = 'Years before present',

y = bquote('Effective Population Size ('~N['e']~')'),

#title = expression(paste(.(pop),' (', italic(.(species)), ') Population History'))

title = bquote(.(pop)~italic(.(species))~population~history)

) +

theme(plot.title = element_text(size = rel(1.25), angle = 00, hjust = 0.5)) +

theme(axis.title.y = element_text(size = rel(1.25), angle = 90)) +

theme(axis.title.x = element_text(size = rel(1.25), angle = 00)) +

theme(axis.text.y = element_text(size = rel(1.25), angle = 00)) +

theme(axis.text.x = element_text(size = rel(1), hjust = 1, angle = 45))

if (save==TRUE) {

ggsave(paste0(outfile,'_skyline',ext),skyline, width = 297, height = 210, units = "mm")

}

else {

return(skyline)

}

#################### END TIM FUNCTIONS ####################

message("")

message("ARGUMENT HINTS!")

message("ARGUMENT 1: EBSP .log.csv FILE")

message("ARGUMENT 2: OUT DIRECTORY")

message("ARGUMENT 3: GENERATION TIME")

message("ARGUMENT 4: REMOVE THE COALESCENCE TIME TO ROOT? (T/F)")

message("")

args = commandArgs(trailingOnly = TRUE) # Set arguments from the command line

ebsp.file = args[1] # Input EBSP .log.csv file

out.dir = args[2] # Output directory

gen.time = args[3] # Generation time (default = 25 [Langergraber et al., 2012])

remove.root = args[4] # Remove the coalescence time to root?

if (is.na(out.dir) == TRUE | is.null(out.dir) == TRUE) {

message('OUT DIRECTORY NOT FOUND')

out.dir = dirname(ebsp.file)

message(paste('OUT DIRECTORY ASSIGNED TO ', out.dir))

}

if (endsWith(out.dir, '/') == FALSE) {

out.dir = paste0(out.dir, '/')

}

if (is.na(gen.time) == TRUE | is.null(gen.time) == TRUE) {

message('GENERATION TIME NOT SPECIFIED')

gen.time = 25

message(paste('DEFAULTING TO ', gen.time))

}

if (is.na(remove.root) == TRUE | is.null(remove.root) == TRUE) {

remove.root = TRUE

}

conv.ebsp.file = trans.log(ebsp.file, t = as.numeric(gen.time), fix.names = T, remove.root = remove.root)

skyline.plot(conv.ebsp.file, save = TRUE, holocene = TRUE, out.of.africa = TRUE)

Message has been deleted

asli.sa...@gmail.com

unread,

Jan 4, 2019, 6:11:35 AM1/4/19

to beast-users

Thannk you. I installed required packages. Should we use all scripts begin as trans.csv and TIM FINCTİONs and messages.

3 Ocak 2019 Perşembe 05:40:12 UTC+3 tarihinde Tim McInerney yazdı:

asli.sa...@gmail.com

unread,

Jan 4, 2019, 6:33:18 AM1/4/19

to beast-users

And before using script? how can you turn EBSP log file to csv. file ? Thanks.

4 Ocak 2019 Cuma 14:11:35 UTC+3 tarihinde asli.sa...@gmail.com yazdı:

Tim McInerney

unread,

Jan 4, 2019, 7:36:27 PM1/4/19

to beast-users

Rscript Skyline_Plots.R <csv_output_from_BEAST> <output_directory> <generation_time> <TRUE or FALSE for trimming last data point>

Are you working with humans? Because script also plots a blue line for the holocene and some grey bars for out-of-Africa. You can delete these within one of the functions

Skyline_Plots.R

Tim McInerney

unread,

Jan 4, 2019, 9:28:12 PM1/4/19

to beast-users

Apologies, I've just realised that you stated you need this for BEAST 2, I wrote this for BEAST 1. Hopefully it still works

carlo pacioni

unread,

Jan 6, 2019, 9:00:40 PM1/6/19

to beast-users

Hi asli.sa, can you be more specific, what's the error message you get?

I have found a while back that sometimes the plotEBSP scripts doesn't parse correctly the log file. If that's your problem try the script attached.

carlo

plotEBSP_CP.R

asli.sa...@gmail.com

unread,

Jan 25, 2019, 9:33:52 AM1/25/19

to beast-users

I am new user in R.

Before using script, how can I load EBSP log file in R. When I load it in R studio (import data set from text reading menut) it says is there valid CSV files? And can not view and load files in R consol.