• allow different colour code for fc > 0 and fc < 0
• add marker outline
• reorder legend levels
• allow label to be colour coded

This version requires the qval() function (used within deseq_qvalue(), line 126 of deseq_polar.R) that I assume must come from the charm package.
However, the charm package is not included in the DESCRIPTION file of this package, neither in Imports nor in Suggests.
Plus, the charm package seems to be unavailable for the latest version of R.

Thanks for the cool package! I am a little confused about the calculation of the r_zscore & axis scale on 2D radial plots- I saw somewhere in the documentation that it represents -log10 pvalue but it doesn't match the p-value of any of the comparisons. how is the r_zscore calculated?

Example where comparisons include 'Lining' and 'SubLining' lead to error in significance grouping

Not sure if this was reported but I think the fc_cutoff argument for polar_coords() doesn't seem to work?

NK_polar_limma <- polar_coords(sampledata = z_meta,
contrast = "group",
pvalues = Pvals_limma_NK2,
expression = z_count,
p_col_suffix = "pvalue",
padj_col_suffix = "padj",
multi_group_prefix = "Overall",
non_sig_name = "Not Significant",
fc_col_suffix = 'log2FoldChange',
significance_cutoff = 0.05,
label_column = NULL,
fc_cutoff = 5,
cutoff_criteria = "padj")

hovertext = "paste(label,
'\nCD56Bright vs CD56DimCD57Neg padj:', format(CD56Bright_CD56DimCD57Neg_padj, digits=3),
'\nCD56DimCD57Neg vs CD56DimCD57Pos padj:', format(CD56DimCD57Neg_CD56DimCD57Pos_padj, digits=3),
'\nCD56DimCD57Pos vs CD56Bright padj:', format(CD56DimCD57Pos_CD56Bright_padj, digits=3),
'\nCD56Bright vs CD56DimCD57Neg logFC:', format(CD56Bright_CD56DimCD57Neg_logFC, digits=3),
'\nCD56DimCD57Neg vs CD56DimCD57Pos logFC:', format(CD56DimCD57Neg_CD56DimCD57Pos_logFC, digits=3),
'\nCD56DimCD57Pos vs CD56Bright logFC:', format(CD56DimCD57Pos_CD56Bright_logFC, digits=3))"

volcano3D(NK_polar_limma,
hover=hovertext)

I set the fc_cutoff to 5 but genes with fc less than 5 still comes up as significant.

Love the tool. Thanks in advance.

I get error message,
Error in boxplot_trio(polar, value = "row.names") :
value/gene is not in rownames(expression)

Have you changed the value, I don't see COBL anywhere in the structure and am not sure what structure to point to?

This is my polar str (shortened). Am I not supposed to be passing in polar. I am using the latest version of Volcano3d on M1 mac.
Using the latest R, within Rstudio. Other commands so far work well.

