From 802d6f8495008b7dcf52998de548a51d766964a0 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Mon, 22 Aug 2022 11:41:42 +0100 Subject: [PATCH] renamed 2 to _v2 --- config/alr.R | 47 +- config/embb.R | 25 +- config/gid.R | 74 +- scripts/Header_TT.R | 28 +- scripts/functions/corr_plot_data.R | 39 +- scripts/functions/dm_om_data.R | 4 +- scripts/functions/plotting_globals.R | 2 +- scripts/plotting/get_plotting_dfs.R | 230 +++--- .../plotting_thesis/appendix_tables.R | 89 ++- .../plotting/plotting_thesis/basic_barplots.R | 753 +++++++----------- .../plotting_thesis/basic_barplots_layout.R | 198 ----- .../basic_barplots_layout_v2.R | 35 +- ...{basic_barplots2.R => basic_barplots_v2.R} | 176 ++-- .../plotting_thesis/corr_plots_thesis.R | 2 +- .../plotting/plotting_thesis/dm_om_plots.R | 8 +- .../plotting_thesis/dm_om_plots_layout.R | 9 +- .../plotting/plotting_thesis/gg_pairs_all.R | 14 +- .../plotting/plotting_thesis/linage_bp_dist.R | 4 +- 18 files changed, 761 insertions(+), 976 deletions(-) mode change 100755 => 100644 scripts/plotting/plotting_thesis/basic_barplots.R delete mode 100644 scripts/plotting/plotting_thesis/basic_barplots_layout.R rename scripts/plotting/plotting_thesis/{basic_barplots2.R => basic_barplots_v2.R} (79%) diff --git a/config/alr.R b/config/alr.R index bfbc350..3543b9a 100644 --- a/config/alr.R +++ b/config/alr.R @@ -26,7 +26,8 @@ aa_plip_dcs_hbond = c(66, 70, 196, 237 aa_plip_dcs_other = aa_plip_dcs[!aa_plip_dcs%in%aa_plip_dcs_hbond] c2 = length(aa_plip_dcs_other) == length(aa_plip_dcs) - length(aa_plip_dcs_hbond) - + + #========== # Arpeggio #=========== @@ -40,6 +41,18 @@ aa_arpeg_dcs_other = aa_arpeg_dcs[!aa_arpeg_dcs%in%c(aa_ligplus_dcs_other c3 = length(aa_arpeg_dcs_other) == length(aa_arpeg_dcs) - ( length(aa_ligplus_dcs_other) + length(aa_plip_dcs_other) ) +####################################################################### +#NEW AFTER ADDING PLP to structure! huh +# ADDED: 18 Aug 2022 +# PLIP server for co factor PLP (CONFUSING!) +#and 2019 lit:lys42, M319, and Y364 : OFFSET is 24 +#K42: K66, Y271:Y295, M319:M343, W89: W113, W203: W227, H209:H233, Q321:Q345 +aa_pos_paper = sort(unique(c(66,70,112,113,164,196,227,233,237,252,254,255,295,342,343,344,345,388))) +plp_pos_paper = sort(unique(c(66, 70, 112, 196, 227, 237, 252, 254, 255, 388))) + +#active_aa_pos = sort(unique(c(aa_pos_paper, active_aa_pos))) +aa_pos_plp = sort(unique(c(plp_pos_paper, 66, 70, 112, 237, 252, 254, 255, 196))) + ####################################################################### #=============== @@ -47,7 +60,9 @@ c3 = length(aa_arpeg_dcs_other) == length(aa_arpeg_dcs) - ( length(aa_ligplus_dc #=============== active_aa_pos = sort(unique(c(aa_ligplus_dcs , aa_plip_dcs - , aa_arpeg_dcs))) + , aa_arpeg_dcs + , aa_pos_paper + , aa_pos_plp))) #================= # Drug binding aa #================= @@ -56,6 +71,12 @@ aa_pos_dcs = sort(unique(c(aa_ligplus_dcs , aa_arpeg_dcs))) aa_pos_drug = aa_pos_dcs + +#=============== +# Co-factor: PLP aa +#=============== +aa_pos_plp = aa_pos_plp + #=============== # Hbond aa #=============== @@ -102,11 +123,25 @@ if ( all(c1, c2, c3) ) { , "\n\nNO other co-factors or ligands present\n") } -####################################################################### +###################################################################### +#NEW +# PLIP server for co factor PLP (CONFUSING!) +#and 2019 lit:lys42, M319, and Y364 : OFFSET is 24 +#K42: K66, Y271:Y295, M319:M343, W89: W113, W203: W227, H209:H233, Q321:Q345 +aa_pos_paper = sort(unique(c(66,70,112,113,164,196,227,233,237,252,254,255,295,342,343,344,345,388))) +plp_pos_paper = sort(unique(c(66, 70, 112, 196, 227, 237, 252, 254, 255, 388))) + +active_aa_pos = sort(unique(c(aa_pos_paper, active_aa_pos))) +aa_pos_plp = sort(unique(c(plp_pos_paper, 66, 70, 112, 237, 252, 254, 255, 196))) + +# add two key residues +#aa_pos_drug = sort(unique(c(319, 364, aa_pos_drug))) +#active_aa_pos = sort(unique(c(319, 364, active_aa_pos, aa_pos_plp))) + # FIXME: these should be populated! -aa_pos_lig1 = NULL +aa_pos_lig1 = aa_pos_plp aa_pos_lig2 = NULL aa_pos_lig3 = NULL -tile_map=data.frame(tile=c("ALR","DPA","CDL","Ca"), - tile_colour=c("green","darkslategrey","navyblue","purple")) +tile_map=data.frame(tile=c("ALR","PLP"), + tile_colour=c("green","darkslategrey")) \ No newline at end of file diff --git a/config/embb.R b/config/embb.R index f7a18e8..a6c3993 100644 --- a/config/embb.R +++ b/config/embb.R @@ -70,30 +70,30 @@ cat("\nNo. of active site residues for gene" ############################################################## aa_pos_emb = sort(unique(c( aa_ligplus_emb - , aa_plip_emb - , aa_arpeg_emb))) + , aa_plip_emb + , aa_arpeg_emb))) aa_pos_drug = aa_pos_emb aa_pos_emb_hbond = sort(unique(c( aa_ligplus_emb_hbond - , aa_plip_emb_hbond))) + , aa_plip_emb_hbond))) aa_pos_ca = sort(unique(c( aa_ligplus_ca - , aa_plip_ca - , aa_arpeg_ca))) + , aa_plip_ca + , aa_arpeg_ca))) aa_pos_cdl = sort(unique(c( aa_ligplus_cdl - , aa_plip_cdl - , aa_arpeg_cdl ))) - + , aa_plip_cdl + , aa_arpeg_cdl ))) + aa_pos_cdl_hbond = sort(unique(c( aa_ligplus_cdl_hbond ))) aa_pos_dsl = sort(unique(c( aa_ligplus_dsl - , aa_plip_dsl - , aa_arpeg_dsl))) + , aa_plip_dsl + , aa_arpeg_dsl))) aa_pos_dsl_hbond = sort(unique(c( aa_ligplus_dsl_hbond - , aa_plip_dsl_hbond))) - + , aa_plip_dsl_hbond))) + cat("\n===================================================" , "\nActive site residues for", gene, "comprise of..." @@ -120,3 +120,4 @@ aa_pos_lig3 = aa_pos_ca #purple tile_map=data.frame(tile=c("EMB","DPA","CDL","Ca"), tile_colour=c("green","darkslategrey","navyblue","purple")) +drug_main_res = c(299 , 302, 303 , 306 , 327 , 592 , 594, 988, 1028) diff --git a/config/gid.R b/config/gid.R index 34ff829..353607c 100644 --- a/config/gid.R +++ b/config/gid.R @@ -17,7 +17,7 @@ aa_ligplus_sam = c(148, 137, 138, 139 , 93, 69, 119, 120 , 220, 219, 118, 223) aa_ligplus_sam_hbond = c(220, 223) - + aa_ligplus_amp = c(123, 125, 213, 214) aa_ligplus_amp_hbond = c(125, 123, 213) @@ -53,19 +53,19 @@ aa_arpeg_amp = c(123, 125, 213) # Active site #============= active_aa_pos = sort(unique(c( -#rna_bind_aa_pos -#, binding_aa_pos -aa_ligplus_sry -, aa_ligplus_sam -, aa_ligplus_amp -, aa_ligplus_rna -, aa_plip_sry -, aa_plip_sam -, aa_plip_amp -, aa_plip_rna -, aa_arpeg_sry -, aa_arpeg_sam -, aa_arpeg_amp + #rna_bind_aa_pos + #, binding_aa_pos + aa_ligplus_sry + , aa_ligplus_sam + , aa_ligplus_amp + , aa_ligplus_rna + , aa_plip_sry + , aa_plip_sam + , aa_plip_amp + , aa_plip_rna + , aa_arpeg_sry + , aa_arpeg_sam + , aa_arpeg_amp ))) ############################################################## @@ -77,42 +77,42 @@ cat("\nNo. of active site residues for gene" ############################################################## aa_pos_sry = sort(unique(c( - aa_ligplus_sry - , aa_plip_sry - , aa_arpeg_sry))) + aa_ligplus_sry + , aa_plip_sry + , aa_arpeg_sry))) aa_pos_drug = aa_pos_sry aa_pos_sry_hbond = sort(unique(c( - aa_ligplus_sry_hbond - , aa_plip_sry_hbond))) + aa_ligplus_sry_hbond + , aa_plip_sry_hbond))) aa_pos_rna = sort(unique(c( - aa_ligplus_rna - , aa_plip_rna))) + aa_ligplus_rna + , aa_plip_rna))) aa_pos_rna_hbond = sort(unique(c( - aa_ligplus_rna_hbond - , aa_plip_rna_hbond))) + aa_ligplus_rna_hbond + , aa_plip_rna_hbond))) aa_pos_sam = sort(unique(c( - aa_ligplus_sam - , aa_plip_sam - , aa_arpeg_sam))) + aa_ligplus_sam + , aa_plip_sam + , aa_arpeg_sam))) aa_pos_sam_hbond = sort(unique(c( - aa_ligplus_sam_hbond - , aa_plip_sam_hbond))) + aa_ligplus_sam_hbond + , aa_plip_sam_hbond))) aa_pos_amp = sort(unique(c( - aa_ligplus_amp - , aa_plip_amp - , aa_arpeg_amp))) + aa_ligplus_amp + , aa_plip_amp + , aa_arpeg_amp))) aa_pos_amp_hbond = sort(unique(c( - aa_ligplus_amp_hbond - , aa_plip_amp_hbond))) + aa_ligplus_amp_hbond + , aa_plip_amp_hbond))) cat("\n===================================================" @@ -129,9 +129,11 @@ cat("\n===================================================" ############################################################## # var for position customisation for plots -aa_pos_lig1 = aa_pos_rna -aa_pos_lig2 = aa_pos_sam +aa_pos_drug +aa_pos_lig1 = aa_pos_sam +aa_pos_lig2 = aa_pos_rna aa_pos_lig3 = aa_pos_amp -tile_map=data.frame(tile=c("GID","DPA","CDL","Ca"), + +tile_map=data.frame(tile=c("SRY","SAM","RNA","AMP"), tile_colour=c("green","darkslategrey","navyblue","purple")) diff --git a/scripts/Header_TT.R b/scripts/Header_TT.R index 26a418d..9f5f735 100755 --- a/scripts/Header_TT.R +++ b/scripts/Header_TT.R @@ -68,8 +68,8 @@ if (!require("ggridges")) { #devtools::install_github("kassambara/ggcorrplot") if (!require ("ggbeeswarm")){ - install.packages("ggbeeswarm") - library(ggbeeswarm) + install.packages("ggbeeswarm") + library(ggbeeswarm) } if (!require("plotly")) { @@ -128,7 +128,7 @@ if(!require("stats4")) { } if(!require("data.table")) { -install.packages("data.table") + install.packages("data.table") library(data.table) } @@ -237,17 +237,17 @@ consurf_palette2 = c("0" = "yellow2" # ) consurf_colours = c( - "0" = rgb(1.00,1.00,0.59) - , "1" = rgb(0.04,0.49,0.51) - , "2" = rgb(0.29,0.69,0.75) - , "3" = rgb(0.65,0.86,0.90) - , "4" = rgb(0.84,0.94,0.94) - , "5" = rgb(1.00,1.00,1.00) - , "6" = rgb(0.98,0.92,0.96) - , "7" = rgb(0.98,0.78,0.86) - , "8" = rgb(0.94,0.49,0.67) - , "9" = rgb(0.63,0.16,0.37) - ) + "0" = rgb(1.00,1.00,0.59) + , "1" = rgb(0.04,0.49,0.51) + , "2" = rgb(0.29,0.69,0.75) + , "3" = rgb(0.65,0.86,0.90) + , "4" = rgb(0.84,0.94,0.94) + , "5" = rgb(1.00,1.00,1.00) + , "6" = rgb(0.98,0.92,0.96) + , "7" = rgb(0.98,0.78,0.86) + , "8" = rgb(0.94,0.49,0.67) + , "9" = rgb(0.63,0.16,0.37) +) ################################################## diff --git a/scripts/functions/corr_plot_data.R b/scripts/functions/corr_plot_data.R index 7bb0970..cd242cc 100644 --- a/scripts/functions/corr_plot_data.R +++ b/scripts/functions/corr_plot_data.R @@ -16,6 +16,7 @@ corr_data_extract <- function(df if ( missing(colnames_to_extract) || missing(colnames_display_key) ){ + # log10maf df$maf2 = log10(df$maf) # can't see otherwise sum(is.na(df$maf2)) @@ -40,7 +41,7 @@ corr_data_extract <- function(df , "consurf_score" , "snap2_score" , "provean_score" , "ligand_affinity_change", "mmcsm_lig" #, "ddg_dynamut", "ddg_encom", "dds_encom", "ddg_mcsm", "ddg_sdm", "ddg_duet" - ) + ) display_common_colnames = c( drug, "dst_mode" , "DUET" , "FoldX" , "DeepDDG", "Dynamut2" @@ -51,7 +52,7 @@ corr_data_extract <- function(df , "ConSurf" , "SNAP2" , "PROVEAN" , "mCSM-lig", "mmCSM-lig" # , "Dynamut" , "ENCoM-DDG" , "mCSM" , "SDM" , "DUET-d" , "ENCoM-DDS" - ) + ) if (tolower(gene)%in%geneL_normal){ colnames_to_extract = c(common_colnames) @@ -59,21 +60,21 @@ corr_data_extract <- function(df corr_df = df[,colnames_to_extract] colnames(corr_df) = display_colnames - } + } if (tolower(gene)%in%geneL_ppi2){ colnames_to_extract = c(common_colnames ,"mcsm_ppi2_affinity", ppi2Dist_colname) display_colnames = c(display_common_colnames,"mCSM-PPI2" , "PPI-Dist") corr_df = df[,colnames_to_extract] colnames(corr_df) = display_colnames - } + } if (tolower(gene)%in%geneL_na){ colnames_to_extract = c(common_colnames,"mcsm_na_affinity", naDist_colname) - display_colnames = c(display_common_colnames, "mCSM-NA", "NA-Dist") + display_colnames = c(display_common_colnames, "mCSM-NA", "NCA-Dist") corr_df = df[,colnames_to_extract] colnames(corr_df) = display_colnames - } + } # [optional] arg: extract_scaled_cols if (extract_scaled_cols){ @@ -98,19 +99,19 @@ corr_data_extract <- function(df # colnames(corr_df)[colnames(corr_df)%in%colnames_to_extract] <- display_colnames # colnames(corr_df) - cat("\nExtracted ncols:", ncol(corr_df) - ,"\nRenaming successful") - - cat("\nSneak peak...") - print(head(corr_df)) - - # Move drug column to the end - last_col = colnames(corr_df[ncol(corr_df)]) - #corr_df_f = corr_df %>% dplyr::relocate(all_of(drug), .after = last_col) - - #return(corr_df_f) - return(corr_df) - + cat("\nExtracted ncols:", ncol(corr_df) + ,"\nRenaming successful") + + cat("\nSneak peak...") + print(head(corr_df)) + + # Move drug column to the end + last_col = colnames(corr_df[ncol(corr_df)]) + #corr_df_f = corr_df %>% dplyr::relocate(all_of(drug), .after = last_col) + + #return(corr_df_f) + return(corr_df) + } } diff --git a/scripts/functions/dm_om_data.R b/scripts/functions/dm_om_data.R index 0bfacb2..a7867ce 100644 --- a/scripts/functions/dm_om_data.R +++ b/scripts/functions/dm_om_data.R @@ -222,12 +222,12 @@ dm_om_wf_lf_data <- function(df } if (tolower(gene)%in%geneL_na){ - colnames_to_extract = c(common_colnames,"mcsm_na_affinity" , "mcsm_na_scaled" , "mcsm_na_outcome" , naDist_colname) + colnames_to_extract = c(common_colnames ,"mcsm_na_affinity" , "mcsm_na_scaled" , "mcsm_na_outcome" , naDist_colname) display_colnames = c(display_common_colnames, "mcsm_na_affinity" , mcsm_na_dn , "mcsm_na_outcome" , na_dist_dn) comb_df_sl = df[, colnames_to_extract] # Rename cols: display names - colnames(comb_df) = display_colnames + colnames(comb_df_sl) = display_colnames # Affinity filtered data: mcsm-na --> naDist_colname comb_df_sl_na = comb_df_sl[comb_df_sl[[na_dist_dn]] 0 ){ -# cat( -# "\n##################################################" -# , "\nSuccessful: get_plotting_dfs.R worked!" -# , "\n###################################################\n") -# } else { -# cat( -# "\n#################################################" -# , "\nFAIL: get_plotting_dfs.R didn't complete fully!Please check" -# , "\n###################################################\n" ) +# cat( +# "\n##################################################" +# , "\nSuccessful: get_plotting_dfs.R worked!" +# , "\n###################################################\n") +# } else { +# cat( +# "\n#################################################" +# , "\nFAIL: get_plotting_dfs.R didn't complete fully!Please check" +# , "\n###################################################\n" ) # } -# -# ######################################################################## -# # clear excess variables: from the global enviornment # -# vars0 = ls(envir = .GlobalEnv)[grepl("curr_*", ls(envir = .GlobalEnv))] -# vars1 = ls(envir = .GlobalEnv)[grepl("^cols_to*", ls(envir = .GlobalEnv))] -# vars2 = ls(envir = .GlobalEnv)[grepl("pivot_cols_*", ls(envir = .GlobalEnv))] -# vars3 = ls(envir = .GlobalEnv)[grepl("expected_*", ls(envir = .GlobalEnv))] -# -# rm( infile_metadata -# , infile_params -# , vars0 -# , vars1 -# , vars2 -# , vars3) +######################################################################## +# clear excess variables: from the global enviornment + +vars0 = ls(envir = .GlobalEnv)[grepl("curr_*", ls(envir = .GlobalEnv))] +vars1 = ls(envir = .GlobalEnv)[grepl("^cols_to*", ls(envir = .GlobalEnv))] +vars2 = ls(envir = .GlobalEnv)[grepl("pivot_cols_*", ls(envir = .GlobalEnv))] +vars3 = ls(envir = .GlobalEnv)[grepl("expected_*", ls(envir = .GlobalEnv))] + +rm( infile_metadata + , infile_params + , vars0 + , vars1 + , vars2 + , vars3) diff --git a/scripts/plotting/plotting_thesis/appendix_tables.R b/scripts/plotting/plotting_thesis/appendix_tables.R index 0452e5c..1afba7a 100644 --- a/scripts/plotting/plotting_thesis/appendix_tables.R +++ b/scripts/plotting/plotting_thesis/appendix_tables.R @@ -38,18 +38,18 @@ angstroms_symbol df3 = merged_df3 cols_to_output = c("mutationinformation" - , "position" - , affinity_dist_colnames[1] - , "ligand_affinity_change" - , "ligand_outcome" - , "mmcsm_lig" - , "mmcsm_lig_outcome" - , affinity_dist_colnames[2] - , "mcsm_ppi2_affinity" - , "mcsm_ppi2_outcome" - , "maf" - , "or_mychisq" - , "pval_fisher") + , "position" + , affinity_dist_colnames[1] + , "ligand_affinity_change" + , "ligand_outcome" + , "mmcsm_lig" + , "mmcsm_lig_outcome" + , affinity_dist_colnames[2] + , "mcsm_ppi2_affinity" + , "mcsm_ppi2_outcome" + , "maf" + , "or_mychisq" + , "pval_fisher") cols_to_output df3_output = df3[, cols_to_output] @@ -65,12 +65,12 @@ colnames(df3_output) df3_output$p_adj_fdr = p.adjust(df3_output$pval_fisher, method = "fdr") df3_output$signif_fdr = df3_output$p_adj_fdr df3_output = dplyr::mutate(df3_output - , signif_fdr = case_when(signif_fdr == 0.05 ~ "." - , signif_fdr <=0.0001 ~ '****' - , signif_fdr <=0.001 ~ '***' - , signif_fdr <=0.01 ~ '**' - , signif_fdr <0.05 ~ '*' - , TRUE ~ 'ns')) + , signif_fdr = case_when(signif_fdr == 0.05 ~ "." + , signif_fdr <=0.0001 ~ '****' + , signif_fdr <=0.001 ~ '***' + , signif_fdr <=0.01 ~ '**' + , signif_fdr <0.05 ~ '*' + , TRUE ~ 'ns')) # rounding df3_output$or_mychisq = round(df3_output$or_mychisq,2) df3_output$p_adj_fdr = round(df3_output$p_adj_fdr,2) @@ -101,17 +101,17 @@ head(df3_output) df_lig = df3_output[df3_output[[LigDist_colname]]% -# dplyr::add_count(eval(parse(text = var_pos))) -# -# class(df4) -# df4 = as.data.frame(df4) -# class(df4) -# nc_change = which(colnames(df4) == "n") -# colnames(df4)[nc_change] <- "pos_count" -# class(df4) -# -# setDT(df4)[, pos_count2 := .N, by = .(eval(parse(text = "position")))] -# class(df4) -# -# all(df4$pos_count==df4$pos_count2) -# -# # %>% -# #group_by(pos_count = position) -# -# # df4 = -# # df4 %>% -# # dplyr::group_by(position) %>% -# # count(position) -#foo2 = df4[, c("mutationinformation", "position", "pos_count")] - -##################################################################### -# ------------------------------ -# bp site site count: ALL -# <10 Ang ligand -# ------------------------------ -posC_all = site_snp_count_bp(plotdf = df3 - , df_colname = "position" - , xaxis_title = "Number of nsSNPs" - , yaxis_title = "Number of Sites" - , subtitle_size = 20) # ------------------------------ # bp site site count: mCSM-lig @@ -532,55 +193,233 @@ common_bp_title = paste0("Sites <", DistCutOff, angstroms_symbol) posC_lig = site_snp_count_bp(plotdf = df3_lig , df_colname = "position" , xaxis_title = "Number of nsSNPs" - , yaxis_title = "Number of Sites"#+ annotate("text", x = 1.5, y = 2.2, label = "Text No. 1") - #, subtitle_text = paste0(common_bp_title, " ligand") + , yaxis_title = "Number of Sites" + , subtitle_colour = "chocolate4" + , subtitle_text = "" , subtitle_size = 8 - , subtitle_colour = subtitle_colour) + , geom_ls = 2.6 + , leg_text_size = 10 + , axis_text_size = 10 + , axis_label_size = 10) + posC_lig # ------------------------------ # bp site site count: ppi2 # < 10 Ang interface # ------------------------------ +if (tolower(gene)%in%geneL_ppi2){ + + posC_ppi2 = site_snp_count_bp(plotdf = df3_ppi2 + , df_colname = "position" + , xaxis_title = "Number of nsSNPs" + , yaxis_title = "Number of Sites" + , subtitle_colour = "chocolate4" + , subtitle_text = "" + , subtitle_size = 8 + , geom_ls = 2.6 + , leg_text_size = 10 + , axis_text_size = 10 + , axis_label_size = 10) + posC_ppi2 +} -posC_ppi2 = site_snp_count_bp(plotdf = df3_ppi2 - , df_colname = "position" - , xaxis_title = "Number of nsSNPs" - , yaxis_title = "Number of Sites" - , subtitle_text = paste0(common_bp_title, " interface") - , subtitle_size = 20 - , subtitle_colour = subtitle_colour) -posC_ppi2 # ------------------------------ -#FIXME: bp site site count: na -# < 10 Ang TBC +# bp site site count: NCA dist +# < 10 Ang nca # ------------------------------ -# posC_na = site_snp_count_bp(plotdf = df3_na -# , df_colname = "position" -# , xaxis_title = "" -# , yaxis_title = "") - - -#=========================== -# output: SITE SNP count: -# all + affinity -#========================== -# my_label_size = 25 -# pos_count_combined_CLP = paste0(outdir_images -# ,tolower(gene) -# ,"_pos_count_PS_AFF.svg") -# -# -# svg(pos_count_combined_CLP, width = 20, height = 5.5) -# print(paste0("plot filename:", pos_count_combined_CLP)) -# -# cowplot::plot_grid(posC_all, posC_lig, posC_ppi2 -# #, posC_na -# , nrow = 1 -# , ncol = 3 -# , labels = "AUTO" -# , label_size = my_label_size) -# -# dev.off() +if (tolower(gene)%in%geneL_na){ + + posC_nca = site_snp_count_bp(plotdf = df3_na + , df_colname = "position" + , xaxis_title = "Number of nsSNPs" + , yaxis_title = "Number of Sites" + , subtitle_colour = "chocolate4" + , subtitle_text = "" + , subtitle_size = 8 + , geom_ls = 2.6 + , leg_text_size = 10 + , axis_text_size = 10 + , axis_label_size = 10) + posC_nca +} #=============================================================== + + +# ------------------------------ +# bp site site count: ALL +# <10 Ang ligand +# ------------------------------ +posC_all = site_snp_count_bp(plotdf = df3 + , df_colname = "position" + , xaxis_title = "Number of nsSNPs" + , yaxis_title = "Number of Sites" + , subtitle_colour = "chocolate4" + , subtitle_text = "All mutations sites" + , subtitle_size = 8 + , geom_ls = 2.6 + , leg_text_size = 10 + , axis_text_size = 10 + , axis_label_size = 10) +posC_all +################################################################## +consurfP = stability_count_bp(plotdf = df3 + , df_colname = "consurf_outcome" + #, leg_title = "ConSurf" + #, label_categories = labels_consurf + , yaxis_title = "Number of nsSNPs" + , leg_position = "top" + , subtitle_text = "ConSurf" + , bar_fill_values = consurf_colours # from globals + , subtitle_colour= "black" + , sts = 10 + , lts = 8 + , ats = 8 + , als = 8 + , ltis = 11 + , geom_ls = 2) + +consurfP + +#################### +# Sensitivity count: Mutations +#################### +table(df3$sensitivity) + +rect_sens=data.frame(mutation_class=c("Resistant","Sensitive") + , tile_colour =c("red","blue") + , numbers = c(table(df3$sensitivity)[[1]], table(df3$sensitivity)[[2]])) + + + +sensP = ggplot(rect_sens, aes(mutation_class, y = 0 + , fill = tile_colour + , label = paste0("n=", numbers) + )) + + geom_tile(width = 1, height = 1) + # make square tiles + geom_label(color = "black", size = 1.7,fill = "white", alpha=0.7) + # add white text in the middle + scale_fill_identity(guide = "none") + # color the tiles with the colors in the data frame + coord_fixed() + # make sure tiles are square + #coord_flip()+ scale_x_reverse() + + # theme_void() # remove any axis markings + theme_nothing() # remove any axis markings +sensP + +# sensP2 = sensP + +# coord_flip() + scale_x_reverse() +# sensP2 + +############################################################## +#=================== +# Stability +#=================== +# duetP +duetP = stability_count_bp(plotdf = df3 + , df_colname = "duet_outcome" + , leg_title = "mCSM-DUET" + #, label_categories = labels_duet + , yaxis_title = "Number of nsSNPs" + , leg_position = "none" + , subtitle_text = "mCSM-DUET" + , bar_fill_values = c("#F8766D", "#00BFC4") + , subtitle_colour= "black" + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5 +) +duetP + +# foldx +foldxP = stability_count_bp(plotdf = df3 + , df_colname = "foldx_outcome" + #, leg_title = "FoldX" + #, label_categories = labels_foldx + , yaxis_title = "" + , leg_position = "none" + , subtitle_text = "FoldX" + , bar_fill_values = c("#F8766D", "#00BFC4") + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5 +) +foldxP + +# deepddg +deepddgP = stability_count_bp(plotdf = df3 + , df_colname = "deepddg_outcome" + #, leg_title = "DeepDDG" + #, label_categories = labels_deepddg + , yaxis_title = "" + , leg_position = "none" + , subtitle_text = "DeepDDG" + , bar_fill_values = c("#F8766D", "#00BFC4") + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5 +) +deepddgP + +# deepddg +dynamut2P = stability_count_bp(plotdf = df3 + , df_colname = "ddg_dynamut2_outcome" + #, leg_title = "Dynamut2" + #, label_categories = labels_ddg_dynamut2_outcome + , yaxis_title = "" + , leg_position = "none" + , subtitle_text = "Dynamut2" + , bar_fill_values = c("#F8766D", "#00BFC4") + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5 +) +dynamut2P + +# provean +proveanP = stability_count_bp(plotdf = df3 + , df_colname = "provean_outcome" + #, leg_title = "PROVEAN" + #, label_categories = labels_provean + , yaxis_title = "Number of nsSNPs" + , leg_position = "none" # top + , subtitle_text = "PROVEAN" + , bar_fill_values = c("#D01C8B", "#F1B6DA") # light pink and deep + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5 +) +proveanP + +# snap2 +snap2P = stability_count_bp(plotdf = df3 + , df_colname = "snap2_outcome" + #, leg_title = "SNAP2" + #, label_categories = labels_snap2 + , yaxis_title = "" + , leg_position = "none" # top + , subtitle_text = "SNAP2" + , bar_fill_values = c("#D01C8B", "#F1B6DA") # light pink and deep + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5) +snap2P + +##################################################################################### diff --git a/scripts/plotting/plotting_thesis/basic_barplots_layout.R b/scripts/plotting/plotting_thesis/basic_barplots_layout.R deleted file mode 100644 index ecbcd8e..0000000 --- a/scripts/plotting/plotting_thesis/basic_barplots_layout.R +++ /dev/null @@ -1,198 +0,0 @@ -# source basic_barplots.R - -#============ -# Plot labels -#============ -tit1 = "Stability outcome" -tit2 = "Affinity outcome" -tit3 = "Conservation outcome" -pt_size = 30 - - -theme_georgia <- function(...) { - theme_gray(base_family = "sans", ...) + - theme(plot.title = element_text(face = "bold")) -} - - -title_theme <- calc_element("plot.title", theme_georgia()) - -pt1 = ggdraw() + - draw_label( - tit1, - fontfamily = title_theme$family, - fontface = title_theme$face, - #size = title_theme$size - size = pt_size - ) - -pt2 = ggdraw() + - draw_label( - tit2, - fontfamily = title_theme$family, - fontface = title_theme$face, - size = pt_size - ) - -pt3 = ggdraw() + - draw_label( - tit3, - fontfamily = title_theme$family, - fontface = title_theme$face, - size = pt_size - ) - -# extract common legend -common_legend_outcome = get_legend(mLigP + - guides(color = guide_legend(nrow = 1)) + - theme(legend.position = "top")) - - - -my_label_size = 25 -#====================== -# Output plot function -#====================== -OutPlotBP = function(x){ -cowplot::plot_grid( - cowplot::plot_grid(pt1, - common_legend_outcome, - cowplot::plot_grid( duetP, foldxP - , deepddgP, dynamut2P - , nrow = 2 - , ncol = 2 - , labels = c("A", "B", "C","D") - , label_size = my_label_size - ) - , ncol = 1 - , rel_heights = c(7, 3, 90)), - - cowplot::plot_grid(pt2, - cowplot::plot_grid(mLigP, mmLigP, ppi2P - , nrow = 1 - , ncol = 3 - , labels = c("E","F", "G") - , label_size = my_label_size - ) - , ncol = 1 - , rel_heights = c(1, 9)), - - cowplot::plot_grid(pt3, - cowplot::plot_grid(consurfP, proveanP, snap2P - , nrow = 1 - , ncol = 3 - , labels = c("H", "I", "J") - , labels_x = 0.2 - , label_size = my_label_size - , rel_widths = c(0.2, 0.2, 0.2) - ) - , ncol = 1 - , rel_heights = c(0.07, 0.93) - ), - - nrow = 3, - rel_heights = c(0.58, 0.25, 0.27), - align = "hv" -) -} - -#===================== -# OutPlot: svg and png -#====================== -#ratio 11.69 by 8.27 -w = 8.27*2 -h = 11.69*2 - -#svg -bp_all_CLP = paste0(outdir_images - ,tolower(gene) - ,"_bp_all_CL.svg") -cat(paste0("plot filename:", bp_all_CLP)) - -svg(bp_all_CLP, width = w, height = h) -OutPlotBP() -dev.off() - -#png -bp_all_CLP_png = paste0(outdir_images - ,tolower(gene) - ,"_bp_all_CL.png") -cat(paste0("plot filename:", bp_all_CLP_png)) - -png(bp_all_CLP_png, width = w, height = h, units = "in", res = 300 ) -OutPlotBP() -dev.off() - -##################################################################### -##################################################################### -# ------------------------------ -# bp site site count: ALL -# <10 Ang ligand -# ------------------------------ - -posC_all = site_snp_count_bp(plotdf = df3 - , df_colname = "position" - , xaxis_title = "Number of nsSNPs" - , yaxis_title = "Number of Sites" - , subtitle_size = 20) - -# ------------------------------ -# bp site site count: mCSM-lig -# < 10 Ang ligand -# ------------------------------ -common_bp_title = paste0("Sites <", DistCutOff, angstroms_symbol) - -posC_lig = site_snp_count_bp(plotdf = df3_lig - , df_colname = "position" - , xaxis_title = "Number of nsSNPs" - , yaxis_title = "Number of Sites"#+ annotate("text", x = 1.5, y = 2.2, label = "Text No. 1") - , subtitle_text = paste0(common_bp_title, " ligand") - , subtitle_size = 20 - , subtitle_colour = subtitle_colour) -# ------------------------------ -# bp site site count: ppi2 -# < 10 Ang interface -# ------------------------------ - -posC_ppi2 = site_snp_count_bp(plotdf = df3_ppi2 - , df_colname = "position" - , xaxis_title = "Number of nsSNPs" - , yaxis_title = "Number of Sites" - , subtitle_text = paste0(common_bp_title, " interface") - , subtitle_size = 20 - , subtitle_colour = subtitle_colour) - -# ------------------------------ -#FIXME: bp site site count: na -# < 10 Ang TBC -# ------------------------------ -# posC_na = site_snp_count_bp(plotdf = df3_na -# , df_colname = "position" -# , xaxis_title = "" -# , yaxis_title = "") - - -#=========================== -# output: SITE SNP count: -# all + affinity -#========================== -my_label_size = 25 -pos_count_combined_CLP = paste0(outdir_images - ,tolower(gene) - ,"_pos_count_PS_AFF.svg") - - -svg(pos_count_combined_CLP, width = 20, height = 5.5) -print(paste0("plot filename:", pos_count_combined_CLP)) - -cowplot::plot_grid(posC_all, posC_lig, posC_ppi2 - #, posC_na - , nrow = 1 - , ncol = 3 - , labels = "AUTO" - , label_size = my_label_size) - -dev.off() - - -#=============================================================== diff --git a/scripts/plotting/plotting_thesis/basic_barplots_layout_v2.R b/scripts/plotting/plotting_thesis/basic_barplots_layout_v2.R index 77062b8..5ae7e02 100644 --- a/scripts/plotting/plotting_thesis/basic_barplots_layout_v2.R +++ b/scripts/plotting/plotting_thesis/basic_barplots_layout_v2.R @@ -10,8 +10,9 @@ mmLigP posC_lig ppi2P posC_ppi2 -peP sensP +peP + #======================== # Common title settings #========================= @@ -157,9 +158,9 @@ ppi2_affT = ggdraw() + ########################################################### p2 = cowplot::plot_grid(cowplot::plot_grid(ppi2_affT, common_legend_outcome, nrow=2), cowplot::plot_grid(ppi2P, posC_ppi2 - , nrow = 1 - , rel_widths = c(1.2,1.8) - , align = "h"), + , nrow = 1 + , rel_widths = c(1.2,1.8) + , align = "h"), nrow = 2, rel_heights = c(1,8) ) @@ -184,7 +185,7 @@ p3 = cowplot::plot_grid(cowplot::plot_grid(peT_allT, nrow = 2 align = "v", axis = "lr", rel_heights = c(1,8) - ), + ), rel_heights = c(1,18), nrow = 2,axis = "lr") p3 @@ -207,7 +208,7 @@ cowplot::plot_grid(p1, p2, p3 , label_size = 12 , rel_widths = c(3,2,2) #, rel_heights = c(1) - ) +) dev.off() ################################################## @@ -234,8 +235,8 @@ h = 3 # dev.off() conCLP = paste0(outdir_images - ,tolower(gene) - ,"_consurf_BP.png") + ,tolower(gene) + ,"_consurf_BP.png") print(paste0("plot filename:", conCLP)) png(conCLP, units = "in", width = w, height = h, res = 300 ) @@ -243,15 +244,27 @@ consurfP dev.off() #================================ -# Sensitivity numbers: geom_tile +# Sensitivity mutation numbers: geom_tile #================================ sensCLP = paste0(outdir_images - ,tolower(gene) - ,"_sensN_tile.png") + ,tolower(gene) + ,"_sensN_tile.png") print(paste0("plot filename:", sensCLP)) png(sensCLP, units = "in", width = 1, height = 1, res = 300 ) sensP dev.off() +#================================ +# Sensitivity SITE numbers: geom_tile +#================================ +sens_siteCLP = paste0(outdir_images + ,tolower(gene) + ,"_sens_siteC_tile.png") + +print(paste0("plot filename:", sens_siteCLP)) +png(sens_siteCLP, units = "in", width = 1, height = 1, res = 300 ) +sens_siteP +dev.off() + ########################################################### diff --git a/scripts/plotting/plotting_thesis/basic_barplots2.R b/scripts/plotting/plotting_thesis/basic_barplots_v2.R similarity index 79% rename from scripts/plotting/plotting_thesis/basic_barplots2.R rename to scripts/plotting/plotting_thesis/basic_barplots_v2.R index b75d049..b6fcfea 100644 --- a/scripts/plotting/plotting_thesis/basic_barplots2.R +++ b/scripts/plotting/plotting_thesis/basic_barplots_v2.R @@ -25,16 +25,21 @@ #============= # Data: Input #============== -#source("~/git/LSHTM_analysis/config/alr.R") -source("~/git/LSHTM_analysis/config/embb.R") -#source("~/git/LSHTM_analysis/config/katg.R") -#source("~/git/LSHTM_analysis/config/gid.R") #source("~/git/LSHTM_analysis/config/pnca.R") +#source("~/git/LSHTM_analysis/config/embb.R") +#source("~/git/LSHTM_analysis/config/gid.R") + +source("~/git/LSHTM_analysis/config/alr.R") +#source("~/git/LSHTM_analysis/config/katg.R") #source("~/git/LSHTM_analysis/config/rpob.R") source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") -source("~/git/LSHTM_analysis/scripts/plotting/plotting_colnames.R") +#source("~/git/LSHTM_analysis/scripts/plotting/plotting_colnames.R") sourced by above +# sanity check +cat("\nSourced plotting cols as well:", length(plotting_cols)) + +#################################################### class(merged_df3) merged_df3 = as.data.frame(merged_df3) @@ -52,7 +57,7 @@ merged_df3 = merged_df3[, !colnames(merged_df3)%in%c("pos_count")] head(merged_df3$pos_count) df3 = merged_df3[, colnames(merged_df3)%in%plotting_cols] - +"nca_distance"%in%colnames(df3) #======= # output @@ -129,25 +134,54 @@ mmLigP # barplot for ppi2 affinity # <10 Ang of interface #------------------------------ -ppi2P = stability_count_bp(plotdf = df3_ppi2 - , df_colname = "mcsm_ppi2_outcome" - #, leg_title = "mCSM-ppi2" - #, label_categories = labels_ppi2 - #, bp_plot_title = paste(common_bp_title, "PP-interface") - - , yaxis_title = "Number of nsSNPs" - , leg_position = "none" - , subtitle_text = "mCSM-ppi2" - , bar_fill_values = c("#F8766D", "#00BFC4") - , subtitle_colour= "black" - , sts = 10 - , lts = 8 - , ats = 12 - , als = 11 - , ltis = 11 - , geom_ls = 2.5 - ) -ppi2P +if (tolower(gene)%in%geneL_ppi2){ + ppi2P = stability_count_bp(plotdf = df3_ppi2 + , df_colname = "mcsm_ppi2_outcome" + #, leg_title = "mCSM-ppi2" + #, label_categories = labels_ppi2 + #, bp_plot_title = paste(common_bp_title, "PP-interface") + + , yaxis_title = "Number of nsSNPs" + , leg_position = "none" + , subtitle_text = "mCSM-ppi2" + , bar_fill_values = c("#F8766D", "#00BFC4") + , subtitle_colour= "black" + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5 + ) + ppi2P +} +#---------------------------- +# barplot for ppi2 affinity +# <10 Ang of interface +#------------------------------ +if (tolower(gene)%in%geneL_na){ + + nca_distP = stability_count_bp(plotdf = df3_na + , df_colname = "mcsm_na_outcome" + #, leg_title = "mCSM-NA" + #, label_categories = + #, bp_plot_title = paste(common_bp_title, "Dist to NA") + + , yaxis_title = "Number of nsSNPs" + , leg_position = "none" + , subtitle_text = "mCSM-NA" + , bar_fill_values = c("#F8766D", "#00BFC4") + , subtitle_colour= "black" + , sts = 10 + , lts = 8 + , ats = 12 + , als = 11 + , ltis = 11 + , geom_ls = 2.5 + ) + nca_distP +} + ##################################################################### # ------------------------------ @@ -173,19 +207,43 @@ posC_lig # bp site site count: ppi2 # < 10 Ang interface # ------------------------------ +if (tolower(gene)%in%geneL_ppi2){ + + posC_ppi2 = site_snp_count_bp(plotdf = df3_ppi2 + , df_colname = "position" + , xaxis_title = "Number of nsSNPs" + , yaxis_title = "Number of Sites" + , subtitle_colour = "chocolate4" + , subtitle_text = "" + , subtitle_size = 8 + , geom_ls = 2.6 + , leg_text_size = 10 + , axis_text_size = 10 + , axis_label_size = 10) + posC_ppi2 +} + +# ------------------------------ +# bp site site count: NCA dist +# < 10 Ang nca +# ------------------------------ +if (tolower(gene)%in%geneL_na){ + + posC_nca = site_snp_count_bp(plotdf = df3_na + , df_colname = "position" + , xaxis_title = "Number of nsSNPs" + , yaxis_title = "Number of Sites" + , subtitle_colour = "chocolate4" + , subtitle_text = "" + , subtitle_size = 8 + , geom_ls = 2.6 + , leg_text_size = 10 + , axis_text_size = 10 + , axis_label_size = 10) + posC_nca +} + -posC_ppi2 = site_snp_count_bp(plotdf = df3_ppi2 - , df_colname = "position" - , xaxis_title = "Number of nsSNPs" - , yaxis_title = "Number of Sites" - , subtitle_colour = "chocolate4" - , subtitle_text = "" - , subtitle_size = 8 - , geom_ls = 2.6 - , leg_text_size = 10 - , axis_text_size = 10 - , axis_label_size = 10) -posC_ppi2 #=============================================================== # PE count rects <- data.frame(x = 1:6, @@ -246,7 +304,7 @@ posC_all = site_snp_count_bp(plotdf = df3 , leg_text_size = 10 , axis_text_size = 10 , axis_label_size = 10) - +posC_all ################################################################## #------------------------------ @@ -290,10 +348,8 @@ consurfP = stability_count_bp(plotdf = df3 consurfP - - #################### -# Sensitivity count +# Sensitivity count: Mutations #################### table(df3$sensitivity) @@ -320,6 +376,36 @@ sensP # sensP2 = sensP + # coord_flip() + scale_x_reverse() # sensP2 +#=============================== +# Sensitivity count: Site +#============================== +table(df3$sensitivity) +#-------- +# embb +#-------- +#rsc = 54 +#ccc = 46 +#ssc = 470 + + +rect_rs_siteC =data.frame(mutation_class=c("A_Resistant sites" + , "B_Common sites" + , "C_Sensitive sites"), + tile_colour =c("red", + "purple", + "blue"), + numbers = c(rsc, ccc, ssc), + order = c(1, 2, 3)) + +rect_rs_siteC$labels = paste0(rect_rs_siteC$mutation_class, "\nn=", rect_rs_siteC$ numbers) + +sens_siteP = ggplot(rect_rs_siteC, aes(mutation_class, y = 0, + fill = tile_colour, + label = paste0("n=", numbers))) + + geom_tile(width = 1, height = 1) + + geom_label(color = "black", size = 1.7,fill = "white", alpha=0.7) + + theme_nothing() +sens_siteP ############################################################## #=================== @@ -360,7 +446,7 @@ foldxP = stability_count_bp(plotdf = df3 , ltis = 11 , geom_ls = 2.5 ) - +foldxP # deepddg deepddgP = stability_count_bp(plotdf = df3 @@ -378,7 +464,7 @@ deepddgP = stability_count_bp(plotdf = df3 , ltis = 11 , geom_ls = 2.5 ) - +deepddgP # deepddg dynamut2P = stability_count_bp(plotdf = df3 @@ -398,7 +484,6 @@ dynamut2P = stability_count_bp(plotdf = df3 ) dynamut2P - # provean proveanP = stability_count_bp(plotdf = df3 , df_colname = "provean_outcome" @@ -415,6 +500,7 @@ proveanP = stability_count_bp(plotdf = df3 , ltis = 11 , geom_ls = 2.5 ) +proveanP # snap2 snap2P = stability_count_bp(plotdf = df3 @@ -431,7 +517,7 @@ snap2P = stability_count_bp(plotdf = df3 , als = 11 , ltis = 11 , geom_ls = 2.5) - +snap2P ############################################################## diff --git a/scripts/plotting/plotting_thesis/corr_plots_thesis.R b/scripts/plotting/plotting_thesis/corr_plots_thesis.R index b9c6bf9..fe47afa 100644 --- a/scripts/plotting/plotting_thesis/corr_plots_thesis.R +++ b/scripts/plotting/plotting_thesis/corr_plots_thesis.R @@ -263,7 +263,7 @@ if (tolower(gene)%in%geneL_ppi2){ # NA affinity #================ if (tolower(gene)%in%geneL_na){ - corr_df_na = corr_df_na[corr_df_na["NA-Dist"]DistCutOff,"mCSM-NA"]=0 +# corr_affinity_df[corr_affinity_df["NCA-Dist"]>DistCutOff,"mCSM-NA"]=0 # } # count 0 @@ -89,10 +91,12 @@ corr_ps_colnames = c(static_cols , "Dynamut2" , aff_dist_cols , "dst_mode") + corr_df_ps = corr_plotdf[, corr_ps_colnames] # Plot #1 plot_corr_df_ps = my_gg_pairs(corr_df_ps, plot_title="Stability estimates") + ########################################################## #================ # Conservation @@ -101,7 +105,7 @@ corr_conservation_cols = c( static_cols , "ConSurf" , "SNAP2" , "PROVEAN" - , aff_dist_cols + #, aff_dist_cols , "dst_mode" ) @@ -166,6 +170,6 @@ png(paste0(outdir_images ,"_CorrC.png"), height =7, width=7, unit="in",res=300) cowplot::plot_grid(ggmatrix_gtable(plot_corr_df_aff), - labels = "C", + labels = "C", label_size = 12) dev.off() diff --git a/scripts/plotting/plotting_thesis/linage_bp_dist.R b/scripts/plotting/plotting_thesis/linage_bp_dist.R index eb7c2ec..bfe6118 100644 --- a/scripts/plotting/plotting_thesis/linage_bp_dist.R +++ b/scripts/plotting/plotting_thesis/linage_bp_dist.R @@ -61,7 +61,7 @@ lin_countP = lin_count_bp(lf_data = lineage_dfL[['lin_lf']] , y_scale_percent = FALSE , y_label = c("Count") ) - +lin_countP #=============================== # lineage SNP diversity count #=============================== @@ -88,7 +88,7 @@ lin_diversityP = lin_count_bp_diversity(lf_data = lineage_dfL[['lin_wf']] , subtitle_text = NULL , sts = 20 , subtitle_colour = "#350E20FF") - +lin_diversityP #============================================= # Output plots: Lineage count and Diversity #=============================================