:- lib(b_real:options_rvar_rmv/2).
:- lib(stoics_lib:at_con/3).
:- lib(stoics_lib:lexi/2).
:- lib(stoics_lib:colour_hex/2).
bio_volcano_plot_defaults( Args, Defs ) :-
bio_analytics:bio_diffex_defaults( Args, DfxDefs ),
Defs = [
clr_down(brandeisblue),
clr_inv(cadetgrey),
clr_line(bazaar),
clr_up(cadmiumred),
ext(pdf),
debug(true),
legend_cnm(regulation),
legend_down("down"),
legend_inv("invariant"),
legend_title(_),
legend_up("up"),
% lim_x_min(_),
% lim_x_max(_),
lim_y_min(0),
% lim_y_max(_),
plot_file(_),
rvar_cv(bvp_cv),
rvar_df(bvp_df),
rvar_gt(bvp_gt),
rvar_rmv(true),
stem(volc_plot),
theme(classic)
| DfxDefs
].
/** bio_volcano_plot( +Mtx, +Opts ).
Create a volcano plot of differential expressed values (x) against p-values (y).
Uses ggplot2 and bio_analytics:bio_diffex/4.
The predicate creates an R data.frame and creates a ggplot2 object from which
typically output files are produced.
Opts
* clr_down(ClrDw=brandeisblue)
bio_colour_hex/2 compatible colour for down regulated entries
* clr_inv(ClrInv=cadetgrey)
colour for invariant entries
* clr_line(ClrLn=bazaar)
colour for cut-off lines
* clr_up(ClrUp=cadmiumred)
colour for up regulated entries
* dir(Dir)
as understood by os_dir_stem_ext/2- see Odir, below too
* debug(Dbg=true)
informational, progress messages
* ext(Ext=pdf)
extension that defines the type of output (passed to os_dir_stem_ext/2)
* legend_cnm(LgCnm='regulation')
the legend frame column name (also appears as legend title, if none is given)
* legend_down(LegUp="down")
token for up-regulation
* legend_inv(LegUp="invariant")
token for invariants
* legend_title(LgTitle)
by default LgCnm is used.
* legend_up(LegUp="up")
token for up-regulation
* lim_x_max(LimXmax)
if LimXmin is given this defaults to -LimXmin. By default ggplot2 decides this.
* lim_x_min(LimXmin)
if LimXmax is given this defaults to -LimXmax. By default ggplot2 decides this.
* lim_y_max(LimXmax)
By default ggplot2 decides this.
* lim_y_min(LimY=0)
this is not used unless LimYmax is given
* odir(Odir)
as understood by os_dir_stem_ext/2- preferred to Dir above
* plot_file(File)
returns the output file
* rvar_cv(Rvar=bvp_cv)
R var for colors vector
* rvar_df(Rvar=bvp_df)
R variable for the 3 column data frame
* rvar_gt(Rvar=bvp_gt)
R variable holding the plot term
* rvar_rmv(RRmv=true)
remove the three R variables from above
* stem(Stem=volc_plot)
stem for output, set to =|false|= if no output is wanted
* theme(Theme=classic)
theme_Theme should be a ggplot() theme
Opts are passed to bio_diffex/4. This predicate also pinches the defaults from there.
Examples, fixme: use iris dataset.
==
?- bio_volcano_plot([]).
==
@author nicos angelopoulos
@version 0.1 2022/12/15
@see bio_diffex/4
@see os_dir_stem_ext/2
@tbd add args for influencing of the gg term and the output-to-file call
*/
bio_volcano_plot( Mtx, Args ) :-
Self = bio_volcano_plot,
<- library("ggplot2"),
options_append( Self, Args, Opts ),
% get the up-down regulated row indices
bio_diffex( Mtx, _, _, [which(dx(Up,Dw))|Opts] ),
debuc( Self, length, [up,down]/[Up,Dw] ),
% construct the data.frame
options( exp_ev_cnm(EvCnm), Opts ),
options( exp_pv_cnm(PvCnm), Opts ),
mtx_column( Mtx, EvCnm, Evs ),
mtx_column( Mtx, PvCnm, Pvs ),
options( legend_cnm(LgCnm), Opts ),
options( legend_inv(LgIv), Opts ),
options( [rvar_cv(Rcv),rvar_df(Rdf),rvar_gt(Rgt)], Opts ),
Rdf <- 'data.frame'(EvCnm=Evs,PvCnm=Pvs,LgCnm=LgIv),
options( legend_up(LgUp), Opts ),
options( legend_down(LgDw), Opts ),
lexi( LgCnm, +(LgCnmAtm) ),
Rdf$LgCnmAtm <- replace( Rdf$LgCnmAtm, Up, LgUp ),
Rdf$LgCnmAtm <- replace( Rdf$LgCnmAtm, Dw, LgDw ),
Rdf$LgCnmAtm <- factor(Rdf$LgCnmAtm, levels=c(LgUp,LgIv,LgDw)),
% build the plot term
options( theme(GGthemeTkn), Opts ),
atom_concat( theme_, GGthemeTkn, GGthemeNm ),
atom_concat( GGthemeNm, '()', GGthemeTerm ),
Rgt <- ggplot(data=Rdf, aes(x=EvCnm, y= -log10(PvCnm), col=LgCnmAtm)) + geom_point() + GGthemeTerm,
options( [clr_down(ClrDw),clr_inv(ClrIv),clr_line(ClrLn),clr_up(ClrUp)], Opts ),
maplist( colour_hex, [ClrDw,ClrIv,ClrLn,ClrUp], [HexDw,HexIv,HexLn,HexUp] ),
options( [exp_ev_cut_let(CutLw),exp_ev_cut_get(CutUp),exp_pv_cut(CutPv)], Opts ),
Rgt <- Rgt + geom_vline(xintercept=c(CutLw,CutUp), col=HexLn)
+ geom_hline(yintercept= -log10(CutPv), col=HexLn),
bio_volcano_gg_limits_term( Rgt, Opts ),
Rcv <- c(HexUp,HexIv,HexDw),
names(Rcv) <- c(+LgUp,+LgIv,+LgDw),
options( legend_title(LgTtl), Opts ),
( var(LgTtl) -> LgTtl = LgCnm; true ),
Rgt <- Rgt + scale_color_manual(+LgTtl,values=Rcv),
% options( plot_save_width(Width), Opts ),
bio_volcano_plot_save( Opts ),
% fixme: keep rvars ?
options_rvar_rmv( [Rcv,Rdf,Rgt], Opts ),
debuc( Self, end, true ).
bio_volcano_plot_save( Opts ) :-
options( stem(false), Opts ),
!.
bio_volcano_plot_save( Opts ) :-
os_dir_stem_ext( File, Opts ),
options( plot_file(File), Opts ),
% <- ggsave( +File, width=Width )
% fixme: allow for arbitary = pairs to be passed ?
<- ggsave( +File ).
bio_volcano_gg_limits_term( Rv, Opts ) :-
bio_volcano_gg_limits_term_x( Xt, Opts ),
bio_volcano_gg_limits_term_y( Yt, Opts ),
bio_volcano_gg_limits_term_conc( Xt, Yt, Rv ).
bio_volcano_gg_limits_term_conc( null, null, _Rv ) :- !.
bio_volcano_gg_limits_term_conc( null, Yt, Rv ) :- !,
Rv <- Rv + coord_cartesian(ylim=Yt).
bio_volcano_gg_limits_term_conc( Xt, null, Rv ) :- !,
Rv <- Rv + coord_cartesian(xlim=Xt).
bio_volcano_gg_limits_term_conc( Xt, Yt, Rv ) :-
Rv <- Rv + coord_cartesian(xlim=Xt, ylim=Yt).
bio_volcano_gg_limits_term_x( Xt, Opts ) :-
bio_volcano_gg_limits_term_x_given( Low, High, Opts ),
!,
Xt = c(Low,High).
bio_volcano_gg_limits_term_x( null, _Opts ).
bio_volcano_gg_limits_term_x_given( Low, High, Opts ) :-
memberchk( lim_x_max(High), Opts ),
memberchk( lim_x_min(Low), Opts ),
!.
bio_volcano_gg_limits_term_x_given( Low, High, Opts ) :-
memberchk( lim_x_max(High), Opts ),
!,
Low is - High.
bio_volcano_gg_limits_term_x_given( Low, High, Opts ) :-
memberchk( lim_x_min(Low), Opts ),
!,
High is - Low.
bio_volcano_gg_limits_term_y( Yt, Opts ) :-
bio_volcano_gg_limits_term_y_given( Low, High, Opts ),
!,
Yt = c(Low,High).
bio_volcano_gg_limits_term_y( null, _Opts ).
bio_volcano_gg_limits_term_y_given( Low, High, Opts ) :-
memberchk( lim_y_max(High), Opts ),
!,
options( lim_y_min(Low), Opts ).