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 ([email protected]), the log + normalized counts (object@data), or the scaled counts ([email protected]). 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 <- [email protected]['DDB_G0267412', ident1] / [email protected]['DDB_G0277853', ident1]
datatoplot2 <- [email protected]['DDB_G0267412', ident2] / [email protected]['DDB_G0277853', ident2]
boxplot(datatoplot1, datatoplot2)
stripchart(datatoplot1,
vertical = TRUE, method = "jitter",
pch = 21, col = "dark green",
add = TRUE, at=1)
stripchart(datatoplot2,
vertical = TRUE, method = "jitter",
pch = 21, col = "dark green",
add = TRUE, at=2)
