diff --git a/scripts/functions/logoP_msa.R b/scripts/functions/logoP_msa.R index e630295..a120d20 100644 --- a/scripts/functions/logoP_msa.R +++ b/scripts/functions/logoP_msa.R @@ -47,13 +47,14 @@ LogoPlotMSA <- function(unified_msa , my_logo_col = "chemistry" , plot_positions , y_breaks - , x_lab_mut = "nsSNP-position" + , x_lab_mut = "nsSNP position" , y_lab_mut , x_ats = 13 # text size , x_tangle = 90 # text angle - , x_axis_offset = 0.05 # dist b/w y-axis and plot start - , x_axis_offset_filtered = 0.2 - , y_axis_offset = 0.05 + , x_axis_offset = 0 # dist b/w y-axis and plot start + , x_axis_offset_filtered = 0 + , y_axis_offset = 0 + , y_axis_increment = 1 , y_ats = 13 , y_tangle = 0 , x_tts = 13 # title size @@ -102,7 +103,7 @@ LogoPlotMSA <- function(unified_msa y_rlow = seq(0, ylim_low, length.out = 3); y_rlow ylim_up <- ceiling(max(data_logo_mut)) + 5; ylim_up - cat("\nY-axis upper limit:", ylim_up) + cat("\nY-axis upper limit:", ylim_up) y_rup = seq(0, ylim_up, by = 2); y_rup }else{ ylim_low = ylim_low + (-0.5) @@ -192,7 +193,7 @@ LogoPlotMSA <- function(unified_msa #--------- # MSA mut #--------- - cat("\n===========================================" + cat("\n===========================================" , "\nGenerated PFM mut: No filtering" , "\n===========================================") @@ -301,7 +302,7 @@ LogoPlotMSA <- function(unified_msa , '\nSelected colour scheme:', my_logo_col , "\nUsing grey theme") - theme_bgc = "grey" + theme_bgc = "white" xfont_bgc = "black" yfont_bgc = "black" xtt_col = "black" @@ -317,10 +318,10 @@ LogoPlotMSA <- function(unified_msa #------------------- # Mutant logo plot #------------------- - p0 = ggseqlogo(plot_mut_edM - , method = msa_method - , col_scheme = my_logo_col - , seq_type = 'auto') + + p0 = ggplot() + geom_logo(plot_mut_edM + , method = msa_method + , col_scheme = my_logo_col + , seq_type = 'auto') + theme(legend.position = leg_pos , legend.direction = leg_dir @@ -335,6 +336,7 @@ LogoPlotMSA <- function(unified_msa , vjust = 0.4 , colour = xfont_bgc) #, axis.text.y = element_blank() + , axis.ticks=element_blank() , axis.text.y = element_text(size = y_ats , angle = y_tangle , hjust = 1 @@ -344,12 +346,24 @@ LogoPlotMSA <- function(unified_msa , colour = xtt_col) , axis.title.y = element_text(size = y_tts , colour = ytt_col) - , plot.background = element_rect(fill = theme_bgc)) + - + , panel.grid=element_blank() + , plot.background = element_rect(fill = theme_bgc, colour=NA) + , panel.background = element_rect(fill = "transparent", colour=NA) + + ) + + labs(y=y_label) + xlab(x_lab_mut) if (missing(plot_positions)){ ed_mut_logo_P = p0 + + scale_y_continuous( + expand = c(0,0), + breaks = seq( + 0, + (y_lim), + by = y_axis_increment + ) + ) + scale_x_discrete(breaks = msa_pos , expand = c(x_axis_offset, 0) , labels = msa_pos @@ -357,55 +371,37 @@ LogoPlotMSA <- function(unified_msa }else{ ed_mut_logo_P = p0 + + scale_y_continuous( + expand = c(0,0)#, + # breaks = seq( + # 0, + # (y_lim), + # by = y_axis_increment + #) + ) + + # scale_x_continuous(expand = c(0,0)) #+ + scale_x_discrete(breaks = i_extract , expand = c(x_axis_offset_filtered, 0) , labels = i_extract , limits = factor(i_extract)) - - } - - if (logo_type == "EDLogo"){ - ed_mut_logo_P = ed_mut_logo_P + - scale_y_continuous(limits = c(ylim_low, ylim_up) - , breaks = ylim_scale - , expand = c(0, y_axis_offset)) + - geom_hline(yintercept = 0 - , linetype = "solid" - , color = "grey" - , size = 1) - } - - if (missing(y_lab_mut)){ - ed_mut_logo_P = ed_mut_logo_P + ylab(y_label) - } else{ - ed_mut_logo_P = ed_mut_logo_P + ylab(y_lab_mut) } cat('\nDone: MSA plot for mutations') - #return(msa_mut_logoP) - PlotlogolasL[['ed_mut_logoP']] <- ed_mut_logo_P - - #--------------------------------- - # Wild-type MSA: gene_fasta file - #--------------------------------- - p1 = ggseqlogo(plot_wt_edM - #, facet = "grid" - , method = msa_method - , col_scheme = my_logo_col - , seq_type = 'aa') + + #### Wild-type MSA: gene_fasta file #### + p1 = ggplot() + geom_logo(plot_wt_edM + #, facet = "grid" + , method = msa_method + , col_scheme = my_logo_col + , seq_type = 'aa') + theme(legend.position = "none" , 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.x = element_blank() + , axis.ticks=element_blank() , axis.text.y = element_blank() , axis.title.x = element_text(size = x_tts @@ -413,51 +409,34 @@ LogoPlotMSA <- function(unified_msa , axis.title.y = element_text(size = y_tts , colour = ytt_col) - , plot.background = element_rect(fill = theme_bgc)) + - + , panel.grid=element_blank() + , plot.background = element_rect(fill = theme_bgc, colour=NA) + , panel.background = element_rect(fill = "transparent", colour=NA) + , plot.margin = margin(r=0,l=0, unit="pt") + + ) + + scale_y_discrete(expand = c(0,0)) + ylab("") + xlab("Wild-type position") if (missing(plot_positions)){ # No y-axis needed - ed_wt_logo_P = p1 + - scale_x_discrete(breaks = wt_pos - , expand = c(x_axis_offset, 0) - , labels = wt_pos - , limits = factor(wt_pos)) + ed_wt_logo_P = p1# + }else{ ed_wt_logo_P = p1 + - scale_x_discrete(breaks = i_extract - , expand = c(x_axis_offset_filtered, 0) - , labels = i_extract - , limits = factor(i_extract)) + scale_x_discrete(expand = c(0, 0), + breaks = i_extract, + #labels = i_extract, + limits = factor(i_extract) + ) } - - cat('\nDone: MSA plot for WT') - #return(msa_wt_logoP) - PlotlogolasL[['ed_wt_logoP']] <- ed_wt_logo_P - - #========================================= - # Output - # Combined plot: logo ED/other logo plot - # customised for ggseqlogo - #========================================= - - cat('\nDone: ed_mut_logoP + ed_wt_logoP') - - # 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) - - LogoED_comb = cowplot::plot_grid(PlotlogolasL[['ed_mut_logoP']] - , PlotlogolasL[['ed_wt_logoP']] - , nrow = 2 - , align = "v" - , axis='lr' - , rel_heights = c(3/4, 1/4)) - - return(LogoED_comb) + cowplot::plot_grid(ed_mut_logo_P + , ed_wt_logo_P + , nrow = 2 + , align = "v" + #, axis='lr' + , rel_heights = c(3/4, 1/4)) } #LogoPlotMSA(unified_msa) diff --git a/scripts/functions/position_count_bp.R b/scripts/functions/position_count_bp.R index d334ae4..a53fe21 100755 --- a/scripts/functions/position_count_bp.R +++ b/scripts/functions/position_count_bp.R @@ -48,7 +48,7 @@ site_snp_count_bp <- function (plotdf #------------------------------------------- # adding column: snpcount for each position #------------------------------------------- - #setDT(plotdf)[, pos_count_check := .N, by = .(eval(parse(text = df_colname)))] + #setDT(plotdf)[, position_count_check := .N, by = .(eval(parse(text = df_colname)))] # from dplyr plotdf = plotdf %>% @@ -57,21 +57,21 @@ site_snp_count_bp <- function (plotdf plotdf = as.data.frame(plotdf) class(plotdf) nc_change = which(colnames(plotdf) == "n") - colnames(plotdf)[nc_change] <- "pos_count" + colnames(plotdf)[nc_change] <- "position_count" class(plotdf) - # if (all(plotdf$pos_count==plotdf$pos_count_check) ){ - # cat("\nPASS: pos_count column created") - # plotdf = plotdf[, !colnames(plotdf)%in%c("pos_count_check")] + # if (all(plotdf$position_count==plotdf$position_count_check) ){ + # cat("\nPASS: position_count column created") + # plotdf = plotdf[, !colnames(plotdf)%in%c("position_count_check")] # }else{ # stop("\nAbort: pos count numbes mismatch from dplyr and data.table") # } cat("\nCumulative nssnp count\n" - , table(plotdf$pos_count)) + , table(plotdf$position_count)) # calculating total no. of mutations - tot_muts = sum(table(plotdf$pos_count)) + tot_muts = sum(table(plotdf$position_count)) # sanity check @@ -92,10 +92,10 @@ site_snp_count_bp <- function (plotdf # creating df: average count of snpcount for each position # created in earlier step #------------------------------------------------------- - # use group by on pos_count + # use group by on position_count snpsBYpos_df <- plotdf %>% dplyr::group_by(eval(parse(text = df_colname))) %>% - dplyr::summarise(snpsBYpos = mean(pos_count)) # changed from summarize! + dplyr::summarise(snpsBYpos = mean(position_count)) # changed from summarize! cat("\nnssnp count per position\n" , table(snpsBYpos_df$snpsBYpos)