LSHTM_analysis/scripts/functions/wideP_consurf3.R

306 lines
12 KiB
R

#!/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_consurf3 <- 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 = "white"
, panel_col_fill = "white"
# 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
, r_margin = 0
, b_margin = 0
, l_margin = 0
, 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
, annotate_ligand_distance = T
, leg_title2 = "Ligand Distance"
, lig_dist_colname = LigDist_colname # from globals
, lig_dist_colours = c("tan", "black")
, 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
, annotate_active_sites = T
, default_gt_clr = "white"
, build_plot_df=FALSE
, debug=FALSE
, ...
){
anno_bar = position_annotation(
plot_df,
bg = panel_col,
...
)
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 (annotate_ligand_distance){
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))
}
###################
# 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"){
cat('\ndoing bar plot')
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, show.legend = FALSE) +
scale_colour_manual(values = point_colours)
}
if (plot_error_bars){
cat('\nplotting error bars')
g0 = g0 +
geom_errorbar(aes(ymin = eval(parse(text = lower_EB_colname))
, ymax = eval(parse(text = upper_EB_colname))
), show.legend = FALSE)
}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 = "transparent")
, plot.background = element_rect(fill = panel_col_fill, color = "transparent")
, panel.grid = element_blank()
#, 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(x = NULL, y = NULL)
# 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 = point_colours , name = leg_title1)
g1_leg
legend1 = get_legend(g1_leg)
out = g1
# x-axis: geom_tiles ~ active sites and others
if(annotate_active_sites){
#aligned = align_plots(out,position_annotation(plot_df))
out=cowplot::plot_grid(
out,
NULL,
anno_bar,
ncol=1,
align='v',
#axis='lr',
rel_heights = c(10,-0.1,1)
)
}
#### add the legends to the output ####
# ...which legends depend on whether the heat bar and position tiles are included
if (annotate_ligand_distance){
cat('\nOutput: Plot + distance heat-bar + other stuff')
legs = cowplot::plot_grid(legend1
, generate_distance_legend(plot_df)
, 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)
}
#wideP_consurf3(small_df3, point_colours = consurf_colours)