--- title: "Elzar Benchmarking Results" author: "RA" output: rmdformats::downcute: self_contained: true thumbnails: false lightbox: true #gallery: false #fig_caption: true #highlight: kate downcute_theme: "chaos" editor_options: chunk_output_type: console --- ## 1.0 Loading Required Libraries ```{r, warning=FALSE, message=FALSE} library(splitstackshape) library(dplyr) library(ggplot2) library(VennDiagram) library(RColorBrewer) library(gridExtra) library(cowplot) library(knitr) library(kableExtra) library(tidyverse) library(Hmisc) library(ggridges) library(venn) library(ggthemes) library(stringr) library(splitstackshape) library(chron) library(magicfor) library(ggtext) library(reshape2) library(plyr) library(formattable) library(data.table) library(ggrepel) library(chron) library(lubridate) require(MASS) require(scales) ``` ```{r} #Setting working directory setwd("/home/asiimwe/projects/run_env/alpha_star_wasp_benchmarking/dataExtractions/") ``` ## Data Loading ```{r} #https://man7.org/linux/man-pages/man2/getrusage.2.html setwd("/home/asiimwe/projects/run_env/alpha_star_wasp_benchmarking/") elzar_runs <- read.table("dataExtractions/benchmarking_stats_elzar.txt",comment.char = "", header = FALSE,fill=TRUE) ``` ## Data cleaning ```{r} elzar_runs$V2 <- paste0(elzar_runs$V2, " ", elzar_runs$V3) elzar_runs$V3 <- NULL elzar_runs_transposed <- as.data.frame(t(elzar_runs)) colnames(elzar_runs_transposed) <- unlist(elzar_runs_transposed[c(1),]) elzar_runs_transposed <- elzar_runs_transposed[-c(1),] k <- 8 nr <- nrow(elzar_runs_transposed) nc <- ncol(elzar_runs_transposed) unames <- unique(names(elzar_runs_transposed)) a <- array(as.matrix(elzar_runs_transposed), c(nr, k, nc/k)) m <- matrix(aperm(a, c(1, 3, 2)),, k, dimnames = list(NULL, unames)) elzar_runs_transposed <- as.data.frame(m, stringsAsFactors = FALSE) #duplicated(elzar_runs_transposed[,1]) #need to filter given some runs were repeated on the same day - filter by time head(elzar_runs_transposed$start_time) elzar_runs_transposed <- cSplit(elzar_runs_transposed, "start_time", " ", direction = "wide", type.convert = "character") elzar_runs_transposed <- cSplit(elzar_runs_transposed, "end_time", " ", direction = "wide", type.convert = "character") #elzar_runs_transposed %>% filter(!duplicated(elzar_runs_transposed$jobname)) # #Capture sample, threads, wall clock and memory # #Wall clock is reported in seconds # time <- "22:13:59.760" # time2 <- "22:23:11.209" # x <- chron(times=time) # y <- chron(times=time2) # y-x # period_to_seconds(hms(y-x)) #551 exact reported wall clock is 551.510 elzar_runs_transposed$cwd_base <- elzar_runs_transposed$cwd elzar_runs_transposed <- cSplit(elzar_runs_transposed, "cwd", "/", direction = "wide", type.convert="as.character") #Filtering to remove prior runs unique(elzar_runs_transposed$cwd_07) runs <- c("STAR_Run1", "STAR_Run2", "STAR_Run3", "STAR_WASP_Run1", "STAR_WASP_Run2", "STAR_WASP_Run3", "WASP_Run1", "WASP_Run2", "WASP_Run3") all(elzar_runs_transposed$cwd_07 %in% runs) #elzar_runs_transposed <- elzar_runs_transposed %>% filter(cwd_07 %in% runs) - prior calls had diff run labels elzar_runs_transposed[,c(11:16)] <- NULL dim(elzar_runs_transposed) dim(elzar_runs_transposed <- na.omit(elzar_runs_transposed)) colnames(elzar_runs_transposed)[6:14] <- c("start_date", "start_time", "end_date", "end_time", "Path", "Run_Group", "Run", "Sample", "Threads") elzar_runs_transposed$wallclock <- as.numeric(as.character(elzar_runs_transposed$wallclock)) elzar_runs_transposed$ru_wallclock <- as.numeric(as.character(elzar_runs_transposed$ru_wallclock)) elzar_runs_transposed$maxrss <- gsub("G", "", elzar_runs_transposed$maxrss) elzar_runs_transposed$maxrss <- as.numeric(as.character(elzar_runs_transposed$maxrss)) elzar_runs_transposed$ru_maxrss <- as.numeric(as.character(elzar_runs_transposed$ru_maxrss)) elzar_runs_transposed$Run <- gsub("_Runs", "", elzar_runs_transposed$Run) elzar_runs_transposed$Threads <- ordered(elzar_runs_transposed$Threads , levels = c("8threads", "16threads", "32threads")) elzar_runs_transposed <- elzar_runs_transposed %>% filter(Sample != "HG00514" & Sample != "NA12878_Small" & Sample != "NA12878_RAMPAGE" & Sample != "NA12878_RAMPAGE_Rep") unique(elzar_runs_transposed$Sample) #Coverting date and time fields from character to data and time respectively # elzar_runs_transposed$start_date <- chron(date=elzar_runs_transposed$start_date) # elzar_runs_transposed$start_time <- chron(times=elzar_runs_transposed$start_time) # elzar_runs_transposed$end_date <- chron(date=elzar_runs_transposed$end_date) # elzar_runs_transposed$end_time <- chron(times=elzar_runs_transposed$end_time) unique(elzar_runs_transposed$Run) dim(STAR <- elzar_runs_transposed %>% filter(Run == "STAR")) #208 we expect 16 * 3 = 48 dim(STAR_WASP <- elzar_runs_transposed %>% filter(Run == "STAR_WASP")) #203 we expect 16 * 3 = 48 dim(WASP <- elzar_runs_transposed %>% filter(Run == "WASP")) #208 we expect 16 * 3 = 48 #Checking call groups: dim(STAR_Run1 <- as.data.frame(STAR %>% filter(Run_Group == "STAR_Run1"))) #48 STAR_Run1 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) #checking that for each sample we have 3 runs, 1 for each thread - looks good STAR_Run1 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #looks good dim(STAR_Run2 <- STAR %>% filter(Run_Group == "STAR_Run2"))#96 STAR_Run2 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) #We have twice the expected number of runs, data called might include a prior run that was conducted on the same day or prior, we shall filter these by date STAR_Run2 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #duplicated dim(STAR_Run3 <- STAR %>% filter(Run_Group == "STAR_Run3"))#64 STAR_Run3 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) #We also have an extra run herein, we shall filter by date as well STAR_Run3 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #some are duplicated #STAR_WASP dim(STAR_WASP_Run1 <- as.data.frame(STAR_WASP %>% filter(Run_Group == "STAR_WASP_Run1"))) #48 STAR_WASP_Run1 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) STAR_WASP_Run1 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n())#looks good dim(STAR_WASP_Run2 <- STAR_WASP %>% filter(Run_Group == "STAR_WASP_Run2")) #96 STAR_WASP_Run2 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) STAR_WASP_Run2 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #we have duplicates dim(STAR_WASP_Run3 <- STAR_WASP %>% filter(Run_Group == "STAR_WASP_Run3"))#59 STAR_WASP_Run3 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) STAR_WASP_Run3 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #we have duplicates #WASP dim(WASP_Run1 <- as.data.frame(WASP %>% filter(Run_Group == "WASP_Run1"))) #48 WASP_Run1 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) WASP_Run1 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n())#looks good dim(WASP_Run2 <- WASP %>% filter(Run_Group == "WASP_Run2")) #71 WASP_Run2 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) WASP_Run2 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #we have duplicates dim(WASP_Run3 <- WASP %>% filter(Run_Group == "WASP_Run3"))#89 WASP_Run3 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) WASP_Run3 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #we have duplicates # Filtering out double calls dim(STAR_Run2 <- as.data.frame(STAR_Run2 %>% filter(start_date >= as.Date("05/05/2022")))) STAR_Run2 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) STAR_Run2 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) dim(STAR_Run3 <- as.data.frame(STAR_Run3 %>% filter(start_date >= as.Date("05/05/2022")))) STAR_Run3 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) STAR_Run3 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) #STAR_WASP dim(STAR_WASP_Run2 <- as.data.frame(STAR_WASP_Run2 %>% filter(start_date >= as.Date("05/05/2022")))) STAR_WASP_Run2 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) STAR_WASP_Run2 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) dim(STAR_WASP_Run3 <- as.data.frame(STAR_WASP_Run3 %>% filter(start_date >= as.Date("05/05/2022")))) STAR_WASP_Run3 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) STAR_WASP_Run3 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) # WASP dim(WASP_Run2 <- as.data.frame(WASP_Run2 %>% filter(start_date >= as.Date("05/05/2022")))) WASP_Run2 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) WASP_Run2 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) dim(WASP_Run3 <- as.data.frame(WASP_Run3 %>% filter(start_date >= as.Date("05/05/2022")))) WASP_Run3 %>% group_by(Sample) %>% dplyr::summarise(cnt = n()) WASP_Run3 %>% group_by(Sample, Threads) %>% dplyr::summarise(cnt = n()) dfs <- list(STAR_Run1, STAR_Run2, STAR_Run3, STAR_WASP_Run1, STAR_WASP_Run2, STAR_WASP_Run3, WASP_Run1, WASP_Run2, WASP_Run3) for(i in dfs){print(class(i))} for(i in dfs){print(nrow(i))} #merging dfs #STAR duplicated(STAR_Run1$jobname) duplicated(STAR_Run2$jobname) duplicated(STAR_Run3$jobname) all(STAR_Run1$jobname %in% STAR_Run2$jobname) all(STAR_Run2$jobname %in% STAR_Run3$jobname) colnames(STAR_Run1)[2:ncol(STAR_Run1)] <- paste0(colnames(STAR_Run1)[2:ncol(STAR_Run1)], "_Run1") colnames(STAR_Run2)[2:ncol(STAR_Run2)] <- paste0(colnames(STAR_Run2)[2:ncol(STAR_Run2)], "_Run2") colnames(STAR_Run3)[2:ncol(STAR_Run3)] <- paste0(colnames(STAR_Run3)[2:ncol(STAR_Run3)], "_Run3") dim(STAR_merged <- inner_join(STAR_Run1, STAR_Run2, by = "jobname")) dim(STAR_merged <- inner_join(STAR_merged, STAR_Run3, by = "jobname")) #STAR_WASP duplicated(STAR_WASP_Run1$jobname) duplicated(STAR_WASP_Run2$jobname) duplicated(STAR_WASP_Run3$jobname) all(STAR_WASP_Run1$jobname %in% STAR_WASP_Run2$jobname) all(STAR_WASP_Run2$jobname %in% STAR_WASP_Run3$jobname) colnames(STAR_WASP_Run1)[2:ncol(STAR_WASP_Run1)] <- paste0(colnames(STAR_WASP_Run1)[2:ncol(STAR_WASP_Run1)], "_Run1") colnames(STAR_WASP_Run2)[2:ncol(STAR_WASP_Run2)] <- paste0(colnames(STAR_WASP_Run2)[2:ncol(STAR_WASP_Run2)], "_Run2") colnames(STAR_WASP_Run3)[2:ncol(STAR_WASP_Run3)] <- paste0(colnames(STAR_WASP_Run3)[2:ncol(STAR_WASP_Run3)], "_Run3") dim(STAR_WASP_merged <- inner_join(STAR_WASP_Run1, STAR_WASP_Run2, by = "jobname")) dim(STAR_WASP_merged <- inner_join(STAR_WASP_merged, STAR_WASP_Run3, by = "jobname")) #WASP duplicated(WASP_Run1$jobname) duplicated(WASP_Run2$jobname) duplicated(WASP_Run3$jobname) all(WASP_Run1$jobname %in% WASP_Run2$jobname) all(WASP_Run2$jobname %in% WASP_Run3$jobname) colnames(WASP_Run1)[2:ncol(WASP_Run1)] <- paste0(colnames(WASP_Run1)[2:ncol(WASP_Run1)], "_Run1") colnames(WASP_Run2)[2:ncol(WASP_Run2)] <- paste0(colnames(WASP_Run2)[2:ncol(WASP_Run2)], "_Run2") colnames(WASP_Run3)[2:ncol(WASP_Run3)] <- paste0(colnames(WASP_Run3)[2:ncol(WASP_Run3)], "_Run3") dim(WASP_merged <- inner_join(WASP_Run1, WASP_Run2, by = "jobname")) dim(WASP_merged <- inner_join(WASP_merged, WASP_Run3, by = "jobname")) # Averaging memory and speed across runs (ru_wallclock_Run1 ru_maxrss_Run1 wallclock_Run1 maxrss_Run1) #STAR STAR_merged <- STAR_merged %>% mutate(average_ru_wallclock = rowMeans(STAR_merged[,c("ru_wallclock_Run1", "ru_wallclock_Run2", "ru_wallclock_Run3")])) #STAR_merged <- STAR_merged %>% mutate(average_ru_wallclock2 = (ru_wallclock_Run1 + ru_wallclock_Run2 + ru_wallclock_Run3)/3) #all(round(STAR_merged$average_ru_wallclock,2) == round(STAR_merged$average_ru_wallclock2,2)) #TRUE STAR_merged <- STAR_merged %>% mutate(average_ru_maxrss = rowMeans(STAR_merged[,c("ru_maxrss_Run1","ru_maxrss_Run2","ru_maxrss_Run3")])) STAR_merged <- STAR_merged %>% mutate(average_wallclock = rowMeans(STAR_merged[,c("wallclock_Run1","wallclock_Run2","wallclock_Run3")])) STAR_merged <- STAR_merged %>% mutate(average_maxrss = rowMeans(STAR_merged[,c("maxrss_Run1","maxrss_Run2","maxrss_Run3")])) #STAR_WASP STAR_WASP_merged <- STAR_WASP_merged %>% mutate(average_ru_wallclock = rowMeans(STAR_WASP_merged[,c("ru_wallclock_Run1", "ru_wallclock_Run2", "ru_wallclock_Run3")])) STAR_WASP_merged <- STAR_WASP_merged %>% mutate(average_ru_maxrss = rowMeans(STAR_WASP_merged[,c("ru_maxrss_Run1","ru_maxrss_Run2","ru_maxrss_Run3")])) STAR_WASP_merged <- STAR_WASP_merged %>% mutate(average_wallclock = rowMeans(STAR_WASP_merged[,c("wallclock_Run1","wallclock_Run2","wallclock_Run3")])) STAR_WASP_merged <- STAR_WASP_merged %>% mutate(average_maxrss = rowMeans(STAR_WASP_merged[,c("maxrss_Run1","maxrss_Run2","maxrss_Run3")])) #WASP WASP_merged <- WASP_merged %>% mutate(average_ru_wallclock = rowMeans(WASP_merged[,c("ru_wallclock_Run1", "ru_wallclock_Run2", "ru_wallclock_Run3")])) WASP_merged <- WASP_merged %>% mutate(average_ru_maxrss = rowMeans(WASP_merged[,c("ru_maxrss_Run1","ru_maxrss_Run2","ru_maxrss_Run3")])) WASP_merged <- WASP_merged %>% mutate(average_wallclock = rowMeans(WASP_merged[,c("wallclock_Run1","wallclock_Run2","wallclock_Run3")])) WASP_merged <- WASP_merged %>% mutate(average_maxrss = rowMeans(WASP_merged[,c("maxrss_Run1","maxrss_Run2","maxrss_Run3")])) # Merging datasets for plotting STAR_merged_subset <- STAR_merged[,c(38:ncol(STAR_merged))] colnames(STAR_merged_subset)[1:3] <- c("Run", "Sample", "Threads") STAR_WASP_merged_subset <- STAR_WASP_merged[,c(38:ncol(STAR_WASP_merged))] colnames(STAR_WASP_merged_subset)[1:3] <- c("Run", "Sample", "Threads") WASP_merged_subset <- WASP_merged[,c(38:ncol(WASP_merged))] colnames(WASP_merged_subset)[1:3] <- c("Run", "Sample", "Threads") colnames(STAR_merged_subset) colnames(STAR_WASP_merged_subset) colnames(WASP_merged_subset) merged_df <- as.data.frame(rbind(STAR_merged_subset, STAR_WASP_merged_subset, WASP_merged_subset)) ``` ### Average Memory Across All Runs ```{r} ## Cleaning up Sample names unique(merged_df$Sample) #Removing NA12878 from all NA12878 derived samples merged_df$Sample <- gsub("NA12878_", "", merged_df$Sample) merged_df$Sample <- gsub("Rep", "R2", merged_df$Sample) merged_df$Sample[merged_df$Sample == "Total"] <- "Total_R1" merged_df$Sample[merged_df$Sample == "PolyA"] <- "PolyA_R1" merged_df$Sample[merged_df$Sample == "Nucleus_PolyA"] <- "Nucleus_PolyA_R1" merged_df$Sample[merged_df$Sample == "Nucleus_nonPolyA"] <- "Nucleus_nonPolyA_R1" order_df <- c("Total_R2" , "Total_R1" , "PolyA_R2" , "PolyA_R1" , "Nucleus_PolyA_R2", "Nucleus_PolyA_R1" , "Nucleus_nonPolyA_R2", "Nucleus_nonPolyA_R1", "NA19240" , "NA19239" , "NA19238" , "HG00733" , "HG00732" , "HG00731", "HG00513" , "HG00512") global_colors <- c("orange", "dodgerblue4", "firebrick4") merged_df %>% ggplot(aes(x = factor(Sample, levels=order_df), y = average_maxrss, group=Run,color=factor(Run))) + geom_point(aes(colour=factor(Run)), size=1) + #geom_smooth(se=FALSE, linetype="dashed", size=0.5) + scale_color_manual(values = global_colors) + facet_wrap(~Threads) + labs(y = paste0("Average Memory Across Runs", "\n", "(maxrss (GB))"), x="") + scale_y_continuous() + theme_bw() + theme(legend.title = element_blank()) + theme(strip.background =element_rect(fill="white", colour = "white"))+ theme(strip.text = element_text(colour = 'black',size=10), strip.text.x = element_markdown(hjust = 0.5)) + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) merged_df %>% #mutate(Sample = fct_reorder2(Sample, Threads, mem_order, .desc = FALSE)) %>% ggplot(aes(x = factor(Sample, levels=order_df), y = average_maxrss, group=Run)) + geom_point(aes(color=factor(Threads),shape=factor(Run),), size=3, alpha=0.9)+ #geom_line()+ scale_shape_manual(values=c("STAR"=3, "WASP"=1, "STAR_WASP"=17))+ #WASP 16 #"STAR"=4, "WASP"=1, "STAR_WASP"=17 #scale_color_manual(values=c('#999999','#E69F00', '#56B4E9'))+ scale_color_manual(values = global_colors) + # facet_wrap(~Threads) + #labs(y = "Mapping Speed (Million of reads/hour", x="") + labs(y = paste0("Average Memory Across Runs", "\n", "(maxrss (GB))"), x="") + #scale_y_continuous(limits = c(30,40)) + theme_bw() + theme(legend.title = element_blank()) + #Base, bw, excel_new(), few, light, lindraw theme(strip.background =element_rect(fill="white", colour = "white"))+ theme(strip.text = element_text(colour = 'black'), strip.text.x = element_markdown(hjust = 0)) + #theme(strip.text.x = element_blank()) + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) #+ merged_df %>% ggplot(aes(x = factor(Sample, levels=order_df), y = (average_ru_maxrss)/1000000, group=Run,color=factor(Run))) + #ru_maxrss is measured in kb, we'll convert it to GB geom_point(aes(colour=factor(Run)), size=1) + scale_color_manual(values = global_colors) + facet_wrap(~Threads) + labs(y = paste0("Average Memory Across Runs", "\n", "(ru_maxrss (GB))"), x="") + scale_y_continuous() + theme_bw() + theme(legend.title = element_blank()) + theme(strip.background =element_rect(fill="white", colour = "white"))+ theme(strip.text = element_text(colour = 'black',size=10), strip.text.x = element_markdown(hjust = 0.5)) + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) ``` ### Average wall clock across all runs ```{r} merged_df %>% ggplot(aes(x = factor(Sample, levels=order_df), y = average_wallclock, group=Run,color=factor(Run))) + geom_point(aes(colour=factor(Run)), size=1) + scale_color_manual(values = global_colors) + facet_wrap(~Threads) + labs(y = paste0("Average Wall Clock Across Runs","\n", "(wallclock)"), x="") + scale_y_continuous() + theme_bw() + theme(legend.title = element_blank()) + theme(strip.background =element_rect(fill="white", colour = "white"))+ theme(strip.text = element_text(colour = 'black',size=10), strip.text.x = element_markdown(hjust = 0.5)) + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) merged_df %>% ggplot(aes(x =factor(Sample, levels=order_df), y = average_ru_wallclock, group=Run,color=factor(Run))) + geom_point(aes(colour=factor(Run)), size=1) + scale_color_manual(values = global_colors) + facet_wrap(~Threads) + facet_wrap(~Threads) + labs(y = paste0("Average Wall Clock Across Runs", "\n", "(ru_wallclock)"), x="") + scale_y_continuous() + theme_bw() + theme(legend.title = element_blank()) + theme(strip.background =element_rect(fill="white", colour = "white"))+ theme(strip.text = element_text(colour = 'black',size=10), strip.text.x = element_markdown(hjust = 0.5)) + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) ``` ### MIllion Reads per Hour ```{r} #Loading pre-processed data from final log results (all runs) data_melted_beta <- read.csv("/home/asiimwe/projects/run_env/alpha_star_wasp_benchmarking/Downstream_Analysis/data_melted_beta.csv") #data_melted_beta <- read.csv(pipe('ssh atlas "cat /home/asiimwe/projects/run_env/alpha_star_wasp_comparison/data_melted_beta.csv"')) data_melted_beta$X <- NULL colnames(data_melted_beta) #"Sample" "Param" "Value" "Threads" "Run" "Mapping_speed_Million_of_reads_per_hour" "Average_input_read_length" #Filtering samples not needed downstream: unique(data_melted_beta$Sample) data_melted_beta$Threads <- as.character(data_melted_beta$Threads ) data_melted_beta$Threads[data_melted_beta$Threads == "8threads" ] <- "8 Threads" data_melted_beta$Threads[data_melted_beta$Threads == "16threads" ] <- "16 Threads" data_melted_beta$Threads[data_melted_beta$Threads == "32threads" ] <- "32 Threads" unique(data_melted_beta$Threads) #Data cleanup str(data_melted_beta) data_melted_beta$Sample <- as.character(data_melted_beta$Sample) data_melted_beta$Param <- as.character(data_melted_beta$Param) data_melted_beta$Value <- gsub("%", "", data_melted_beta$Value) #Reordering threads to appear in the correct order (8, 16, 32) data_melted_beta$Threads <- as.factor(data_melted_beta$Threads) data_melted_beta$Threads <- ordered(data_melted_beta$Threads , levels = c("8 Threads", "16 Threads", "32 Threads")) ## Cleaning up Sample names unique(data_melted_beta$Sample) #Removing NA12878 from all NA12878 derived samples data_melted_beta$Sample <- gsub("NA12878_", "", data_melted_beta$Sample) data_melted_beta$Sample <- gsub("Rep", "R2", data_melted_beta$Sample) data_melted_beta$Sample[data_melted_beta$Sample == "Total"] <- "Total_R1" data_melted_beta$Sample[data_melted_beta$Sample == "PolyA"] <- "PolyA_R1" data_melted_beta$Sample[data_melted_beta$Sample == "Nucleus_PolyA"] <- "Nucleus_PolyA_R1" data_melted_beta$Sample[data_melted_beta$Sample == "Nucleus_nonPolyA"] <- "Nucleus_nonPolyA_R1" ``` ```{r} data_melted_beta_wall_clock <- data_melted_beta %>% filter(Param == "Wall_Clock") # Extracting number of input reads per sample based on STAR runs - should be same inwput reads for each sample/fastq files num_input_reads_per_sample <- unique(data_melted_beta %>% filter(Run == "STAR") %>% dplyr::select("Sample", "Number_of_input_reads")) colnames(num_input_reads_per_sample)[2] <- "Number_of_input_reads_initial" data_melted_beta_wall_clock <- inner_join(data_melted_beta_wall_clock, num_input_reads_per_sample, by = "Sample") ### Number of input reads per sample: num_input_reads_per_sample %>% ggplot(aes(x=factor(Sample, levels=order_df), y=Number_of_input_reads_initial)) + geom_bar(stat = "identity", fill="gray60") + theme_light() + scale_y_continuous(expand=c(0,0), limits=c(0, max(num_input_reads_per_sample$Number_of_input_reads_initial)+1000000))+ scale_x_discrete(expand=c(0,0)) + labs(y="Number of Input Reads", x = "") + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) # Plotting by average input read length input_read_len_per_sample <- unique(data_melted_beta %>% filter(Run == "STAR") %>% dplyr::select("Sample", "Average_input_read_length")) input_read_len_per_sample %>% ggplot(aes(x=factor(Sample, levels=order_df), y=Average_input_read_length)) + geom_bar(stat = "identity", fill="gray60") + theme_light() + scale_y_continuous(expand=c(0,0), limits = c(0, max(input_read_len_per_sample$Average_input_read_length)+10))+ scale_x_discrete(expand=c(0,0)) + labs(y="Average Input Read Length", x = "") + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) min(num_input_reads_per_sample$Number_of_input_reads_initial)#56254714 max(num_input_reads_per_sample$Number_of_input_reads_initial)#128402941 mean(num_input_reads_per_sample$Number_of_input_reads_initial)#86852506 median(num_input_reads_per_sample$Number_of_input_reads_initial)#92285172 data_melted_beta_wall_clock$Unique_ID <- 1:nrow(data_melted_beta_wall_clock) data_melted_beta_wall_clock$Unique_ID <- as.character(data_melted_beta_wall_clock$Unique_ID) data_melted_beta_wall_clock$Wall_Clock <- data_melted_beta_wall_clock$Value #maintaining a copy of the original variable num_input_reads <- unique(data_melted_beta_wall_clock[,c(1,9)]) merged_df <- inner_join(merged_df, num_input_reads, by = "Sample") #Converting wall clock to hours (reported wall clock is in seconds) merged_df <- merged_df %>% mutate(wall_clock_seconds_in_hrs = (average_wallclock/(60*60))) #0.09966657 0.19093528 0.12747778 merged_df <- merged_df %>% mutate(reads_per_hour = (Number_of_input_reads_initial/wall_clock_seconds_in_hrs)) # Plotting wall clock logticks <- function(datavar,type) { minimum <- 1/10^abs(floor(log10(min(datavar, na.rm=TRUE)))) maximum <- 1*10^abs(floor(log10(max(datavar, na.rm=TRUE)))+1) multiple <- floor(log10(maximum/minimum)) yourtickvector <- c() if (type=="breaks") { yourtickvector <- c(minimum) for (x in seq(0,multiple)) { andadd <- seq(minimum*10^x,minimum*10^(x+1),minimum*10^x)[-1] yourtickvector <- c(yourtickvector,andadd) } } else if (type=="labels") { for (x in seq(0,multiple)) { andadd <- c(minimum*10^x,rep("",8)) yourtickvector <- c(yourtickvector,andadd) } yourtickvector <- c(yourtickvector,minimum*10^multiple) } return(yourtickvector) } merged_df %>% ggplot(aes(x = factor(Sample, levels=order_df), y = reads_per_hour/1000000, group=Run, color=Run)) + geom_point(aes(colour=factor(Run)), size=1) + scale_color_manual(values = global_colors) + facet_wrap(~Threads) + labs(y = "Million Reads/Hour", x = "")+ theme_bw() + theme(legend.title = element_blank()) + #scale_y_continuous(trans='log10') + theme(strip.background = element_rect(fill="white", colour = "white"))+ theme(strip.text = element_text(colour = 'black'), strip.text.x = element_markdown(hjust = 0.5, size=12)) + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) + scale_y_log10( breaks = logticks(merged_df$reads_per_hour,"breaks"), labels = logticks(merged_df$reads_per_hour,"labels") ) merged_df %>% ggplot(aes(x = factor(Sample, levels=order_df), y = reads_per_hour/1000000, group=Run, color=Run)) + geom_point(aes(color=factor(Run),shape=factor(Run),fill=factor(Run)), size=2)+#,alpha=.8 scale_shape_manual(values=c(3, 17, 16))+ facet_wrap(~Threads) + scale_color_manual(values=global_colors)+ labs(y = paste0("Average Speed Across Runs", "\n", "(Million Reads/Hour)"), x="")+ theme_bw() + theme(legend.title = element_blank())+ scale_y_continuous(trans='log10') + theme(strip.background =element_rect(fill="white", colour = "white"))+ theme(strip.text = element_text(colour = 'black'), strip.text.x = element_markdown(hjust = 0.5, size=12)) + theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5, size=10)) + scale_y_log10( breaks = logticks(merged_df$reads_per_hour,"breaks"), labels = logticks(merged_df$reads_per_hour,"labels") ) ``` ### Prepairing data for Manuscript summary table ```{r} #STAR merged_df_STAR <- merged_df %>% filter(Run == "STAR") merged_df_STAR_8 <- merged_df_STAR %>% filter(Threads == "8threads") merged_df_STAR_8[,c(1,3,4,5)] <- NULL colnames(merged_df_STAR_8) <- c("Sample", "Average Wall Clock 8 Threads STAR (sec)", "Average Memory 8 Threads STAR", "Number of Input Reads", "Average Wall Clock 8 Threads (hrs) STAR", "Reads Per Hour 8 Threads STAR") merged_df_STAR_8 <- merged_df_STAR_8[,c(1, 4, 2, 5, 3, 6)] merged_df_STAR_16 <- merged_df_STAR %>% filter(Threads == "16threads") merged_df_STAR_16[,c(1,3,4,5)] <- NULL colnames(merged_df_STAR_16) <- c("Sample", "Average Wall Clock 16 Threads STAR (sec)", "Average Memory 16 Threads STAR", "Number of Input Reads", "Average Wall Clock 16 Threads (hrs) STAR", "Reads Per Hour 16 Threads STAR") merged_df_STAR_16 <- merged_df_STAR_16[,c(1, 2, 5, 3, 6)] merged_df_STAR_32 <- merged_df_STAR %>% filter(Threads == "32threads") merged_df_STAR_32[,c(1,3,4,5)] <- NULL colnames(merged_df_STAR_32) <- c("Sample", "Average Wall Clock 32 Threads STAR (sec)", "Average Memory 32 Threads STAR", "Number of Input Reads", "Average Wall Clock 32 Threads (hrs) STAR", "Reads Per Hour 32 Threads STAR") merged_df_STAR_32 <- merged_df_STAR_32[,c(1, 2, 5, 3, 6)] #STAR_WASP merged_df_STAR_WASP <- merged_df %>% filter(Run == "STAR_WASP") merged_df_STAR_WASP_8 <- merged_df_STAR_WASP %>% filter(Threads == "8threads") merged_df_STAR_WASP_8[,c(1,3,4,5)] <- NULL colnames(merged_df_STAR_WASP_8) <- c("Sample", "Average Wall Clock 8 Threads STAR_WASP (sec)", "Average Memory 8 Threads STAR_WASP", "Number of Input Reads", "Average Wall Clock 8 Threads (hrs) STAR_WASP", "Reads Per Hour 8 Threads STAR_WASP") merged_df_STAR_WASP_8 <- merged_df_STAR_WASP_8[,c(1, 2, 5, 3, 6)] merged_df_STAR_WASP_16 <- merged_df_STAR_WASP %>% filter(Threads == "16threads") merged_df_STAR_WASP_16[,c(1,3,4,5)] <- NULL colnames(merged_df_STAR_WASP_16) <- c("Sample", "Average Wall Clock 16 Threads STAR_WASP (sec)", "Average Memory 16 Threads STAR_WASP", "Number of Input Reads", "Average Wall Clock 16 Threads (hrs) STAR_WASP", "Reads Per Hour 16 Threads STAR_WASP") merged_df_STAR_WASP_16 <- merged_df_STAR_WASP_16[,c(1, 2, 5, 3, 6)] merged_df_STAR_WASP_32 <- merged_df_STAR_WASP %>% filter(Threads == "32threads") merged_df_STAR_WASP_32[,c(1,3,4,5)] <- NULL colnames(merged_df_STAR_WASP_32) <- c("Sample", "Average Wall Clock 32 Threads STAR_WASP (sec)", "Average Memory 32 Threads STAR_WASP", "Number of Input Reads", "Average Wall Clock 32 Threads (hrs) STAR_WASP", "Reads Per Hour 32 Threads STAR_WASP") merged_df_STAR_WASP_32 <- merged_df_STAR_WASP_32[,c(1, 2, 5, 3, 6)] #WASP merged_df_WASP <- merged_df %>% filter(Run == "WASP") merged_df_WASP_8 <- merged_df_WASP %>% filter(Threads == "8threads") merged_df_WASP_8[,c(1,3,4,5)] <- NULL colnames(merged_df_WASP_8) <- c("Sample", "Average Wall Clock 8 Threads WASP (sec)", "Average Memory 8 Threads WASP", "Number of Input Reads", "Average Wall Clock 8 Threads (hrs) WASP", "Reads Per Hour 8 Threads WASP") merged_df_WASP_8 <- merged_df_WASP_8[,c(1, 2, 5, 3, 6)] merged_df_WASP_16 <- merged_df_WASP %>% filter(Threads == "16threads") merged_df_WASP_16[,c(1,3,4,5)] <- NULL colnames(merged_df_WASP_16) <- c("Sample", "Average Wall Clock 16 Threads WASP (sec)", "Average Memory 16 Threads WASP", "Number of Input Reads", "Average Wall Clock 16 Threads (hrs) WASP", "Reads Per Hour 16 Threads WASP") merged_df_WASP_16 <- merged_df_WASP_16[,c(1, 2, 5, 3, 6)] merged_df_WASP_32 <- merged_df_WASP %>% filter(Threads == "32threads") merged_df_WASP_32[,c(1,3,4,5)] <- NULL colnames(merged_df_WASP_32) <- c("Sample", "Average Wall Clock 32 Threads WASP (sec)", "Average Memory 32 Threads WASP", "Number of Input Reads", "Average Wall Clock 32 Threads (hrs) WASP", "Reads Per Hour 32 Threads WASP") merged_df_WASP_32 <- merged_df_WASP_32[,c(1, 2, 5, 3, 6)] latex_merged_df<- merged_df_STAR_8 %>% inner_join(merged_df_STAR_8) %>% inner_join(merged_df_STAR_16, by = "Sample") %>% inner_join(merged_df_STAR_32, by = "Sample") %>% inner_join(merged_df_STAR_WASP_8, by = "Sample") %>% inner_join(merged_df_STAR_WASP_16, by = "Sample") %>% inner_join(merged_df_STAR_WASP_32, by = "Sample") %>% inner_join(merged_df_WASP_8, by = "Sample") %>% inner_join(merged_df_WASP_16, by = "Sample") %>% inner_join(merged_df_WASP_32, by = "Sample") write.csv(latex_merged_df, file= "/home/asiimwe/projects/run_env/alpha_star_wasp_benchmarking/Downstream_Analysis/elzar_run_results_averaged.csv") ```