From bdbc97c40a1d19fd60e6bed4d129521f3c8f239b Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Wed, 3 Aug 2022 18:56:14 +0100 Subject: [PATCH] fix many plot functions to stop them using the "g=ggplot()" pattern, which annoyingly throws away lots of useful data that RShiny needs for clickable plots. Also split the "flame bar" for ligand distance out into separate functions in generate_distance_colour_map.R. This can now be easily incorporated into any "wide" graph showing all positions. --- scripts/functions/bp_subcolours.R | 88 ++- scripts/functions/consurfP.R | 691 ++++++++---------- .../functions/generate_distance_colour_map.R | 121 +++ scripts/functions/lineage_dist.R | 80 +- scripts/functions/logoP_or.R | 193 ++--- scripts/functions/logoP_snp.R | 358 +++++---- scripts/functions/position_count_bp.R | 11 +- scripts/functions/wideP_consurf.R | 506 +++++++++++++ 8 files changed, 1323 insertions(+), 725 deletions(-) create mode 100644 scripts/functions/generate_distance_colour_map.R create mode 100644 scripts/functions/wideP_consurf.R diff --git a/scripts/functions/bp_subcolours.R b/scripts/functions/bp_subcolours.R index 18089c9..5adc071 100755 --- a/scripts/functions/bp_subcolours.R +++ b/scripts/functions/bp_subcolours.R @@ -2,6 +2,7 @@ # 1b: Define function: coloured barplot by subgroup # LINK: https://stackoverflow.com/questions/49818271/stacked-barplot-with-colour-gradients-for-each-bar ######################################################### +source("~/git/LSHTM_analysis/scripts/functions/generate_distance_colour_map.R") ColourPalleteMulti = function(df, group, subgroup){ @@ -35,6 +36,7 @@ ColourPalleteMulti = function(df, group, subgroup){ bp_stability_hmap <- function(plotdf = merged_df3 , xvar_colname = "position" + , yvar_colname = 'duet_scaled' #FIXME: temp, remove #, bar_col_colname = "group" , stability_colname = "duet_scaled" , stability_outcome_colname = "duet_outcome" @@ -46,8 +48,27 @@ bp_stability_hmap <- function(plotdf = merged_df3 , my_pts = 20 # plot-title size , my_xlab = "Position" , my_ylab = "No. of nsSNPs" + + # Custom 2: x-axis: geom tiles ~ lig distance + #, A_xvar_lig = T + , lig_dist_colname = LigDist_colname # from globals + , tpos0 = 0 # 0 is a magic number that does my sensible default + , tW0 = 1 + , tH0 = 0.2 + + + ) { + ################################################ + # Custom 2: x-axis geom tiles ~ lig distance + ################################################ + + #========================= + # Build data with colours + # ~ ligand distance + #========================= + plotdf = generate_distance_colour_map(plotdf, yvar_colname = stability_colname, debug=TRUE) # order the df by position and ensure it is a factor plotdf = plotdf[order(plotdf[[xvar_colname]]), ] @@ -74,38 +95,45 @@ bp_stability_hmap <- function(plotdf = merged_df3 cat("\nNo. of sub colours generated:", length(subcols_bp_hmap)) + #------------------------------- # Generate the subcols barplot #------------------------------- - - #g = ggplot(plotdf, aes(x = factor(position, ordered = T))) - g = ggplot(plotdf, aes_string(x = xvar_colname - # , ordered = T) - )) - - - OutWidePlot = g + geom_bar(aes(fill = group) - , colour = "grey") + - - scale_fill_manual( values = subcols_bp_hmap - , guide = "none") + - - theme( axis.text.x = element_text(size = my_xaxls - , angle = 90 - , hjust = 1 - , vjust = 0.4) - , axis.text.y = element_text(size = my_yaxls - , angle = 0 + cowplot::plot_grid( + ggplot(plotdf, aes_string(x = xvar_colname + # , ordered = T) + )) + + geom_bar(aes(fill = group) + , colour = "grey") + + + scale_fill_manual( values = subcols_bp_hmap + , guide = "none") + + + theme( axis.text.x = element_text(size = my_xaxls + , angle = 90 , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = my_xaxts) - , axis.title.y = element_text(size = my_yaxts ) - , plot.title = element_text(size = my_pts - , hjust = 0.5)) + - - labs(title = p_title - , x = my_xlab - , y = my_ylab) - - return(OutWidePlot) + , vjust = 0.4) + , axis.text.y = element_text(size = my_yaxls + , angle = 0 + , hjust = 1 + , vjust = 0) + , axis.title.x = element_text(size = my_xaxts) + , axis.title.y = element_text(size = my_yaxts ) + , plot.title = element_text(size = my_pts + , hjust = 0.5)) + + geom_tile(aes(, tpos0 # heat-mapped distance tiles along the bot + , width = tW0 + , height = tH0) + , fill = plotdf$ligD_colours + , colour = plotdf$ligD_colours + , linetype = "blank") + + + labs(title = p_title + , x = my_xlab + , y = my_ylab), + generate_distance_legend(plotdf, yvar_colname = stability_colname), + ncol = 2, + #align = "hv", + rel_widths = c(9/10, 0.4/10) + ) } diff --git a/scripts/functions/consurfP.R b/scripts/functions/consurfP.R index d806e0b..7be5d9c 100644 --- a/scripts/functions/consurfP.R +++ b/scripts/functions/consurfP.R @@ -19,7 +19,7 @@ wideP_consurf <- function(plotdf , plot_error_bars = T , upper_EB_colname = "consurf_ci_upper" , lower_EB_colname = "consurf_ci_lower" - + , plot_type = "point" # default is point , point_colours , p_size = 2 @@ -30,7 +30,7 @@ wideP_consurf <- function(plotdf , "9": "Conserved") , panel_col = "black" , panel_col_fill = "black" - + # axes title and label sizes , x_axls = 12 # x-axis label size , y_axls = 15 # y-axis label size @@ -41,20 +41,20 @@ wideP_consurf <- function(plotdf , xlab = "" , ylab = "" , pts = 20 - + # plot margins , t_margin = 0.5 , r_margin = 0.5 , b_margin = 1 , l_margin = 1 , unit_margin = "cm" - + # Custom 1: x-axis: text colour , xtext_colour_aa = F , xtext_colour_aa1 = active_aa_pos , xtext_colour_aa2 = aa_pos_drug , xtext_colours = c("purple", "brown", "black") - + # Custom 2: x-axis: geom tiles ~ lig distance , A_xvar_lig = T , leg_title2 = "Ligand Distance" @@ -63,7 +63,7 @@ wideP_consurf <- function(plotdf , tpos0 = 0 # 0 is a magic number that does my sensible default , tW0 = 1 , tH0 = 0.3 - + # Custom 3: x-axis: geom tiles ~ active sites and ligand , A_xvar_aa = F , aa_pos_drug = NULL @@ -72,34 +72,34 @@ wideP_consurf <- function(plotdf , tH = 0.2 , active_aa_pos = NULL , active_aa_colour = "brown" - + , aa_pos_lig1 = NULL , aa_colour_lig1 = "blue" , tpos1 = 0 - + , aa_pos_lig2 = NULL , aa_colour_lig2 = "cyan" , tpos2 = 0 - + , aa_pos_lig3 = NULL , aa_colour_lig3 = "cornflowerblue" , tpos3 = 0 - + , default_gt_clr = "white" , debug=FALSE -){ - + ){ + if(missing(point_colours)){ temp_cols = colorRampPalette(c("seagreen", "sienna3"))(30) point_colours = temp_cols }else{ point_colours = point_colours } - + ############################### # custom 1: x-axis text colour ############################## - + if (xtext_colour_aa){ positionF <- levels(as.factor(plotdf[[xvar_colname]])) length(positionF) @@ -110,42 +110,42 @@ wideP_consurf <- function(plotdf }else{ aa_pos_colours = default_xtc } - + ################################################ # Custom 2: x-axis geom tiles ~ lig distance ################################################ - + #========================= # Build data with colours # ~ ligand distance #========================= if (A_xvar_lig){ cat("\nAnnotating x-axis ~", lig_dist_colname, "requested...") - + #------------------------------------- # round column values: to colour by #-------------------------------------- #plotdf = plotdf[order(plotdf[[lig_dist_colname]]),] plotdf['lig_distR'] = round(plotdf[[lig_dist_colname]], digits = 0) head(plotdf['lig_distR']) - + #------------------------------------- # ligand distance range, min, max, etc #-------------------------------------- lig_min = min(round(plotdf[[lig_dist_colname]]), na.rm = T); lig_min lig_max = max(round(plotdf[[lig_dist_colname]]), na.rm = T); lig_max lig_mean = round(mean(round(plotdf[[lig_dist_colname]]), na.rm = T)); lig_mean - + #------------------------------------- # Create mapping colour key #-------------------------------------- # sorting removes NA, so that n_colours == length(ligD_valsR) n_colours = length(sort(unique(round(plotdf[[lig_dist_colname]], digits = 0)))); n_colours - + lig_cols = colorRampPalette(lig_dist_colours)(n_colours); lig_cols ligD_valsR = sort(unique(round(plotdf[[lig_dist_colname]], digits = 0))); ligD_valsR length(ligD_valsR) - + if (n_colours == length(ligD_valsR)) { cat("\nStarting: mapping b/w" , lig_dist_colname @@ -156,17 +156,17 @@ wideP_consurf <- function(plotdf , "No. of colours: ", n_colours , "\nValues to map:", length(ligD_valsR)) } - + ligDcolKey <- data.frame(ligD_colours = lig_cols , lig_distR = ligD_valsR); ligDcolKey names(ligDcolKey) cat("\nSuccessful: Mapping b/w", lig_dist_colname, "and colours") - + #------------------------------------- # merge colour key with plotdf #-------------------------------------- plotdf = merge(plotdf, ligDcolKey, by = 'lig_distR') - + plotdf_check = as.data.frame(cbind(position = plotdf[[xvar_colname]] , ligD = plotdf[[lig_dist_colname]] , ligDR = plotdf$lig_distR @@ -174,24 +174,24 @@ wideP_consurf <- function(plotdf } else{ plotdf = plotdf } - + ############################################### # Custom 3: x-axis geom tiles ~ active sites ################################################ - + #========================== # Build Data with colours # ~ on active sites #========================== - + if(A_xvar_aa) { cat("\nAnnotation for xvar requested. Building colours for annotation...") - + aa_colour_colname = "bg_all" aa_colour_colname1 = "col_bg1" aa_colour_colname2 = "col_bg2" aa_colour_colname3 = "col_bg3" - + #-------------------------------------------------- # column colour 0: Active site + drug binding sites #-------------------------------------------------- @@ -201,7 +201,7 @@ wideP_consurf <- function(plotdf , active_aa_colour, default_gt_clr )) plotdf[[aa_colour_colname]] cat("\nColumn created 'bg_all':", length(plotdf[[aa_colour_colname]])) - + #------------------------------------------------ # column colour 1: Ligand 1 + drug binding sites #------------------------------------------------ @@ -219,19 +219,19 @@ wideP_consurf <- function(plotdf plotdf[[aa_colour_colname2]] = plotdf[[aa_colour_colname1]] plotdf[[aa_colour_colname2]] = ifelse(plotdf[[xvar_colname]]%in%aa_pos_lig2 , aa_colour_lig2, plotdf[[aa_colour_colname1]]) - + #------------------------------------------------ # column colour 3: Ligand 3 #------------------------------------------------ plotdf[[aa_colour_colname3]] = plotdf[[aa_colour_colname2]] plotdf[[aa_colour_colname3]] = ifelse(plotdf[[xvar_colname]]%in%aa_pos_lig3 , aa_colour_lig3, plotdf[[aa_colour_colname2]]) - + } ################### # start plot ################### - + #------------------- # x and y axis # range, scale, etc @@ -239,30 +239,12 @@ wideP_consurf <- function(plotdf my_xlim = length(unique(plotdf[[xvar_colname]])); my_xlim ymin = min(plotdf[[yvar_colname]]); ymin ymax = max(plotdf[[yvar_colname]]); ymax - # DEBUG - # x_factor = as.factor(plotdf[[xvar_colname]]) - # y_colour_factor = as.factor(plotdf[[yvar_colourN_colname]]) - - # cat('XXXX THIS:', str(x_factor)) - # - # g = ggplot(plotdf, aes_string(x = x_factor) - # , y = yvar_colname - # , colour = y_colour_factor) - + g = ggplot(plotdf, aes_string(x = sprintf("factor(%s)", xvar_colname) , y = yvar_colname , colour = sprintf("factor(%s)", yvar_colourN_colname) - )) - # g = ggplot(plotdf, aes_string(x = "position") - # , y = "consurf_score" - # , colour = "consurf_colour_rev") - # g = ggplot(plotdf, aes_string(x = "position") - # , y = "consurf_score" - # , colour = "consurf_colour_rev") - - # DEBUG - #g0=g - + )) + "if SPECIAL do SPECIAL THING, otherwise do NORMAL THING" if (plot_type == "bar"){ g0 = g + @@ -273,23 +255,21 @@ wideP_consurf <- function(plotdf coord_cartesian(xlim = c(1, my_xlim) , ylim = c(ymin, ymax) , clip = "off") + - geom_point(size = p_size) + - scale_colour_manual(values = point_colours) - # , labels = leg_labels - # , name = leg_title1) +geom_point(size = p_size) + +scale_colour_manual(values = point_colours) } - + if (plot_error_bars){ g0 = g0 + geom_errorbar(aes(ymin = eval(parse(text = lower_EB_colname)) , ymax = eval(parse(text = upper_EB_colname)) - )) + )) }else{ - + g0 = g0 - + } - + #--------------------- # add axis formatting #--------------------- @@ -299,327 +279,304 @@ wideP_consurf <- function(plotdf , vjust = 0.4 , face = "bold" , colour = aa_pos_colours) - , axis.text.y = element_text(size = y_axts - , angle = 0 - , hjust = 1 - , vjust = 0) - , axis.title.x = element_text(size = x_axls) - , axis.title.y = element_text(size = y_axls ) - , panel.background = element_rect(fill = panel_col_fill, color = panel_col) - , panel.grid.major = element_line(color = "black") - , panel.grid.minor = element_line(color = "black") - , plot.title = element_text(size = pts - , hjust = 0.5) - , plot.margin = margin(t = t_margin - , r = r_margin - , b = b_margin - , l = l_margin - , unit = unit_margin))+ - guides(colour = guide_legend(title = "ConsurfXXXX")) + - - labs(title = ptitle - , x = xlab - , y = ylab) - - #------------------ - #Extract legend1 - #------------------ - # yayy - g1_leg = ggplot(plotdf, aes_string(x = sprintf("factor(%s)" - , xvar_colname) )) - g1_leg = g1_leg + geom_bar(); g1_leg - g1_leg = g1_leg + geom_bar(aes_string(fill = sprintf("factor(%s)" - , yvar_colourN_colname))) - - g1_leg = g1_leg + scale_fill_manual(values = consurf_palette2 , name = leg_title1) - g1_leg - - legend1 = get_legend(g1_leg) - - - ##################################################### - #============================================ - # x-axis: geom_tiles ~ ligand distance - #============================================ - #------- - # plot - #------- - if(A_xvar_lig){ # 0 is a magic number that does my sensible default - if (tpos0 == 0){ - tpos0 = ymin-0.5 - } - if (tpos1 == 0){ - tpos1 = ymin-0.65 - } - if (tpos2 == 0){ - tpos2 = ymin-0.75 - } - if (tpos3 == 0){ - tpos3 = ymin-0.85 - } - - - cat("\nColouring x-axis aa based on", lig_dist_colname - , "\nNo. of colours:", n_colours) - - g2 = g1 + geom_tile(aes(, tpos0 - , width = tW0 - , height = tH0) - , fill = plotdf$ligD_colours - , colour = plotdf$ligD_colours - , linetype = "blank") - - #cat("Nrows of plot df", length(plotdf$ligD_colours)) - out = g2 - # - # #------------------ - # # Extract legend2 - # #------------------ - # labels = seq(lig_min, lig_max, len = 5); labels - # labelsD = round(labels, digits = 0); labelsD - # - # g2_leg = g1 + - # geom_tile(aes(fill = .data[[lig_dist_colname]]) - # , colour = "white") + - # scale_fill_gradient2(midpoint = lig_mean - # , low = "green" - # , mid = "yellow" - # , high = "red" - # , breaks = labels - # #, n.breaks = 11 - # #, minor_breaks = c(2, 4, 6, 8, 10) - # , limits = c(lig_min, lig_max) - # , labels = labelsD - # , name = leg_title2) - # - # legend2 = get_legend(g2_leg) - # - # }else{ - # out = g1 - # } - ###################################################### - #------------------ - # Extract legend2 - #------------------ - labels = seq(lig_min, lig_max, len = 5); labels - labelsD = round(labels, digits = 0); labelsD - g2_leg = ggplot(plotdf, aes_string(x = sprintf("factor(%s)", xvar_colname) - , y = yvar_colname) - ) + - geom_tile(aes(fill = .data[[lig_dist_colname]]) - , colour = "white") + - scale_fill_gradient2(midpoint = lig_mean - , low = "green" - , mid = "yellow" - , high = "red" - , breaks = labels - #, n.breaks = 11 - #, minor_breaks = c(2, 4, 6, 8, 10) - , limits = c(lig_min, lig_max) - , labels = labelsD - , name = leg_title2) - - legend2 = get_legend(g2_leg) - - }else{ - out = g1 - } - - ##################################################### - # #============================================ - # # x-axis: geom_tiles ~ ligand distance - # #============================================ - # if(A_xvar_lig){ - # cat("\nColouring x-axis aa based on", lig_dist_colname - # , "\nNo. of colours:", n_colours) + , axis.text.y = element_text(size = y_axts + , angle = 0 + , hjust = 1 + , vjust = 0) + , axis.title.x = element_text(size = x_axls) + , axis.title.y = element_text(size = y_axls ) + , panel.background = element_rect(fill = panel_col_fill, color = panel_col) + , panel.grid.major = element_line(color = "black") + , panel.grid.minor = element_line(color = "black") + , plot.title = element_text(size = pts + , hjust = 0.5) + , plot.margin = margin(t = t_margin + , r = r_margin + , b = b_margin + , l = l_margin + , unit = unit_margin))+ +guides(colour = guide_legend(title = "ConsurfXXXX")) + + +labs(title = ptitle + , x = xlab + , y = ylab) + +#------------------ +#Extract legend1 +#------------------ +# yayy +g1_leg = ggplot(plotdf, aes_string(x = sprintf("factor(%s)" + , xvar_colname) )) +g1_leg = g1_leg + geom_bar(); g1_leg +g1_leg = g1_leg + geom_bar(aes_string(fill = sprintf("factor(%s)" + , yvar_colourN_colname))) + +g1_leg = g1_leg + scale_fill_manual(values = consurf_palette2 , name = leg_title1) +g1_leg + +legend1 = get_legend(g1_leg) + + +##################################################### +#============================================ +# x-axis: geom_tiles ~ ligand distance +#============================================ +#------- +# plot +#------- +if(A_xvar_lig){ # 0 is a magic number that does my sensible default + if (tpos0 == 0){ + tpos0 = ymin-0.5 + } + if (tpos1 == 0){ + tpos1 = ymin-0.65 + } + if (tpos2 == 0){ + tpos2 = ymin-0.75 + } + if (tpos3 == 0){ + tpos3 = ymin-0.85 + } + + + cat("\nColouring x-axis aa based on", lig_dist_colname + , "\nNo. of colours:", n_colours) + + g2 = g1 + geom_tile(aes(, tpos0 + , width = tW0 + , height = tH0) + , fill = plotdf$ligD_colours + , colour = plotdf$ligD_colours + , linetype = "blank") + + #cat("Nrows of plot df", length(plotdf$ligD_colours)) + out = g2 # - # g2 = g1 + geom_tile(aes(, tpos0 - # , width = tW0 - # , height = tH0) - # , fill = plotdf$ligD_colours - # , colour = plotdf$ligD_colours - # , linetype = "blank") - # - # #cat("Nrows of plot df", length(plotdf$ligD_colours)) - # out = g2 - # }else{ + # #------------------ + # # Extract legend2 + # #------------------ + # labels = seq(lig_min, lig_max, len = 5); labels + # labelsD = round(labels, digits = 0); labelsD + # + # g2_leg = g1 + + # geom_tile(aes(fill = .data[[lig_dist_colname]]) + # , colour = "white") + + # scale_fill_gradient2(midpoint = lig_mean + # , low = "green" + # , mid = "yellow" + # , high = "red" + # , breaks = labels + # #, n.breaks = 11 + # #, minor_breaks = c(2, 4, 6, 8, 10) + # , limits = c(lig_min, lig_max) + # , labels = labelsD + # , name = leg_title2) + # + # legend2 = get_legend(g2_leg) + # + # }else{ # out = g1 # } - # - ############################################################################################ - #============================================== - # x-axis: geom_tiles ~ active sites and others - #============================================== - if(A_xvar_aa){ - #tpos = 0 - #tW = 1 - #tH = 0.2 - - #--------------------- - # Add2plot: 3 ligands - #--------------------- - if (all(!is.null(active_aa_pos) && - !is.null(aa_pos_drug) && - !is.null(aa_pos_lig1) && !is.null(aa_pos_lig2) && !is.null(aa_pos_lig3))) { - if (debug){ - cat("\n\nAnnotating xvar with active, drug binding, and Lig 1&2&3 sites") - cat("\nCreating column colours, column name:", aa_colour_colname3) - - cat("\nDoing Plot with 3 ligands") - } - out = out + geom_tile(aes(,tpos3 - , width = tW - , height = tH ) - , fill = plotdf[[aa_colour_colname3]] - , colour = plotdf[[aa_colour_colname3]] - , linetype = "solid") + - geom_tile(aes(, tpos2 - , width = tW - , height = tH ) - , fill = plotdf[[aa_colour_colname2]] - , colour = plotdf[[aa_colour_colname2]] - , linetype = "solid")+ - - geom_tile(aes(, tpos1 - , width = tW - , height = tH) - , fill = plotdf[[aa_colour_colname1]] - , colour = plotdf[[aa_colour_colname1]] - , linetype = "solid") - if (debug){ - cat("\nDone Plot with 3 ligands") - } + ###################################################### + #------------------ + # Extract legend2 + #------------------ + labels = seq(lig_min, lig_max, len = 5); labels + labelsD = round(labels, digits = 0); labelsD + g2_leg = ggplot(plotdf, aes_string(x = sprintf("factor(%s)", xvar_colname) + , y = yvar_colname) + ) + +geom_tile(aes(fill = .data[[lig_dist_colname]]) + , colour = "white") + +scale_fill_gradient2(midpoint = lig_mean + , low = "green" + , mid = "yellow" + , high = "red" + , breaks = labels + #, n.breaks = 11 + #, minor_breaks = c(2, 4, 6, 8, 10) + , limits = c(lig_min, lig_max) + , labels = labelsD + , name = leg_title2) + +legend2 = get_legend(g2_leg) + +}else{ + out = g1 +} +#============================================== +# x-axis: geom_tiles ~ active sites and others +#============================================== +if(A_xvar_aa){ + #tpos = 0 + #tW = 1 + #tH = 0.2 + + #--------------------- + # Add2plot: 3 ligands + #--------------------- + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + !is.null(aa_pos_lig1) && !is.null(aa_pos_lig2) && !is.null(aa_pos_lig3))) { + if (debug){ + cat("\n\nAnnotating xvar with active, drug binding, and Lig 1&2&3 sites") + cat("\nCreating column colours, column name:", aa_colour_colname3) + + cat("\nDoing Plot with 3 ligands") } - #--------------------- - # Add2plot: 2 ligands - #--------------------- - if (all(!is.null(active_aa_pos) && - !is.null(aa_pos_drug) && - !is.null(aa_pos_lig1) && !is.null(aa_pos_lig2) && is.null(aa_pos_lig3))) { - if (debug){ + out = out + geom_tile(aes(,tpos3 + , width = tW + , height = tH ) + , fill = plotdf[[aa_colour_colname3]] + , colour = plotdf[[aa_colour_colname3]] + , linetype = "solid") + +geom_tile(aes(, tpos2 + , width = tW + , height = tH ) +, fill = plotdf[[aa_colour_colname2]] +, colour = plotdf[[aa_colour_colname2]] +, linetype = "solid")+ + +geom_tile(aes(, tpos1 + , width = tW + , height = tH) +, fill = plotdf[[aa_colour_colname1]] +, colour = plotdf[[aa_colour_colname1]] +, linetype = "solid") +if (debug){ + cat("\nDone Plot with 3 ligands") +} + } + #--------------------- + # Add2plot: 2 ligands + #--------------------- + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + !is.null(aa_pos_lig1) && !is.null(aa_pos_lig2) && is.null(aa_pos_lig3))) { + if (debug){ cat("\n\nAnnotating xvar with active, drug binding, and Lig 1&2 sites") cat("\nCreating column colours, column name:", aa_colour_colname2) - + cat("\nDoing Plot with 2 ligands") - } - out = out + - geom_tile(aes(, tpos2 - , width = tW - , height = tH) - , fill = plotdf[[aa_colour_colname2]] - , colour = plotdf[[aa_colour_colname2]] - , linetype = "solid")+ - geom_tile(aes(, tpos1 - , width = tW - , height = tH) - , fill = plotdf[[aa_colour_colname1]] - , colour = plotdf[[aa_colour_colname1]] - , linetype = "solid") - if (debug){ - cat("\nDone Plot with 2 ligands") - } } - - #--------------------- - # Add2plot: 1 ligand - #--------------------- - if (all(!is.null(active_aa_pos) && - !is.null(aa_pos_drug) && - !is.null(aa_pos_lig1) && is.null(aa_pos_lig2) && is.null(aa_pos_lig3))) { - if (debug){ + out = out + + geom_tile(aes(, tpos2 + , width = tW + , height = tH) + , fill = plotdf[[aa_colour_colname2]] + , colour = plotdf[[aa_colour_colname2]] + , linetype = "solid")+ +geom_tile(aes(, tpos1 + , width = tW + , height = tH) +, fill = plotdf[[aa_colour_colname1]] +, colour = plotdf[[aa_colour_colname1]] +, linetype = "solid") +if (debug){ + cat("\nDone Plot with 2 ligands") +} + } + + #--------------------- + # Add2plot: 1 ligand + #--------------------- + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + !is.null(aa_pos_lig1) && is.null(aa_pos_lig2) && is.null(aa_pos_lig3))) { + if (debug){ cat("\n\nAnnotating xvar with active, drug binding, and Lig 1 sites") cat("\nCreating column colours, column name:", aa_colour_colname1) - + cat("\nDoing Plot with 1 ligands") - } - out = out + - geom_tile(aes(, tpos1 - , width = tW - , height = tH) - , fill = plotdf[[aa_colour_colname1]] - , colour = plotdf[[aa_colour_colname1]] - , linetype = "solid") - - cat("\nDone Plot with 1 ligand") - } - - #----------------------------------- - # Add2plot:NO ligands - # No Ligs: Just drug and active site - # DEFAULT: A_xvar_aa == TRUE - #---------------------------------- - if (all(!is.null(active_aa_pos) && - !is.null(aa_pos_drug) && - is.null(aa_pos_lig1) && - is.null(aa_pos_lig2) && - is.null(aa_pos_lig3))) { - if (debug){ + out = out + + geom_tile(aes(, tpos1 + , width = tW + , height = tH) + , fill = plotdf[[aa_colour_colname1]] + , colour = plotdf[[aa_colour_colname1]] + , linetype = "solid") + + cat("\nDone Plot with 1 ligand") + + } + + #----------------------------------- + # Add2plot:NO ligands + # No Ligs: Just drug and active site + # DEFAULT: A_xvar_aa == TRUE + #---------------------------------- + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + is.null(aa_pos_lig1) && + is.null(aa_pos_lig2) && + is.null(aa_pos_lig3))) { + if (debug){ cat("\n\nAnnotating xvar with active and drug binding sites") cat("\nCreating column colours, column name:", aa_colour_colname) cat("\nDoing Plot with 0 ligands: active and drug site only") - } - out = out + geom_tile(aes(, tpos3 - , width = tW - , height = tH) - , fill = plotdf[[aa_colour_colname]] - , colour = plotdf[[aa_colour_colname]] - , linetype = "solid") - if (debug){ - cat("\nDone Plot with: Active and Drug sites") - } } - }else{ - cat("\nNo annotation for additional ligands on xvar requested") - } - #============================================== - if (A_xvar_lig){ - legs = cowplot::plot_grid(legend1 - , legend2 - , ncol = 1 - , align = "hv" - , rel_heights = c(2/4,3/4)) - - out2 = cowplot::plot_grid( out + theme(legend.position = "none") - , legs - , ncol = 2 - , align = "hv" - , rel_widths = c(9/10, 0.4/10) - ) - }else{ - out2 = cowplot::plot_grid( out + theme(legend.position = "none") - , legend1 - , ncol = 2 - , align = "hv" - , rel_widths = c(9/10, 0.5/10) - ) + out = out + geom_tile(aes(, tpos3 + , width = tW + , height = tH) + , fill = plotdf[[aa_colour_colname]] + , colour = plotdf[[aa_colour_colname]] + , linetype = "solid") + if (debug){ + cat("\nDone Plot with: Active and Drug sites") + } } - #============================================== - - - #============================================== - # if (A_xvar_lig){ - # legs = grid.arrange(legend1 - # , legend2 - # , ncol = 1 - # , heights = c(3/4,1)) - # - # out2 = grid.arrange( out + theme(legend.position = "none") - # , legs - # , ncol = 2 - # , widths = c(9/10, 0.5/10) - # ) - # }else{ - # out2 = grid.arrange( out + theme(legend.position = "none") - # , legend1 - # , ncol = 2 - # , widths = c(9/10, 0.5/10) - # ) - # } - #============================================== - return(out2) - #return(out2) - +}else{ + cat("\nNo annotation for additional ligands on xvar requested") +} +#============================================== +if (A_xvar_lig){ + legs = cowplot::plot_grid(legend1 + , legend2 + , ncol = 1 + , align = "hv" + , rel_heights = c(2/4,3/4)) + + out2 = cowplot::plot_grid( out + theme(legend.position = "none") + , legs + , ncol = 2 + , align = "hv" + , rel_widths = c(9/10, 0.4/10) + ) +}else{ + out2 = cowplot::plot_grid( out + theme(legend.position = "none") + , legend1 + , ncol = 2 + , align = "hv" + , rel_widths = c(9/10, 0.5/10) + ) +} +#============================================== + + +#============================================== +# if (A_xvar_lig){ +# legs = grid.arrange(legend1 +# , legend2 +# , ncol = 1 +# , heights = c(3/4,1)) +# +# out2 = grid.arrange( out + theme(legend.position = "none") +# , legs +# , ncol = 2 +# , widths = c(9/10, 0.5/10) +# ) +# }else{ +# out2 = grid.arrange( out + theme(legend.position = "none") +# , legend1 +# , ncol = 2 +# , widths = c(9/10, 0.5/10) +# ) +# } +#============================================== +return(out2) +#return(out2) + } ############################################################# diff --git a/scripts/functions/generate_distance_colour_map.R b/scripts/functions/generate_distance_colour_map.R new file mode 100644 index 0000000..bf5974a --- /dev/null +++ b/scripts/functions/generate_distance_colour_map.R @@ -0,0 +1,121 @@ +# takes a dataframe and returns the same dataframe with two extra columns for colours and position +generate_distance_colour_map = function(plotdf, + xvar_colname = "position", + yvar_colname = 'duet_scaled', + lig_dist_colname = "ligand_distance", + lig_dist_colours = c("green", "yellow", "orange", "red"), + #tpos0 = 0, + #tpos1 = 0, + #tpos2 = 0, + #tpos3 = 0, + debug = FALSE +) +{ + #------------------- + # x and y axis + # range, scale, etc + #------------------- + my_xlim = length(unique(plotdf[[yvar_colname]])); my_xlim + ymin = min(plotdf[[yvar_colname]]); ymin + ymax = max(plotdf[[yvar_colname]]); ymax + + #if (tpos0 == 0){ + # tpos0 = ymin-0.5 + #} + #if (tpos1 == 0){ + # tpos1 = ymin-0.65 + #} + #if (tpos2 == 0){ + # tpos2 = ymin-0.75 + #} + #if (tpos3 == 0){ + # tpos3 = ymin-0.85 + #} + if (debug) { + cat("\nAnnotating x-axis ~", lig_dist_colname, "requested...") + } + #------------------------------------- + # round column values: to colour by + #-------------------------------------- + #plotdf = plotdf[order(plotdf[[lig_dist_colname]]),] + plotdf['lig_distR'] = round(plotdf[[lig_dist_colname]], digits = 0) + #head(plotdf['lig_distR']) + + #------------------------------------- + # ligand distance range, min, max, etc + #-------------------------------------- + lig_min = min(round(plotdf[[lig_dist_colname]]), na.rm = T); lig_min + lig_max = max(round(plotdf[[lig_dist_colname]]), na.rm = T); lig_max + lig_mean = round(mean(round(plotdf[[lig_dist_colname]]), na.rm = T)); lig_mean + + #------------------------------------- + # Create mapping colour key + #-------------------------------------- + # sorting removes NA, so that n_colours == length(ligD_valsR) + n_colours = length(sort(unique(round(plotdf[[lig_dist_colname]], digits = 0)))); n_colours + + lig_cols = colorRampPalette(lig_dist_colours)(n_colours); lig_cols + ligD_valsR = sort(unique(round(plotdf[[lig_dist_colname]], digits = 0))); ligD_valsR + length(ligD_valsR) + + if (debug) { + if (n_colours == length(ligD_valsR)) { + cat("\nStarting: mapping b/w" + , lig_dist_colname + , "and colours") + }else{ + cat("\nCannot start mapping b/w", lig_dist_colname, "and colours..." + , "\nLength mismatch:" + , "No. of colours: ", n_colours + , "\nValues to map:", length(ligD_valsR)) + } + } + + ligDcolKey <- data.frame(ligD_colours = lig_cols + , lig_distR = ligD_valsR); ligDcolKey + names(ligDcolKey) + if (debug) { + cat("\nSuccessful: Mapping b/w", lig_dist_colname, "and colours") + } + #------------------------------------- + # merge colour key with plotdf + #-------------------------------------- + plotdf = merge(plotdf, ligDcolKey, by = 'lig_distR') + + plotdf_check = as.data.frame(cbind(position = plotdf[[xvar_colname]] + , ligD = plotdf[[lig_dist_colname]] + , ligDR = plotdf$lig_distR + , ligD_cols = plotdf$ligD_colours)) + return(plotdf) +} + +generate_distance_legend = function(plotdf, + yvar_colname, + xvar_colname = 'position', + lig_dist_colname = "ligand_distance", + legend_title = "Ligand\nDistance" + ) +{ + # build legend for ligand distance "heat bar" + lig_min = min(round(plotdf[[lig_dist_colname]]), na.rm = T); lig_min + lig_max = max(round(plotdf[[lig_dist_colname]]), na.rm = T); lig_max + lig_mean = round(mean(round(plotdf[[lig_dist_colname]]), na.rm = T)); lig_mean + + labels = seq(lig_min, lig_max, len = 5); labels + labelsD = round(labels, digits = 0); labelsD + + get_legend(ggplot(plotdf, aes_string(x = sprintf("factor(%s)", xvar_colname) + , y = yvar_colname)) + + + geom_tile(aes(fill = .data[[lig_dist_colname]]) + , colour = "white") + + scale_fill_gradient2(midpoint = lig_mean + , low = "green" + , mid = "yellow" + , high = "red" + , breaks = labels + , limits = c(lig_min, lig_max) + , labels = labelsD + , name = legend_title) + ) +} diff --git a/scripts/functions/lineage_dist.R b/scripts/functions/lineage_dist.R index 0ea221f..5f5db56 100644 --- a/scripts/functions/lineage_dist.R +++ b/scripts/functions/lineage_dist.R @@ -16,7 +16,7 @@ lineage_distP <- function(plotdf , all_lineages = F , use_lineages = c("L1", "L2", "L3", "L4") , with_facet = F - , facet_wrap_var = "" + , facet_wrap_var = "" # FIXME: document what this is for , fill_categ = "mutation_info_labels" , fill_categ_cols = c("#E69F00", "#999999") , my_ats = 15 # axis text size @@ -30,46 +30,44 @@ lineage_distP <- function(plotdf , leg_label = "") { - -if(!all_lineages){ - plotdf = plotdf[plotdf[[y_axis]]%in%use_lineages,] -} - -LinDistP = ggplot(plotdf, aes_string(x = x_axis - , y = y_axis))+ - - geom_density_ridges(aes_string(fill = fill_categ) - , scale = 3 - , size = 0.3 - , alpha = 0.8) + - scale_x_continuous(expand = c(0.01, 0.01)) + - #coord_cartesian( xlim = c(-1, 1)) + - scale_fill_manual(values = fill_categ_cols) + - theme(axis.text.x = element_text(size = my_ats - , angle = 90 - , hjust = 1 - , vjust = 0.4) - , axis.text.y = element_text(size = my_ats) - , axis.title.x = element_text(size = my_ats) - , axis.title.y = element_blank() - , strip.text = element_text(size = my_strip_ts) - , legend.text = element_text(size = my_leg_ts) - , legend.title = element_text(size = my_leg_title) - , legend.position = c(0.8, 0.9)) + - labs(x = x_lab - , fill = leg_label) - -if (with_facet){ - - # used reformulate or make as formula - #fwv = reformulate(facet_wrap_var) - fwv = as.formula(paste0("~", facet_wrap_var)) - LinDistP = LinDistP + + if(!all_lineages){ + plotdf = plotdf[plotdf[[y_axis]]%in%use_lineages,] + } + + ggplot(plotdf, aes_string(x = x_axis + , y = y_axis))+ + + geom_density_ridges(aes_string(fill = fill_categ) + , scale = 3 + , size = 0.3 + , alpha = 0.8) + + scale_x_continuous(expand = c(0.01, 0.01)) + + #coord_cartesian( xlim = c(-1, 1)) + + scale_fill_manual(values = fill_categ_cols) + + theme(axis.text.x = element_text(size = my_ats + , angle = 90 + , hjust = 1 + , vjust = 0.4) + , axis.text.y = element_text(size = my_ats) + , axis.title.x = element_text(size = my_ats) + , axis.title.y = element_blank() + , strip.text = element_text(size = my_strip_ts) + , legend.text = element_text(size = my_leg_ts) + , legend.title = element_text(size = my_leg_title) + , legend.position = c(0.8, 0.9)) + + labs(x = x_lab + , fill = leg_label) + + + # FIXME: This didn't work BEFORE i fixed the ggplot() assignment thing!!! + if (with_facet){ + + # used reformulate or make as formula + #fwv = reformulate(facet_wrap_var) + fwv = as.formula(paste0("~", facet_wrap_var)) + facet_wrap(fwv) + - theme(legend.position = leg_pos_wf - , legend.direction = leg_dir_wf) -} - -return(LinDistP) + theme(legend.position = leg_pos_wf + , legend.direction = leg_dir_wf) + } } diff --git a/scripts/functions/logoP_or.R b/scripts/functions/logoP_or.R index 8db127f..650bf8a 100644 --- a/scripts/functions/logoP_or.R +++ b/scripts/functions/logoP_or.R @@ -1,5 +1,5 @@ # Input: - # Data: +# Data: # plot_df: merged_df3 containing the OR column to use as y-axis or any other relevant column # x_axis_colname = "position" @@ -7,19 +7,19 @@ # symbol_colname = "mutant_type" # y_axis_log = F # log_value = log10 - # if used, y-axis label has "Log" appended to it +# if used, y-axis label has "Log" appended to it # my_logo_col = c("chemistry", "hydrophobicity", "clustalx", "taylor") - # --> if clustalx and taylor, set variable to black bg + white font - # --> if chemistry and hydrophobicity, then grey bg + black font +# --> if clustalx and taylor, set variable to black bg + white font +# --> if chemistry and hydrophobicity, then grey bg + black font # rm_empty_y = F - # option to remove empty positions i.e positions with no assocaited y-val +# option to remove empty positions i.e positions with no assocaited y-val # y_axis_log = F - # option to use log scale - # FIXME Minor bug: if used with rm_empty_y, sometimes the labels are too small to render(!?) - # so positions appear empty despite having y-vals +# option to use log scale +# FIXME Minor bug: if used with rm_empty_y, sometimes the labels are too small to render(!?) +# so positions appear empty despite having y-vals # ...other params @@ -39,48 +39,52 @@ # logo data: OR #================== LogoPlotCustomH <- function(plot_df - , x_axis_colname = "position" - , y_axis_colname = "or_mychisq" - , symbol_colname = "mutant_type" - , my_logo_col = "chemistry" - , rm_empty_y = F - , y_axis_log = F - , log_value = log10 - , y_axis_increment = 5 - , x_lab = "Position" - , y_lab = "Odds Ratio" - , x_ats = 12 # text size - , x_tangle = 90 # text angle - , y_ats = 22 - , y_tangle = 0 - , x_tts = 20 # title size - , y_tts = 23 - , leg_pos = "none" # can be top, left, right and bottom or c(0.8, 0.9) - , leg_dir = "horizontal" #can be vertical or horizontal - , leg_ts = 15 # leg text size - , leg_tts = 16 # leg title size - ) + , x_axis_colname = "position" + , y_axis_colname = "or_mychisq" + , symbol_colname = "mutant_type" + , my_logo_col = "chemistry" + , rm_empty_y = F + , y_axis_log = F + , log_value = log10 + , y_axis_increment = 5 + , x_lab = "Position" + , y_lab = "Odds Ratio" + , x_ats = 12 # text size + , x_tangle = 90 # text angle + , y_ats = 22 + , y_tangle = 0 + , x_tts = 20 # title size + , y_tts = 23 + , leg_pos = "none" # can be top, left, right and bottom or c(0.8, 0.9) + , leg_dir = "horizontal" #can be vertical or horizontal + , leg_ts = 15 # leg text size + , leg_tts = 16 # leg title size + , tpos0 = 0 # 0 is a magic number that does my sensible default + , tW0 = 1 + , tH0 = 0.3 +) { - ################################# - # Data processing for logo plot - ################################# - - if (rm_empty_y){ - plot_df = plot_df[!is.na(plot_df[y_axis_colname]),] - cat("\nRemoving empty positions...\n") - }else{ - plot_df = plot_df - } - - y_max = max(plot_df['or_mychisq'], na.rm = T) - cat("\nRemoving y scale incremenet:", y_axis_increment) - y_lim = round_any(y_max, y_axis_increment, f = ceiling) - - #------------------- - # logo data: LogOR - #------------------- - if (y_axis_log){ + ################################# + # Data processing for logo plot + ################################# + plot_df = generate_distance_colour_map(plot_df, yvar_colname = y_axis_colname, debug=TRUE) + + if (rm_empty_y){ + plot_df = plot_df[!is.na(plot_df[y_axis_colname]),] + cat("\nRemoving empty positions...\n") + }else{ + plot_df = plot_df + } + + y_max = max(plot_df['or_mychisq'], na.rm = T) + cat("\nRemoving y scale incremenet:", y_axis_increment) + y_lim = round_any(y_max, y_axis_increment, f = ceiling) + + #------------------- + # logo data: LogOR + #------------------- + if (y_axis_log){ log_colname = paste0("log_", y_axis_colname) plot_df[log_colname] = log_value(plot_df[y_axis_colname]) @@ -95,8 +99,8 @@ LogoPlotCustomH <- function(plot_df #y_lim = round_any(y_max, y_axis_increment, f = ceiling) - } else { - + } else { + #------------------- # logo data: OR #------------------- @@ -105,7 +109,7 @@ LogoPlotCustomH <- function(plot_df logo_dfP_wf = as.matrix(logo_df_plot %>% spread(x_axis_colname, y_axis_colname, fill = 0.0)) } - + class(logo_dfP_wf) rownames(logo_dfP_wf) = logo_dfP_wf[,1] @@ -116,7 +120,7 @@ LogoPlotCustomH <- function(plot_df colnames(logo_dfP_wf) position_or = as.numeric(colnames(logo_dfP_wf)) - + ###################################### # Generating plots with given y_axis ##################################### @@ -130,7 +134,7 @@ LogoPlotCustomH <- function(plot_df xtt_col = "white" ytt_col = "white" } - + if (my_logo_col %in% c('chemistry', 'hydrophobicity')) { cat('\nSelected colour scheme:', my_logo_col , "\nUsing grey theme") @@ -141,57 +145,60 @@ LogoPlotCustomH <- function(plot_df xtt_col = "black" ytt_col = "black" } - - p0 = ggseqlogo(logo_dfP_wf - , method = "custom" - , col_scheme = my_logo_col - , seq_type = "aa") + + + if (y_axis_log){ + + if (grepl("Log", y_lab)){ + y_lab = y_lab + + }else{ + y_lab = paste("Log", y_lab) + } + } + ggseqlogo(logo_dfP_wf + , method = "custom" + , col_scheme = my_logo_col + , seq_type = "aa") + #ylab("my custom height") + theme(axis.text.x = element_text(size = x_ats - , angle = x_tangle - , hjust = 1 - , vjust = 0.4 - , colour = xfont_bgc) - , axis.text.y = element_text(size = y_ats - , angle = y_tangle - , hjust = 1 - , vjust = 0 - , colour = yfont_bgc) - , axis.title.x = element_text(size = x_tts + , angle = x_tangle + , hjust = 1 + , vjust = 0.4 + , colour = xfont_bgc) + , axis.text.y = element_text(size = y_ats + , angle = y_tangle + , hjust = 1 + , vjust = 0 + , colour = yfont_bgc) + , axis.title.x = element_text(size = x_tts , colour = xtt_col) - , axis.title.y = element_text(size = y_tts + , axis.title.y = element_text(size = y_tts , colour = ytt_col) - , legend.title = element_text(size = leg_tts + , legend.title = element_text(size = leg_tts , colour = ytt_col) - , legend.text = element_text(size = leg_ts) + , legend.text = element_text(size = leg_ts) - , legend.position = leg_pos - , legend.direction = leg_dir - , plot.background = element_rect(fill = theme_bgc))+ + , legend.position = leg_pos + , legend.direction = leg_dir + , plot.background = element_rect(fill = theme_bgc))+ scale_x_discrete(x_lab #, breaks , labels = position_or - , limits = factor(1:length(position_or))) - - LogoPlot = p0 + scale_y_continuous(y_lab + , limits = factor(1:length(position_or))) + + + scale_y_continuous(y_lab , breaks = seq(0, (y_lim), by = y_axis_increment) #, labels = seq(0, (y_lim), by = y_axis_increment) - , limits = c(0, y_lim)) - if (y_axis_log){ - - if (grepl("Log", y_lab)){ - y_lab = y_lab - - }else{ - y_lab = paste("Log", y_lab) - } - - LogoPlot = p0 + ylab(y_lab) + , limits = c(0, y_lim)) + + ylab(y_lab) + + geom_tile(aes(plot_df$position, tpos0 # heat-mapped distance tiles along the bot + , width = tW0 + , height = tH0) + , fill = plot_df$ligD_colours + , colour = plot_df$ligD_colours + , linetype = "blank") + #LogoPlot = p0 + ylab(y_lab) + #return(LogoPlot) - } - - - return(LogoPlot) - } diff --git a/scripts/functions/logoP_snp.R b/scripts/functions/logoP_snp.R index 984aaad..d48df52 100644 --- a/scripts/functions/logoP_snp.R +++ b/scripts/functions/logoP_snp.R @@ -1,7 +1,7 @@ ########################a########################################################### # Input: - # Data - # plot_df: merged_df3 containing the OR column to use as y-axis or any other relevant column +# Data +# plot_df: merged_df3 containing the OR column to use as y-axis or any other relevant column # x_axis_colname = "position" # symbol_mut_colname = "mutant_type" @@ -15,7 +15,7 @@ # ...other params # Returns: Logo plot from combined data containing all nsSNPs per position. - # Helps to see the overview of SNP diversity +# Helps to see the overview of SNP diversity # TODO: SHINY # select/drop down: omit_snp_count @@ -31,38 +31,47 @@ # NOTE: my_logo_col LogoPlotSnps <- function(plot_df - , x_axis_colname = "position" - , symbol_mut_colname = "mutant_type" - , symbol_wt_colname = "mutant_type" - , omit_snp_count = c(0) # can be 1, 2, etc. - , my_logo_col = "chemistry" - , x_lab = "Position" - , y_lab = "Count" - , x_ats = 14 # text size - , x_tangle = 90 # text angle - , y_ats = 22 - , y_tangle = 0 - , x_tts = 20 # title size - , y_tts = 23 - , leg_pos = "none" # can be top, left, right and bottom or c(0.8, 0.9) - , leg_dir = "horizontal" #can be vertical or horizontal - , leg_ts = 20 # leg text size - , leg_tts = 16 # leg title size - ) + , x_axis_colname = "position" + , symbol_mut_colname = "mutant_type" + , symbol_wt_colname = "mutant_type" + , omit_snp_count = c(0) # can be 1, 2, etc. + , my_logo_col = "chemistry" + , x_lab = "Position" + , y_lab = "Count" + , x_ats = 14 # text size + , x_tangle = 90 # text angle + , y_ats = 22 + , y_tangle = 0 + , x_tts = 20 # title size + , y_tts = 23 + , leg_pos = "none" # can be top, left, right and bottom or c(0.8, 0.9) + , leg_dir = "horizontal" #can be vertical or horizontal + , leg_ts = 20 # leg text size + , leg_tts = 16 # leg title size + , tpos0 = 0 # 0 is a magic number that does my sensible default + , tW0 = 1 + , tH0 = 0.2 + , debug=FALSE + +) { # handle funky omit_snp_count. DOES NOT WORK YET if (class(omit_snp_count) != "numeric"){ omit_snp_count <- as.numeric(unlist(str_extract_all(omit_snp_count, regex("[0-9]+")))) } - ############################################ - # Data processing for logo plot for nsSNPS - ############################################ - setDT(plot_df)[, mut_pos_occurrence := .N, by = .(eval(parse(text=x_axis_colname)))] + ############################################ + # Data processing for logo plot for nsSNPS + ############################################ + # Generate "ligand distance" colour map + plot_df = generate_distance_colour_map(plot_df, debug=TRUE) + + setDT(plot_df)[, mut_pos_occurrence := .N, by = .(eval(parse(text=x_axis_colname)))] + if (debug) { table(plot_df[[x_axis_colname]]) table(plot_df$mut_pos_occurrence) - + } max_mut = max(table(plot_df[[x_axis_colname]])) # Subset Data as specified by user @@ -73,13 +82,14 @@ LogoPlotSnps <- function(plot_df my_data_snp = plot_df u = unique(my_data_snp[[x_axis_colname]]) max_mult_mut = max(table(my_data_snp[[x_axis_colname]])) - - cat("\nNo filtering requested:" - , "\nTotal no. of nsSNPs:", sum(table(plot_df$mut_pos_occurrence)) - , "\nTotal no. of nsSNPs omitted:", sum(table(plot_df$mut_pos_occurrence)[omit_snp_count]) - , "\nDim of data:", dim(my_data_snp) - , "\nNo. of positions:", length(u) - , "\nMax no. of muts at any position:", max_mult_mut) + if (debug) { + cat("\nNo filtering requested:" + , "\nTotal no. of nsSNPs:", sum(table(plot_df$mut_pos_occurrence)) + , "\nTotal no. of nsSNPs omitted:", sum(table(plot_df$mut_pos_occurrence)[omit_snp_count]) + , "\nDim of data:", dim(my_data_snp) + , "\nNo. of positions:", length(u) + , "\nMax no. of muts at any position:", max_mult_mut) + } } else { my_data_snp = subset(plot_df, !(mut_pos_occurrence%in%omit_snp_count) ) @@ -88,19 +98,20 @@ LogoPlotSnps <- function(plot_df got_rows = sum(table(my_data_snp$mut_pos_occurrence)) u = unique(my_data_snp[[x_axis_colname]]) max_mult_mut = max(table(my_data_snp[[x_axis_colname]])) - - if (got_rows == exp_nrows) { - cat("\nPass: Position with the stated nsSNP frequency filtered:", omit_snp_count - , "\nTotal no. of nsSNPs:", sum(table(plot_df$mut_pos_occurrence)) - , "\nTotal no. of nsSNPs omitted:", sum(table(plot_df$mut_pos_occurrence)[omit_snp_count]) - , "\nDim of subsetted data:", dim(my_data_snp) - , "\nNo. of positions:", length(u) - , "\nMax no. of muts at any position:", max_mult_mut) - } else { - - cat("\nFAIL:Position with the stated nsSNP frequency COULD NOT be filtered..." - , "\nExpected:",exp_nrows - , "\nGot:", got_rows ) + if (debug) { + if (got_rows == exp_nrows) { + cat("\nPass: Position with the stated nsSNP frequency filtered:", omit_snp_count + , "\nTotal no. of nsSNPs:", sum(table(plot_df$mut_pos_occurrence)) + , "\nTotal no. of nsSNPs omitted:", sum(table(plot_df$mut_pos_occurrence)[omit_snp_count]) + , "\nDim of subsetted data:", dim(my_data_snp) + , "\nNo. of positions:", length(u) + , "\nMax no. of muts at any position:", max_mult_mut) + } else { + + cat("\nFAIL:Position with the stated nsSNP frequency COULD NOT be filtered..." + , "\nExpected:",exp_nrows + , "\nGot:", got_rows ) + } } } @@ -116,17 +127,21 @@ LogoPlotSnps <- function(plot_df tab_mt = unclass(tab_mt) if (is.matrix(tab_mt)){ - cat("\nPASS: Mutant matrix successfully created..." - #, "\nRownames of mutant matrix:", rownames(tab_mt) - #, "\nColnames of mutant matrix:", colnames(tab_mt) - ) + if (debug) { + cat("\nPASS: Mutant matrix successfully created..." + #, "\nRownames of mutant matrix:", rownames(tab_mt) + #, "\nColnames of mutant matrix:", colnames(tab_mt) + ) + } } else{ tab_mt = as.matrix(tab_mt, rownames = T) if (is.matrix(tab_mt)){ - cat("\nCreating mutant matrix..." - #, "\nRownames of mutant matrix:", rownames(tab_mt) - #, "\nColnames of mutant matrix:", colnames(tab_mt) - ) + if (debug) { + cat("\nCreating mutant matrix..." + #, "\nRownames of mutant matrix:", rownames(tab_mt) + #, "\nColnames of mutant matrix:", colnames(tab_mt) + ) + } } } @@ -146,18 +161,19 @@ LogoPlotSnps <- function(plot_df wt tab_wt = table(wt[[symbol_wt_colname]], wt[[x_axis_colname]]); tab_wt # should all be 1 - - if ( identical(colnames(tab_mt), colnames(tab_wt) ) && identical(ncol(tab_mt), ncol(tab_wt)) ){ - - cat("\nPASS: Wild type matrix successfully created" - , "\nDim of wt matrix:", dim(tab_wt) - , "\nDim of mutant matrix:", dim(tab_mt) - , "\n" - #, "\nRownames of mutant matrix:", rownames(tab_wt) - #, "\nColnames of mutant matrix:", colnames(tab_wt) - ) + if (debug) { + if ( identical(colnames(tab_mt), colnames(tab_wt) ) && identical(ncol(tab_mt), ncol(tab_wt)) ){ + + cat("\nPASS: Wild type matrix successfully created" + , "\nDim of wt matrix:", dim(tab_wt) + , "\nDim of mutant matrix:", dim(tab_mt) + , "\n" + #, "\nRownames of mutant matrix:", rownames(tab_wt) + #, "\nColnames of mutant matrix:", colnames(tab_wt) + ) + } } - + ###################################### # Generating plots for muts and wt ##################################### @@ -184,132 +200,100 @@ LogoPlotSnps <- function(plot_df xtt_col = "black" ytt_col = "black" } - + ##################################### # Generating logo plots for nsSNPs ##################################### - - #------------------- - # Mutant logo plot - #------------------- - p0 = ggseqlogo(tab_mt - , method = 'custom' - , col_scheme = my_logo_col - , seq_type = 'aa') + - - theme(text=element_text(family="FreeSans"))+ - theme(axis.text.x = element_blank()) + - theme_logo()+ - - scale_x_continuous(breaks = 1:ncol(tab_mt) - , expand = c(0.01,0) - , labels = colnames(tab_mt))+ - scale_y_continuous(breaks = 0:(max_mult_mut-1) - , labels = c(1:max_mult_mut) - , limits = c(0, max_mult_mut)) + - #xlab(x_lab) + - ylab(y_lab) - - cat('\nDone: p0') - - # further customisation - mut_logo_p = p0 + theme(legend.position = leg_pos - , legend.direction = leg_dir - #, legend.title = element_blank() - , legend.title = element_text(size = leg_tts - , colour = ytt_col) - , legend.text = element_text(size = leg_ts) - - , axis.text.x = element_text(size = x_ats - , angle = x_tangle - , hjust = 1 - , vjust = 0.4 - , colour = xfont_bgc) - #, axis.text.y = element_blank() - , axis.text.y = element_text(size = y_ats - , angle = y_tangle - , hjust = 1 - , vjust = -1.0 - , colour = yfont_bgc) - , axis.title.x = element_text(size = x_tts - , colour = xtt_col) - , axis.title.y = element_text(size = y_tts - , colour = ytt_col) - - , plot.background = element_rect(fill = theme_bgc)) - - cat('\nDone: mut_logo_p') - #return(mut_logo_p) - LogoPlotL[['mut_logoP']] <- mut_logo_p - + cowplot::plot_grid( + #------------------- + # Mutant logo plot + #------------------- + ggseqlogo(tab_mt + , method = 'custom' + , col_scheme = my_logo_col + , seq_type = 'aa') + + + theme(text=element_text(family="FreeSans"))+ + theme(axis.text.x = element_blank()) + + theme_logo()+ + + scale_x_continuous(breaks = 1:ncol(tab_mt) + , expand = c(0.01,0) + , labels = colnames(tab_mt))+ + + scale_y_continuous(breaks = 0:(max_mult_mut-1) + , labels = c(1:max_mult_mut) + , limits = c(0, max_mult_mut)) + + # FIXME: currently broken, possibly due to ggseqlogo() not working in the + # way standard ggplot2() does + geom_tile(aes(plot_df$position, tpos0 # heat-mapped distance tiles along the bottom. + , width = tW0 + , height = tH0) + , fill = plot_df$ligD_colours + , colour = plot_df$ligD_colours + , linetype = "blank") + + ylab(y_lab) + + theme(legend.position = leg_pos + , legend.direction = leg_dir + , legend.title = element_text(size = leg_tts + , colour = ytt_col) + , legend.text = element_text(size = leg_ts) + + , axis.text.x = element_text(size = x_ats + , angle = x_tangle + , hjust = 1 + , vjust = 0.4 + , colour = xfont_bgc) + , axis.text.y = element_text(size = y_ats + , angle = y_tangle + , hjust = 1 + , vjust = -1.0 + , colour = yfont_bgc) + , axis.title.x = element_text(size = x_tts + , colour = xtt_col) + , axis.title.y = element_text(size = y_tts + , colour = ytt_col) + + , plot.background = element_rect(fill = theme_bgc) + ) + ,ggseqlogo(tab_wt + , method = 'custom' + , col_scheme = my_logo_col + , seq_type = 'aa') + + + theme(text = element_text(family="FreeSans"))+ + theme(axis.text.x = element_blank() + , axis.text.y = element_blank()) + + theme_logo()+ + + scale_x_continuous(breaks = 1:ncol(tab_wt) + , expand = c(0.01,0) + , labels = colnames(tab_wt))+ + + xlab(x_lab) + + + theme(legend.position = "none" + , legend.direction = leg_dir + #, legend.title = element_blank() + , legend.title = element_text(size = y_tts + , colour = ytt_col) + , legend.text = element_text(size = leg_ts) + , axis.text.x = element_text(size = x_ats + , angle = x_tangle + , hjust = 1 + , vjust = 0.4 + , colour = xfont_bgc) + , axis.text.y = element_blank() + , axis.title.x = element_text(size = x_tts + , colour = xtt_col) + , axis.title.y = element_text(size = y_tts + , colour = ytt_col) + , plot.background = element_rect(fill = theme_bgc) + ) + , nrow = 2 + , align = "v" + , rel_heights = c(3/4, 1/4)) #------------------ # Wild logo plot #------------------ - p1 = ggseqlogo(tab_wt - , method = 'custom' - , col_scheme = my_logo_col - , seq_type = 'aa') + - - theme(text = element_text(family="FreeSans"))+ - theme(axis.text.x = element_blank() - , axis.text.y = element_blank()) + - theme_logo()+ - - scale_x_continuous(breaks = 1:ncol(tab_wt) - , expand = c(0.01,0) - , labels = colnames(tab_wt))+ - - xlab(x_lab) - - cat('\nDone: p1') - - # further customisation - wt_logo_p = p1 + - - theme(legend.position = "none" - , legend.direction = leg_dir - #, legend.title = element_blank() - , legend.title = element_text(size = y_tts - , colour = ytt_col) - , legend.text = element_text(size = leg_ts) - - , axis.text.x = element_text(size = x_ats - , angle = x_tangle - , hjust = 1 - , vjust = 0.4 - , colour = xfont_bgc) - , axis.text.y = element_blank() - - , axis.title.x = element_text(size = x_tts - , colour = xtt_col) - , axis.title.y = element_text(size = y_tts - , colour = ytt_col) - - , plot.background = element_rect(fill = theme_bgc)) - - cat('\nDone: wt_logo_p') - #return(wt_logo_p) - LogoPlotL[['wt_logoP']] <- wt_logo_p - - #========================================= - # Output - # Combined plot: logo_mutliple_muts.svg - #========================================= - #suppressMessages( require(cowplot) ) - #plot_grid(p1, p3, ncol = 1, align = 'v') - cat('\nDone: mut_logo_p + wt_logo_p') - - # colour scheme: https://rdrr.io/cran/ggseqlogo/src/R/col_schemes.r - #cat("\nOutput plot:", LogoSNPs_comb, "\n") - #svg(LogoSNPs_combined, width = 32, height = 10) - LogoPlotL[['wt_logoP']] <- wt_logo_p - - LogoSNPs_comb = cowplot::plot_grid(LogoPlotL[['mut_logoP']] - , LogoPlotL[['wt_logoP']] - , nrow = 2 - , align = "v" - , rel_heights = c(3/4, 1/4)) - - return(LogoSNPs_comb) - } diff --git a/scripts/functions/position_count_bp.R b/scripts/functions/position_count_bp.R index 2907b4b..c08f2d5 100755 --- a/scripts/functions/position_count_bp.R +++ b/scripts/functions/position_count_bp.R @@ -105,9 +105,9 @@ site_snp_count_bp <- function (plotdf # not sure if to use with sort or directly my_x = sort(unique(snpsBYpos_df$snpsBYpos)) - g = ggplot(snpsBYpos_df, aes(x = snpsBYpos)) - OutPlot_pos_count = g + geom_bar(aes (alpha = 0.5) - , show.legend = FALSE) + + ggplot(snpsBYpos_df, aes(x = snpsBYpos)) + + geom_bar(aes (alpha = 0.5) + , show.legend = FALSE) + scale_x_continuous(breaks = unique(snpsBYpos_df$snpsBYpos)) + geom_label(stat = "count", aes(label = ..count..) , color = "black" @@ -126,14 +126,11 @@ site_snp_count_bp <- function (plotdf , colour = title_colour) , plot.subtitle = element_text(size = subtitle_size , hjust = 0.5 - , colour = subtitle_colour)) + - + , colour = subtitle_colour)) + labs(title = bp_plot_title , subtitle = subtitle_text , x = xaxis_title , y = yaxis_title) - - return(OutPlot_pos_count) } ######################################################################## diff --git a/scripts/functions/wideP_consurf.R b/scripts/functions/wideP_consurf.R new file mode 100644 index 0000000..e32f446 --- /dev/null +++ b/scripts/functions/wideP_consurf.R @@ -0,0 +1,506 @@ +#!/usr/bin/env Rscript + +######################################################### +# TASK: function for wide plot +#with consurf score and error bars +#position numbers coloured by +# - ligand distance +# - active site residues +######################################################### + +#========================================================== +# wideP(): +# input args +#========================================================== +wideP_consurf <- function(plot_df + , xvar_colname = "position" + , yvar_colname = "consurf_score" + , yvar_colourN_colname = "consurf_colour_rev" # num from 0-1 + , plot_error_bars = T + , upper_EB_colname = "consurf_ci_upper" + , lower_EB_colname = "consurf_ci_lower" + + , plot_type = "point" # default is point + , point_colours + , p_size = 2 + , leg_title1 = "" + , leg_labels = c("0": "Insufficient Data" + , "1": "Variable" + , "2", "3", "4", "5", "6", "7", "8" + , "9": "Conserved") + , panel_col = "black" + , panel_col_fill = "black" + + # axes title and label sizes + , x_axls = 12 # x-axis label size + , y_axls = 15 # y-axis label size + , x_axts = 12 # x-axis text size + , y_axts = 12 # y-axis text size + , default_xtc = "black" # x-axis text colour + , ptitle = "" + , xlab = "" + , ylab = "" + , pts = 20 + + # plot margins + , t_margin = 0.5 + , r_margin = 0.5 + , b_margin = 1 + , l_margin = 1 + , unit_margin = "cm" + + # Custom 1: x-axis: text colour + , xtext_colour_aa = F + , xtext_colour_aa1 = active_aa_pos + , xtext_colour_aa2 = aa_pos_drug + , xtext_colours = c("purple", "brown", "black") + + # Custom 2: x-axis: geom tiles ~ lig distance + , A_xvar_lig = T + , leg_title2 = "Ligand Distance" + , lig_dist_colname = LigDist_colname # from globals + , lig_dist_colours = c("green", "yellow", "orange", "red") + , tpos0 = 0 # 0 is a magic number that does my sensible default + , tW0 = 1 + , tH0 = 0.3 + + # Custom 3: x-axis: geom tiles ~ active sites and ligand + , A_xvar_aa = F + , aa_pos_drug = NULL + , drug_aa_colour = "purple" + , tW = 1 + , tH = 0.2 + , active_aa_pos = NULL + , active_aa_colour = "brown" + + , aa_pos_lig1 = NULL + , aa_colour_lig1 = "blue" + , tpos1 = 0 + + , aa_pos_lig2 = NULL + , aa_colour_lig2 = "cyan" + , tpos2 = 0 + + , aa_pos_lig3 = NULL + , aa_colour_lig3 = "cornflowerblue" + , tpos3 = 0 + + , default_gt_clr = "white" + , debug=FALSE +){ + + if(missing(point_colours)){ + temp_cols = colorRampPalette(c("seagreen", "sienna3"))(30) + point_colours = temp_cols + }else{ + point_colours = point_colours + } + + ############################### + # custom 1: x-axis text colour + ############################## + + if (xtext_colour_aa){ + positionF <- levels(as.factor(plot_df[[xvar_colname]])) + length(positionF) + aa_pos_colours = ifelse(positionF%in%xtext_colour_aa1, xtext_colours[1] + , ifelse(positionF%in%xtext_colour_aa2 + , xtext_colours[2] + , xtext_colours[3])) + }else{ + aa_pos_colours = default_xtc + } + + ################################################ + # Custom 2: x-axis geom tiles ~ lig distance + ################################################ + + #========================= + # Build data with colours + # ~ ligand distance + #========================= + if (A_xvar_lig){ + cat("\nAnnotating x-axis ~", lig_dist_colname, "requested...") + + #------------------------------------- + # round column values: to colour by + #-------------------------------------- + #plot_df = plot_df[order(plot_df[[lig_dist_colname]]),] + plot_df['lig_distR'] = round(plot_df[[lig_dist_colname]], digits = 0) + head(plot_df['lig_distR']) + + #------------------------------------- + # ligand distance range, min, max, etc + #-------------------------------------- + lig_min = min(round(plot_df[[lig_dist_colname]]), na.rm = T); lig_min + lig_max = max(round(plot_df[[lig_dist_colname]]), na.rm = T); lig_max + lig_mean = round(mean(round(plot_df[[lig_dist_colname]]), na.rm = T)); lig_mean + + #------------------------------------- + # Create mapping colour key + #-------------------------------------- + # sorting removes NA, so that n_colours == length(ligD_valsR) + n_colours = length(sort(unique(round(plot_df[[lig_dist_colname]], digits = 0)))); n_colours + + lig_cols = colorRampPalette(lig_dist_colours)(n_colours); lig_cols + ligD_valsR = sort(unique(round(plot_df[[lig_dist_colname]], digits = 0))); ligD_valsR + length(ligD_valsR) + + if (n_colours == length(ligD_valsR)) { + cat("\nStarting: mapping b/w" + , lig_dist_colname + , "and colours") + }else{ + cat("\nCannot start mapping b/w", lig_dist_colname, "and colours..." + , "\nLength mismatch:" + , "No. of colours: ", n_colours + , "\nValues to map:", length(ligD_valsR)) + } + + ligDcolKey <- data.frame(ligD_colours = lig_cols + , lig_distR = ligD_valsR); ligDcolKey + names(ligDcolKey) + cat("\nSuccessful: Mapping b/w", lig_dist_colname, "and colours") + + #------------------------------------- + # merge colour key with plot_df + #-------------------------------------- + plot_df = merge(plot_df, ligDcolKey, by = 'lig_distR') + + plot_df_check = as.data.frame(cbind(position = plot_df[[xvar_colname]] + , ligD = plot_df[[lig_dist_colname]] + , ligDR = plot_df$lig_distR + , ligD_cols = plot_df$ligD_colours)) + } else{ + plot_df = plot_df + } + + ############################################### + # Custom 3: x-axis geom tiles ~ active sites + ################################################ + + #========================== + # Build Data with colours + # ~ on active sites + #========================== + + if(A_xvar_aa) { + cat("\nAnnotation for xvar requested. Building colours for annotation...") + + aa_colour_colname = "bg_all" + aa_colour_colname1 = "col_bg1" + aa_colour_colname2 = "col_bg2" + aa_colour_colname3 = "col_bg3" + + #-------------------------------------------------- + # column colour 0: Active site + drug binding sites + #-------------------------------------------------- + plot_df[[aa_colour_colname]] = ifelse(plot_df[[xvar_colname]]%in%aa_pos_drug + , drug_aa_colour + , ifelse(plot_df[[xvar_colname]]%in%active_aa_pos + , active_aa_colour, default_gt_clr )) + plot_df[[aa_colour_colname]] + cat("\nColumn created 'bg_all':", length(plot_df[[aa_colour_colname]])) + + #------------------------------------------------ + # column colour 1: Ligand 1 + drug binding sites + #------------------------------------------------ + cat("\nAssigning colours to drug binding and ligand-1 binding residues") + plot_df[[aa_colour_colname1]] = plot_df[[aa_colour_colname]] + plot_df[[aa_colour_colname1]] = ifelse(plot_df[[xvar_colname]]%in%aa_pos_lig1 + , aa_colour_lig1, plot_df[[aa_colour_colname]]) + #------------------------------------------------ + # column colour 2: Ligand 2 + #------------------------------------------------ + plot_df[[aa_colour_colname2]] = plot_df[[aa_colour_colname1]] + plot_df[[aa_colour_colname2]] = ifelse(plot_df[[xvar_colname]]%in%aa_pos_lig2 + , aa_colour_lig2, plot_df[[aa_colour_colname1]]) + + #------------------------------------------------ + # column colour 3: Ligand 3 + #------------------------------------------------ + plot_df[[aa_colour_colname3]] = plot_df[[aa_colour_colname2]] + plot_df[[aa_colour_colname3]] = ifelse(plot_df[[xvar_colname]]%in%aa_pos_lig3 + , aa_colour_lig3, plot_df[[aa_colour_colname2]]) + + } + ################### + # start plot + ################### + + #------------------- + # x and y axis + # range, scale, etc + #------------------- + my_xlim = length(unique(plot_df[[xvar_colname]])); my_xlim + ymin = min(plot_df[[yvar_colname]]); ymin + ymax = max(plot_df[[yvar_colname]]); ymax + + g = ggplot(plot_df, aes_string(x = sprintf("factor(%s)", xvar_colname) + , y = yvar_colname + , colour = sprintf("factor(%s)", yvar_colourN_colname) + )) + + "if SPECIAL do SPECIAL THING, otherwise do NORMAL THING" + if (plot_type == "bar"){ + g0 = g + + geom_bar(stat = "identity") + } + else{ + g0 = g + + coord_cartesian(xlim = c(1, my_xlim) + , ylim = c(ymin, ymax) + , clip = "off") + + geom_point(size = p_size) + + scale_colour_manual(values = point_colours) + } + + if (plot_error_bars){ + g0 = g0 + + geom_errorbar(aes(ymin = eval(parse(text = lower_EB_colname)) + , ymax = eval(parse(text = upper_EB_colname)) + )) + }else{ + + g0 = g0 + + } + + #--------------------- + # add axis formatting + #--------------------- + g1 = g0 + theme( axis.text.x = element_text(size = x_axts + , angle = 90 + , hjust = 1 + , vjust = 0.4 + , face = "bold" + , colour = aa_pos_colours) + , axis.text.y = element_text(size = y_axts + , angle = 0 + , hjust = 1 + , vjust = 0) + , axis.title.x = element_text(size = x_axls) + , axis.title.y = element_text(size = y_axls ) + , panel.background = element_rect(fill = panel_col_fill, color = panel_col) + , panel.grid.major = element_line(color = "black") + , panel.grid.minor = element_line(color = "black") + , plot.title = element_text(size = pts + , hjust = 0.5) + , plot.margin = margin(t = t_margin + , r = r_margin + , b = b_margin + , l = l_margin + , unit = unit_margin))+ + guides(colour = guide_legend(title = "ConsurfXXXX")) + + + labs(title = ptitle + , x = xlab + , y = ylab) + + #------------------ + #Extract legend1 + #------------------ + # yayy + g1_leg = ggplot(plot_df, aes_string(x = sprintf("factor(%s)" + , xvar_colname) )) + g1_leg = g1_leg + geom_bar(); g1_leg + g1_leg = g1_leg + geom_bar(aes_string(fill = sprintf("factor(%s)" + , yvar_colourN_colname))) + + g1_leg = g1_leg + scale_fill_manual(values = consurf_palette2 , name = leg_title1) + g1_leg + + legend1 = get_legend(g1_leg) + + + ##################################################### + #============================================ + # x-axis: geom_tiles ~ ligand distance + #============================================ + #------- + # plot + #------- + if(A_xvar_lig){ # 0 is a magic number that does my sensible default + if (tpos0 == 0){ + tpos0 = ymin-0.5 + } + if (tpos1 == 0){ + tpos1 = ymin-0.65 + } + if (tpos2 == 0){ + tpos2 = ymin-0.75 + } + if (tpos3 == 0){ + tpos3 = ymin-0.85 + } + + + cat("\nColouring x-axis aa based on", lig_dist_colname + , "\nNo. of colours:", n_colours) + + g2 = g1 + geom_tile(aes(, 2 +# g2 = g1 + geom_tile(aes(, tpos0 + + , width = tW0 + , height = tH0) + , fill = plot_df$ligD_colours + , colour = plot_df$ligD_colours + , linetype = "blank") + + #cat("Nrows of plot df", length(plot_df$ligD_colours)) + out = g2 + + }else{ + out = g1 + } + #============================================== + # x-axis: geom_tiles ~ active sites and others + #============================================== + if(A_xvar_aa){ + #tpos = 0 + #tW = 1 + #tH = 0.2 + + #--------------------- + # Add2plot: 3 ligands + #--------------------- + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + !is.null(aa_pos_lig1) && !is.null(aa_pos_lig2) && !is.null(aa_pos_lig3))) { + if (debug){ + cat("\n\nAnnotating xvar with active, drug binding, and Lig 1&2&3 sites") + cat("\nCreating column colours, column name:", aa_colour_colname3) + + cat("\nDoing Plot with 3 ligands") + } + out = out + geom_tile(aes(,tpos3 + , width = tW + , height = tH ) + , fill = plot_df[[aa_colour_colname3]] + , colour = plot_df[[aa_colour_colname3]] + , linetype = "solid") + + geom_tile(aes(, tpos2 + , width = tW + , height = tH ) + , fill = plot_df[[aa_colour_colname2]] + , colour = plot_df[[aa_colour_colname2]] + , linetype = "solid")+ + + geom_tile(aes(, tpos1 + , width = tW + , height = tH) + , fill = plot_df[[aa_colour_colname1]] + , colour = plot_df[[aa_colour_colname1]] + , linetype = "solid") + if (debug){ + cat("\nDone Plot with 3 ligands") + } + } + #--------------------- + # Add2plot: 2 ligands + #--------------------- + + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + !is.null(aa_pos_lig1) && !is.null(aa_pos_lig2) && is.null(aa_pos_lig3))) { + if (debug){ + cat("\n\nAnnotating xvar with active, drug binding, and Lig 1&2 sites") + cat("\nCreating column colours, column name:", aa_colour_colname2) + + cat("\nDoing Plot with 2 ligands") + } + out = out + + geom_tile(aes(, tpos2 + , width = tW + , height = tH) + , fill = plot_df[[aa_colour_colname2]] + , colour = plot_df[[aa_colour_colname2]] + , linetype = "solid")+ + geom_tile(aes(, tpos1 + , width = tW + , height = tH) + , fill = plot_df[[aa_colour_colname1]] + , colour = plot_df[[aa_colour_colname1]] + , linetype = "solid") + if (debug){ + cat("\nDone Plot with 2 ligands") + } + } + + #--------------------- + # Add2plot: 1 ligand + #--------------------- + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + !is.null(aa_pos_lig1) && is.null(aa_pos_lig2) && is.null(aa_pos_lig3))) { + if (debug){ + cat("\n\nAnnotating xvar with active, drug binding, and Lig 1 sites") + cat("\nCreating column colours, column name:", aa_colour_colname1) + + cat("\nDoing Plot with 1 ligands") + } + out = out + + geom_tile(aes(, tpos1 + , width = tW + , height = tH) + , fill = plot_df[[aa_colour_colname1]] + , colour = plot_df[[aa_colour_colname1]] + , linetype = "solid") + + cat("\nDone Plot with 1 ligand") + + } + + #----------------------------------- + # Add2plot:NO ligands + # No Ligs: Just drug and active site + # DEFAULT: A_xvar_aa == TRUE + #---------------------------------- + if (all(!is.null(active_aa_pos) && + !is.null(aa_pos_drug) && + is.null(aa_pos_lig1) && + is.null(aa_pos_lig2) && + is.null(aa_pos_lig3))) { + if (debug){ + cat("\n\nAnnotating xvar with active and drug binding sites") + cat("\nCreating column colours, column name:", aa_colour_colname) + cat("\nDoing Plot with 0 ligands: active and drug site only") + } + out = out + geom_tile(aes(, tpos3 + , width = tW + , height = tH) + , fill = plot_df[[aa_colour_colname]] + , colour = plot_df[[aa_colour_colname]] + , linetype = "solid") + if (debug){ + cat("\nDone Plot with: Active and Drug sites") + } + } + }else{ + cat("\nNo annotation for additional ligands on xvar requested") + } + if (A_xvar_lig){ + legs = cowplot::plot_grid(legend1 + , generate_distance_legend(plot_df, yvar_colname = yvar_colname) + , ncol = 1 + , align = "hv" + , rel_heights = c(2/4,3/4)) + + out2 = cowplot::plot_grid( out + theme(legend.position = "none") + , legs + , ncol = 2 + , align = "hv" + , rel_widths = c(9/10, 0.4/10) + ) + }else{ + out2 = cowplot::plot_grid( out + theme(legend.position = "none") + , legend1 + , ncol = 2 + , align = "hv" + , rel_widths = c(9/10, 0.5/10) + ) + } + return(out2) + #return(out2) + +}