# Violinplot of gene expression

I have plotted the log normalized expression of two genes by violonplot for 4 clusters. I have links to my pictures and Seurat object too. My problem is this; in violin plot I can not see the mean or any centennial tendencies so that I don't know if two genes is expressing higher or lower in contrast to each other in each cluster.

violonplot vlnPlot(seuset,c(""))

Actually when I am thinking I see I need a violin or box plot showing the ratio of these genes to each other in each cluster and how nice to have mean also in the plot

Somebody knows how I can do that?

• Please, remember to add the code you use to make it easier to provide the accurate advise to help you. But do you want to see the mean of the cluster or to see the differences of genes between clusters?
– llrs
Jun 15 '18 at 15:26

You want to (1) see the mean for each gene, and also to (2) calculate a ratio of expression levels of two genes, then compare it between clusters.

(1) First, notice that vlnPlot() is deprecated. Use VlnPlot().

You can try using the parameter do.sort=T:

VlnPlot(object=seuset, features.plot=c("DDB_G0267412", "DDB_G0277853"), do.sort=T)

Alternatively, you can return the ggplot2 object and then plot the means etc:

p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)
p


(2) There are a few problems with calculating the ratio of gene expression levels. Firstly, what do you mean by gene expression level and how do you measure it? Surely you can have the read counts, but how do you interpret them? Is 30 counts a high level? Is it much more than 60 counts, or is it roughly the same? The same applies to the calculated ratios and the differences between them, even if we ignore amplification, gene length and other biases. And it is very hard to interpret ratios if the reference can also change. The problem is somewhat similar to RT-qPCR, where people use a set of reference genes whose expression has previously been shown to be invariant under the conditions.

In theory, you could use the raw counts (object@raw.data), the log + normalized counts (object@data), or the scaled counts (object@scale.data). The raw counts are biased by sequencing-depth, and the ratio of log or scaled values are not easily interpretable or intuitive. Therefore the library-size normalized (non-log) values seem to be the best.

Additionally, you could calculate the ratio of two genes either (a) for each cell (paired), or (b) for each group. (a) is problematic, because of the zero values: you will have many NaN and Inf values, which cannot be removed without biasing the data.

Here are some functions for retrieving and plotting data from the object:

ident1 <- WhichCells(seuset, ident=1)
ident2 <- WhichCells(seuset, ident=2)

FetchData(seuset, vars.all='DDB_G0267412', cells.use=ident1, use.raw=T)

datatoplot1 <- seuset@raw.data['DDB_G0267412', ident1] / seuset@raw.data['DDB_G0277853', ident1]
datatoplot2 <- seuset@raw.data['DDB_G0267412', ident2] / seuset@raw.data['DDB_G0277853', ident2]

boxplot(datatoplot1, datatoplot2)

stripchart(datatoplot1,
vertical = TRUE, method = "jitter",
pch = 21, col = "dark green",