This commit is contained in:
Tanushree Tunstall 2022-08-25 17:57:58 +01:00
parent cd76a4b919
commit afa9166ca8
2 changed files with 73 additions and 94 deletions

View file

@ -47,13 +47,14 @@ LogoPlotMSA <- function(unified_msa
, my_logo_col = "chemistry" , my_logo_col = "chemistry"
, plot_positions , plot_positions
, y_breaks , y_breaks
, x_lab_mut = "nsSNP-position" , x_lab_mut = "nsSNP position"
, y_lab_mut , y_lab_mut
, x_ats = 13 # text size , x_ats = 13 # text size
, x_tangle = 90 # text angle , x_tangle = 90 # text angle
, x_axis_offset = 0.05 # dist b/w y-axis and plot start , x_axis_offset = 0 # dist b/w y-axis and plot start
, x_axis_offset_filtered = 0.2 , x_axis_offset_filtered = 0
, y_axis_offset = 0.05 , y_axis_offset = 0
, y_axis_increment = 1
, y_ats = 13 , y_ats = 13
, y_tangle = 0 , y_tangle = 0
, x_tts = 13 # title size , x_tts = 13 # title size
@ -301,7 +302,7 @@ LogoPlotMSA <- function(unified_msa
, '\nSelected colour scheme:', my_logo_col , '\nSelected colour scheme:', my_logo_col
, "\nUsing grey theme") , "\nUsing grey theme")
theme_bgc = "grey" theme_bgc = "white"
xfont_bgc = "black" xfont_bgc = "black"
yfont_bgc = "black" yfont_bgc = "black"
xtt_col = "black" xtt_col = "black"
@ -317,7 +318,7 @@ LogoPlotMSA <- function(unified_msa
#------------------- #-------------------
# Mutant logo plot # Mutant logo plot
#------------------- #-------------------
p0 = ggseqlogo(plot_mut_edM p0 = ggplot() + geom_logo(plot_mut_edM
, method = msa_method , method = msa_method
, col_scheme = my_logo_col , col_scheme = my_logo_col
, seq_type = 'auto') + , seq_type = 'auto') +
@ -335,6 +336,7 @@ LogoPlotMSA <- function(unified_msa
, vjust = 0.4 , vjust = 0.4
, colour = xfont_bgc) , colour = xfont_bgc)
#, axis.text.y = element_blank() #, axis.text.y = element_blank()
, axis.ticks=element_blank()
, axis.text.y = element_text(size = y_ats , axis.text.y = element_text(size = y_ats
, angle = y_tangle , angle = y_tangle
, hjust = 1 , hjust = 1
@ -344,12 +346,24 @@ LogoPlotMSA <- function(unified_msa
, colour = xtt_col) , colour = xtt_col)
, axis.title.y = element_text(size = y_tts , axis.title.y = element_text(size = y_tts
, colour = ytt_col) , 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) xlab(x_lab_mut)
if (missing(plot_positions)){ if (missing(plot_positions)){
ed_mut_logo_P = p0 + 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 scale_x_discrete(breaks = msa_pos
, expand = c(x_axis_offset, 0) , expand = c(x_axis_offset, 0)
, labels = msa_pos , labels = msa_pos
@ -357,38 +371,25 @@ LogoPlotMSA <- function(unified_msa
}else{ }else{
ed_mut_logo_P = p0 + 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 scale_x_discrete(breaks = i_extract
, expand = c(x_axis_offset_filtered, 0) , expand = c(x_axis_offset_filtered, 0)
, labels = i_extract , labels = i_extract
, limits = factor(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') cat('\nDone: MSA plot for mutations')
#return(msa_mut_logoP) #### Wild-type MSA: gene_fasta file ####
PlotlogolasL[['ed_mut_logoP']] <- ed_mut_logo_P p1 = ggplot() + geom_logo(plot_wt_edM
#---------------------------------
# Wild-type MSA: gene_fasta file
#---------------------------------
p1 = ggseqlogo(plot_wt_edM
#, facet = "grid" #, facet = "grid"
, method = msa_method , method = msa_method
, col_scheme = my_logo_col , col_scheme = my_logo_col
@ -396,16 +397,11 @@ LogoPlotMSA <- function(unified_msa
theme(legend.position = "none" theme(legend.position = "none"
, legend.direction = leg_dir , legend.direction = leg_dir
#, legend.title = element_blank()
, legend.title = element_text(size = leg_tts , legend.title = element_text(size = leg_tts
, colour = ytt_col) , colour = ytt_col)
, legend.text = element_text(size = leg_ts) , legend.text = element_text(size = leg_ts)
, axis.text.x = element_blank()
, axis.text.x = element_text(size = x_ats , axis.ticks=element_blank()
, angle = x_tangle
, hjust = 1
, vjust = 0.4
, colour = xfont_bgc)
, axis.text.y = element_blank() , axis.text.y = element_blank()
, axis.title.x = element_text(size = x_tts , axis.title.x = element_text(size = x_tts
@ -413,51 +409,34 @@ LogoPlotMSA <- function(unified_msa
, axis.title.y = element_text(size = y_tts , axis.title.y = element_text(size = y_tts
, colour = ytt_col) , 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") ylab("") + xlab("Wild-type position")
if (missing(plot_positions)){ if (missing(plot_positions)){
# No y-axis needed # No y-axis needed
ed_wt_logo_P = p1 + ed_wt_logo_P = p1# +
scale_x_discrete(breaks = wt_pos
, expand = c(x_axis_offset, 0)
, labels = wt_pos
, limits = factor(wt_pos))
}else{ }else{
ed_wt_logo_P = p1 + ed_wt_logo_P = p1 +
scale_x_discrete(breaks = i_extract scale_x_discrete(expand = c(0, 0),
, expand = c(x_axis_offset_filtered, 0) breaks = i_extract,
, labels = i_extract #labels = i_extract,
, limits = factor(i_extract)) limits = factor(i_extract)
)
} }
cowplot::plot_grid(ed_mut_logo_P
cat('\nDone: MSA plot for WT') , ed_wt_logo_P
#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 , nrow = 2
, align = "v" , align = "v"
, axis='lr' #, axis='lr'
, rel_heights = c(3/4, 1/4)) , rel_heights = c(3/4, 1/4))
return(LogoED_comb)
} }
#LogoPlotMSA(unified_msa) #LogoPlotMSA(unified_msa)

View file

@ -48,7 +48,7 @@ site_snp_count_bp <- function (plotdf
#------------------------------------------- #-------------------------------------------
# adding column: snpcount for each position # 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 # from dplyr
plotdf = plotdf %>% plotdf = plotdf %>%
@ -57,21 +57,21 @@ site_snp_count_bp <- function (plotdf
plotdf = as.data.frame(plotdf) plotdf = as.data.frame(plotdf)
class(plotdf) class(plotdf)
nc_change = which(colnames(plotdf) == "n") nc_change = which(colnames(plotdf) == "n")
colnames(plotdf)[nc_change] <- "pos_count" colnames(plotdf)[nc_change] <- "position_count"
class(plotdf) class(plotdf)
# if (all(plotdf$pos_count==plotdf$pos_count_check) ){ # if (all(plotdf$position_count==plotdf$position_count_check) ){
# cat("\nPASS: pos_count column created") # cat("\nPASS: position_count column created")
# plotdf = plotdf[, !colnames(plotdf)%in%c("pos_count_check")] # plotdf = plotdf[, !colnames(plotdf)%in%c("position_count_check")]
# }else{ # }else{
# stop("\nAbort: pos count numbes mismatch from dplyr and data.table") # stop("\nAbort: pos count numbes mismatch from dplyr and data.table")
# } # }
cat("\nCumulative nssnp count\n" cat("\nCumulative nssnp count\n"
, table(plotdf$pos_count)) , table(plotdf$position_count))
# calculating total no. of mutations # calculating total no. of mutations
tot_muts = sum(table(plotdf$pos_count)) tot_muts = sum(table(plotdf$position_count))
# sanity check # sanity check
@ -92,10 +92,10 @@ site_snp_count_bp <- function (plotdf
# creating df: average count of snpcount for each position # creating df: average count of snpcount for each position
# created in earlier step # created in earlier step
#------------------------------------------------------- #-------------------------------------------------------
# use group by on pos_count # use group by on position_count
snpsBYpos_df <- plotdf %>% snpsBYpos_df <- plotdf %>%
dplyr::group_by(eval(parse(text = df_colname))) %>% 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" cat("\nnssnp count per position\n"
, table(snpsBYpos_df$snpsBYpos) , table(snpsBYpos_df$snpsBYpos)