new("volc3d", df = list(scaled = structure(list(KO_D7 = c(-1.05648244015833,

• 11.4184810397571, 5.99822138275276, 4.45723462152403, 3.50925761487716,
  0.727161529954877, 8.30547288694249, 18.7912362592882, 4.28282258335944
  ), col = c("red", "cyan", "green3", "grey60", "cyan", "red",
  "red", "cyan", "cyan", "red", "cyan", "red", "cyan", "red",
  "grey60", "cyan", "cyan", "cyan", "cyan", "red", "cyan",
  5L, 2L, 2L, 2L, 1L, 5L, 6L, 5L, 5L, 2L, 2L, 2L, 1L, 5L, 2L,
  2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 5L, 2L, 5L, 4L, 5L,
  5L, 1L, 2L, 5L, 1L, 5L, 1L, 5L, 5L, 5L, 5L, 1L), levels = c("ns",
  "WT_D0+", "WT_D0+WT_D7+", "WT_D7+", "K+WT_D7+", "K+", "K+WT_D0+"
  ), class = "factor")), class = "data.frame", row.names = c("chrX.100402221.100442717.-.Awat2",
  "chrX.100494290.100516125.+.Igbp1", "chrX.100622905.100625907.-.Pdzd11",
  "chrX.99975605.100400760.+.Eda")), type = "voom_polar"), outcome = structure(c(2L,
  2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L), levels = c("KO_D7", "WT_D0",
  "WT_D7"), class = "factor"), data = structure(c(3, 3.70043971814109,
  3, 3.16992500144231, 3.4594316186373, 2.32192809488736, 2, 2,
  2.32192809488736, 6.90689059560852, 7.20945336562895, 6.4262647547021,
  6.02236781302845, 8.36632221424582, 9.08480838780436, 8.74819284958946,
  4.95419631038687), dim = c(9L, 418L), dimnames = list(c("sample1",
  "sample2", "sample3", "sample4", "sample5", "sample6", "sample7",
  "sample8", "sample9"), c("chrX.100402221.100442717.-.Awat2",
  "chrX.100494290.100516125.+.Igbp1", "chrX.100622905.100625907.-.Pdzd11",
  "chrX.99136060.99149022.+.Efnb1", "chrX.99465733.99470735.-.Pja1",
  "chrX.99975605.100400760.+.Eda"))), pvals = structure(c(0.00114784428547215,
  0.0360467209510449, 0.00580588842816712, 0.0784652760874108,
  1.04301150682251e-14, 1.52336054334604e-10, 0.5307512359096), dim = c(418L,
  4L), dimnames = list(NULL, c("Pvals_overall", "", "", ""))),
  padj = structure(c(0.00164879350971602, 0.0444187468827374,
  0.0078539202685238, 0.0931775153538003, 4.40085120120768e-07,
  2.37806623555533e-13, 8.14045790350539e-10, NA), dim = c(418L,
  4L), dimnames = list(NULL, c("LRTpadj", "", "", ""))), pcutoff = 0.05,
  scheme = c("grey60", "red", "gold2", "green3", "cyan", "blue",
  "purple"), labs = c("ns", "WT_D0+", "WT_D0+WT_D7+", "WT_D7+",
  "K+WT_D7+", "K+", "K+WT_D0+"))

Hi guys,

I am using your volcano3D package to plot my own data. However, when I do the layout as you did on your website, I can't make my volcano plot as large as yours. I tried to set the height from 720 to 850 but still not large enough. Could you give me some advice on it?

Thanks!

Hi,
Thank you for this useful package.

I am a novice to R and was wondering how could I influence the division into significance groups? As I understand, this is somehow determined by pairwise comparisons to the minimum group mean expression as a reference?
There are genes that are classified and colored as enriched in one group, although their angel sets them far away from the relevant axis (e.g. green genes which reside far from the AdultA axis and near the L3A axis) - Is there a way I could control/change that?
For example, in my data, I have 3 groups: L3A, Adult A, and DevN. The gene 'Ama' is classified as DevN+ L3A+, but the expression is much higher in the L3A group compared to the two other groups (although the difference between DevN and AdultA might also be significant, AdultA being the minimum group in this case).

Thanks for your time.

This is my code:

metadata <- read.csv("metadata.csv")
pvalues <- read.csv("pvalues.csv")
expression = read.csv("expression.csv")
genes <- c(pvalues[,1])
rownames(pvalues) <- genes
pvalues <- pvalues[,2:10]
rownames(expression) = genes
expression <- expression[,2:55]
check.rows <- apply(expression, 1, sd)
expression <- expression[check.rows>0, ]
pvalues <- pvalues[check.rows>0, ]

polar_d <- polar_coords(sampledata = metadata,
contrast = "Pathotype",
pvalues = pvalues,
expression = expression,
p_col_suffix = "pvalue",
padj_col_suffix = "padj",
fc_col_suffix = "log2FoldChange",
non_sig_name = "Not Significant",
significance_cutoff = 0.05,
label_column = NULL,
cutoff_criteria = "padj",
fc_cutoff = 1)

Radial ggplot:

p1 <- radial_ggplot(polar = polar_d,
fc_or_zscore = 'fc',
label_rows = c("Ama"),
arrow_length = 1.1,
marker_size = 2.5,
legend_size = 10) +
theme(legend.position = "right")
p1

when I was ploting volcano3D using radial_ggplot, I found some gene labels at the left edge of plot showed imcompletely. Can any parameters change this?

Hi Kat!

There is a minor bug on the color ordering of the volcano3D.
Look at this example:

If you look look at the legend you'll see categories are ordered alphabetically while "Not significant" should rather go last.
This impacts the volcano visualization as if you tip the image the gray dots will overlay the green ones:

