Question: Hello everyone,

Could anyone help me understand if there is any way I can add p values to violin plots produced by Seurat? Currently, I am trying to plot several genes across all clusters, broken up by responders vs non-responders and this is my code:

vp_case1 <- function (gene_signature, file_name){
plot_case1 <- function(signature){
VlnPlot(pbmc, features = signature,
pt.size = 0.1,
group.by = "Response"
) + stat_compare_means(comparisons = pbmc$Response, label = "p.signif")
}

map(gene_signature, plot_case1) %>% cowplot::plot_grid(plotlist = .)
file_name <- paste0(file_name, "_r.png")
ggsave(file_name, width = 14, height = 8)
}

vp_case1(gene_sig, "gene_sig")

where gene_sig is a vector of a few genes of interest.

However, I get the following error:

Warning messages:
1: Computation failed in `stat_signif()`:
not enough 'y' observations
2: Computation failed in `stat_signif()`:
not enough 'y' observations
.
.
.

How could I do this better?

Thank you!

Answer: Hi,

The problem that you got is related with the fact that you need to provide a list of comparisons to be made to the function `stat_compare_means()`, and you are providing a character/factor variable. So, if you adjust your code to:

vp_case1 <- function(gene_signature, file_name, test_sign, y_max){
plot_case1 <- function(signature){
VlnPlot(pbmc, features = signature,
pt.size = 0.1,
group.by = "Response",
y.max = y_max, # add the y-axis maximum value - otherwise p-value hidden
) + stat_compare_means(comparisons = test_sign, label = "p.signif")
}
map(gene_signature, plot_case1) %>% cowplot::plot_grid(plotlist = .)
file_name <- paste0(file_name, "_r.png")
ggsave(file_name, width = 14, height = 8)
}

In this adjusted code I added 2 extra parameters:

+ `test_sign`: that is a list of a vector of characters. So, in your case you have a `Response` variable that should assume 2 or more levels. Let's say that your `Response` variable can be `CTRL` or `TRT` (control or treatment). Thus `test_sign = list(c("CTRL", "TRT"))`.

+ `y_max`: by default the function `VlnPLot` will adjust the y-axis to the maximum expression value *per* comparison made. What this means is that if you're comparing the gene A between CTRL and TRT and the max. value is 7, the maximum y-axis scale will be 7. When you apply the `stat_compare_means()` function, the p-value will be added above, but since the maximum y-axis value is 7, you'll not see it. Therefore, you can adjust the `y_max` to one more unit, let's say 8, and therefore you'll be available to see the p-value. Though as you can think this is not the best option since often you'll have distinct y maximum values, and therefore setting the same `y_max` will not be the best option. In that case what you can do is just a loop instead a `map()` giving a `y_max` a numeric vector. Let me know if you have problems regarding this issue.

Here is an example:

## Follow the Seurat tutorial until UMAP at: https://satijalab.org/seurat/v3.2/pbmc3k_tutorial.html

## create fake 'Reponse' variable
pbmc@meta.data$Response <- rep(c("CTRL", "TRT"), times = c(1500, 2638-1500))

vp_case1 <- function(gene_signature, file_name, test_sign, y_max){
plot_case1 <- function(signature){
VlnPlot(pbmc, features = signature,
pt.size = 0.1,
group.by = "Response",
y.max = y_max, # add the y-axis maximum value - otherwise p-value hidden
) + stat_compare_means(comparisons = test_sign, label = "p.signif")
}
purrr::map(gene_signature, plot_case1) %>% cowplot::plot_grid(plotlist = .)
file_name <- paste0(file_name, "_r.png")
ggsave(file_name, width = 14, height = 8)
}

gene_sig <- c("NKG7", "PF4")
comparisons <- list(c("CTRL", "TRT"))
vp_case1(gene_signature = gene_sig, file_name = "Desktop/gene_sig", test_sign = comparisons, y_max = 7)

Result:

![enter image description here][1]

[1]: https://i.ibb.co/GvLFPs8/gene-sig-r.png

I hope this answers your question,

António