so you don't get what you would see with radial_plotly:

Please let me know if this is not clear!

When I plot the volcano3D, I observed the legend was away from the main figure. Is there any way to reduce their distance? Below are my codes:

radial_ggplot(polar = polar,
marker_size = 3,
label_size = 8,
text_coords = 0.2,
arrow_length = 1,
marker_alpha = 1,
axis_title_size = 12,
axis_label_size = 8,
legend_size = 15) +
theme(legend.position = "right")+
coord_fixed(ratio = 1, xlim = c(-3, 3),ylim = c(-3, 3))

Idea to create a function to extract information about specific significant populations from the polar object.

Hi,

I love this 3D volcano plot and want to show my RNAseq study result. I tried to set the log2 fold change in the deseq_polar function, but I only see 'pcutoff' is available in this function as a threshold, what can I do to change the default foldchange value to make the highlighted gene of interest to be significant (those genes past the cutoff, not the default FC cutoff)?

Best,
Shuai

Hi!
First, thanks for the awesome package.
I've just tested it, and replicated the vignette you explained here without any problems.
https://towardsdatascience.com/interactive-differential-expression-analysis-with-volcano3d-cf6c37b5aa6a

Then I used my data, with the disclaimer that I tested it on methylation data, and everything is working fine except the volcano3D function.

I get the following message:

Error in data.frame(x = rep(c(max(r_breaks) * cospi(0:200/100), NA), n_z_breaks),  : 
  arguments imply differing number of rows: 0, 404, 1
In addition: Warning message:
In max(z_breaks) : no non-missing arguments to max; returning -Inf

Could it be because of my data-specific values?
instead of fold change I'm using estimate which then ranges usually from 0.0 to 0.35.
and the values in the matrix range from 0.0 to 1.0.

Thank you very much!

Jasko

Hi Kat,

Why not to turn polar_coords' cutoff_criteria to "padj" by default?
I know this is not a bug really, but it's something Myles and I found a bit weird/misleading.

Wouldn't be more intuitive to set:

 significance_cutoff = 0.05,
 fc_cutoff = 0,
 cutoff_criteria = "padj"

by default?
I'm saying this just because while I was using the function I took for granted pvalues where already adjusted but then the volcano plot wasn't look like I was expecting and it took me a while to find out this was the problem.

Hi,

First of all, thanks for an excellent package.

I'm struggling a bit with getting a shiny app set up. While I can produce both volcano3d and radial_plotly plots, I can't get the app to output boxplots based on clicks and I also get this warning:

Warning: The 'plotly_click' event tied a source ID of 'heatmap' is not registered. In order to obtain this event data, please add event_register(p, 'plotly_click') to the plot (p) that you wish to obtain event data from.

Any help would be much appreciated.

Cheers,
Lasse

I'm currently searching for a way of using volcano3D on single cell data analyzed with Seurat. Volcano3D looks like a fun and interesting method for DE analysis representation but I don't know which is the right way to do it from a single cell seurat object? EnhancedVolcano works like a charm with DE analysis by Seurat and I wanted to know if there is a similar way to get volcano3d working with this kind of data.

There seems to be an issue with colour_code_labels=TRUE for the radial_plotly function where the labels that are assigned "green3" and "gold2" (F+ and M+F+, respectively) are not colour coded accordingly but are grey.

The colour_code_labels=TRUE works fine for ggplot2.

It may be because plotly does not recognise all colours as described here but I imagine it should recognise given points are coloured appropriately? https://community.rstudio.com/t/plotly-colors-setting-dont-work-in-subplot/72073

Hello,
I'm a bit confused about which p-values are used for the pcutoff value set in polar_coords. When I change the cutoff I get more/fewer genes accordingly but the actual value doesn't match the returned one-way or pairwise p-values (for example, if I set the p-cutoff to 0.01 the gene with the highest returned one-way pvalue is 7E-05 even though there are many genes that come back a nonsignificant that have a p-value between 7E-05 and 0.01)
Thanks

Nice job about the R package !
I have ran all the pipeline ,well, no error shown .
But ,i cannot understand the radial plot result.
Well, here my IMQ+ DEGs were more closed to the SHP099 axis ???
the Normal+ DEGs were more closed to SHP099 axis ???
I have no idea to understand this ,could you help me ,please .