FeaturePlot issues and suggestions in Seurat3

[jeremycfd]

I've noticed a few problems with the Seurat3 FeaturePlot function that I wanted to provide feedback on:

1. When using the split.by option, FeaturePlot correctly separates according to the factor of interest; however, it seems that each sub-plot scales the color (corresponding to feature expression) separately. This only became obvious to me when I plotted a feature that is not expressed in any cells in one of the sub-plots. The default behavior with facet_wrap() is to maintain the color scale across all subplots, and I think that should probably be the default here as well (and also provide the option to display that scale if possible); otherwise, I fear a lot of people will be inadvertently mislead by these plots.

2. The ncol option in FeaturePlot doesn't seem to work as far as I can tell, or perhaps just not how I envision it working. Setting it to any value doesn't seem to change the resulting plot. I think most people would appreciate being able to control this in a way similar to nrow and ncol options in ggplot2's facet_wrap(), at least if there is only a single feature being plotted, which currently can be done by tacking on the facet_wrap() to the FeaturePlot object.

3. I understand if you don't want to implement this as a feature, but do you have any suggestions on how to edit the font size for different aspects of the FeaturePlots through ggplot() options? I've dug through the FP object but haven't been able to find where those sorts of things are defined.

Thanks!

[andrewwbutler]

Hi,

Thanks for the feedback!

1. I agree. The split.by option still needs a bit of work regarding scaling/legends. Will post back here when we have any updates.
2. I'm not sure what issue you're having with ncol. Using pbmc_small, this example shows changing from 3 to 1 column. Could you provide an example where it breaks?

FeaturePlot(object = pbmc_small, features = c("MS4A1", "LYZ", "CD3D"), ncol = 3)
FeaturePlot(object = pbmc_small, features = c("MS4A1", "LYZ", "CD3D"), ncol = 1)

3. When possible, all plots that Seurat v3 produces are made using the ggplot2. Things like font size are controlled through the theme framework. https://ggplot2.tidyverse.org/reference/theme.html So for example, to change the axis text size:

FeaturePlot(object = pbmc_small, features = "MS4A1") + theme(axis.text = element_text(size = 20))

[bumproo]

Hi, I'm writing to followup on point1 from the original message because I have been struggling with the same issue. I am working on making FeaturePlots for a single gene split.by three experimental conditions. In one of my conditions the gene is not expressed so the scale in that condition is just 0. Unfortunately the color scale for that condition sets 0 as the max color in the plot.

the code looks like this: ############# p <- FeaturePlot(seuobj,reduction="umap",features=features, min.cutoff=0,max.cutoff=6,cols=c("grey","red"), split.by = "Condition", combine = FALSE)

cowplot::plot_grid(plotlist = p) ############# Can anyone suggest a way to enforce the same color scale (0 to 6) on all 3 panels of the plot? thanks!

[jeremycfd]

I ended up exporting the data from the Seurat object and plotting using ggplot2

[vertesy]

The ncol still does not work when

FeaturePlot(combined.obj, features = 'TOP2A', ncol = 2 , split.by = 'experiment');

...where experiment splits the data into 4 subsets.

[lukedow]

Hi, I was wondering if there has been any resolution of the scaling and legend issues in FeaturePlot when using "split.by"? I really like the whole package and the ability to show expression across experimental conditions, but has anyone figured out how to ensure the plots are scaled together (without exporting to ggplot2)? Thanks!

[alexander-7]

Hi, first of all @satijalab thanks a lot for the great package (Seurat v3), which I am using a lot! I also really like the functionality of the "split.by" option of the FeaturePlot. However, due to the problems with the scaling/legend, I ended up subsetting the data after RunUMAP and use the same embedding for different subsets to plot the expression. However, since you manually need to adapt the scale for every plot, this is quite tedious... Are there any news regarding the scaling/legend issue referred to above? Thanks a lot!

[alexander-7]

I wrote an R-Function as workaround for the scaling issue of the 'split.by' option of Feature Plot. Hope this is helpful:

`plot_seperate_features <- function(object, gene, limit = NULL,ident=NULL,pt.size=0.5,revers.order = F,numcol=NULL, colors =c('#DCD9DE','#870052'), ...){

if(is.null(ident)){ Idents(object) <- 'combined' } else{ Idents(object) <- object[[ident]] } p <- vector("list", length = (length(gene)*length(levels(Idents(object))))) for(g in 1:(length(gene))){ if(is.null(limit)==T){ for(i in 1:length(levels(Idents(object)))){ j <- i+((g-1)*length(levels(Idents(object)))) p[[j]] <- FeaturePlot(object, cells = WhichCells(object,idents = levels(Idents(object))[i]), features = gene[g],order = T,cols = colors, pt.size = pt.size, sort.cell = T,...) + scale_color_gradient(low = colors[1], high = colors[2], limits = c(0,ceiling(max(object@assays$RNA@data[gene[g],]))),oob = squish) + ggtitle(paste(levels(Idents(object))[i],gene[g],sep=", ")) } } else{ for(i in 1:length(levels(Idents(object)))){ j <- i+((g-1)*length(levels(Idents(object)))) p[[j]] <- FeaturePlot(object, cells = WhichCells(object,idents = levels(Idents(object))[i]), features = gene[g],order = T,cols = colors, pt.size = pt.size, sort.cell = T,...) + scale_color_gradient(low = colors[1], high = colors[2], limits = limit,oob = squish) + ggtitle(paste(levels(Idents(object))[i],gene[g],sep=", ")) } } } #correct axis in all plots y_axis_max <- 1:length(p) y_axis_min <- 1:length(p) x_axis_max <- 1:length(p) x_axis_min <- 1:length(p) for(i in 1:length(p)){ x_axis_max[i] <- max(layer_scales(p[[i]])$x$range$range) x_axis_min[i] <- min(layer_scales(p[[i]])$x$range$range) y_axis_max[i] <- max(layer_scales(p[[i]])$y$range$range) y_axis_min[i] <- min(layer_scales(p[[i]])$y$range$range) } for(i in 1:length(p)){ p[[i]] <- p[[i]]+coord_equal(xlim=c(min(x_axis_min),max(x_axis_max)),ylim=c(min(y_axis_min),max(y_axis_max))) } if(revers.order==T) p <- rev(p) gridExtra::grid.arrange(grobs = p,ncol=numcol) }`

[arielyh]

Hi @andrewwbutler, The "ncol = " setting works fine when there is no "split.by" option. But seems like the "split.by" suppress the "ncol".

[hemadelon]

I wrote an R-Function as workaround for the scaling issue of the 'split.by' option of Feature Plot. Hope this is helpful:

code

This is really great, thanks!

[greensii]

a simpler workaround for the issue with consistent color scaling when using FeaturePlot + split.by is to use combine=FALSE and patchwork, then add theme/scale elements at the end that will be applied to all plots:

patchwork::wrap_plots( FeaturePlot( seurat_obj, features=gene, split.by = "treatment", combine=FALSE)) & theme_minimal() & scale_color_gradient(low = "gray", high = "blue", limits = c(0,5))

[ghar1821]

Following up on point 1 as I recently have had to deal with the disparate colour scales across different condition/group. I want the plots for each gene to have the same colour scale across all conditions for comparison. Setting one global min max limit for all plots is not ideal as I want to compare each gene by itself across all conditions rather than all genes across all conditions if that makes sense.

So I ended up subsetting the data for each condition, calculate the min/max for each gene, draw the feature plot for each gene and condition, and override the colour scaling using limits. Then use patchwork to organise everything into a single plot.

It's rather painful as my dataset is huge and I have just enough RAM to store each condition after removing a bunch of variables. Not ideal, but it got the job done.

My workaround below (timepoint is my condition). Hopefully it is useful for future users.

samples <- c("D0", "D1", "D2")
# Store data per time-point
dat_per_tp <- lapply(samples, function(samp) {
    dat_tp <- subset(dat_sub, subset = timepoint == samp)
    return(dat_tp)   
})
names(dat_per_tp) <- samples

# The genes to plot.
to_plots <- c("Cd3", "Pdgfra", "Cd14")

plots <- lapply(to_plots, function(to_plot) {
    
    # Compute new colour scale
    assay_dat <- GetAssayData(subset(x = dat_sub, features = c(to_plot)))
    
    min_cnt <- min(assay_dat)
    max_cnt <- max(assay_dat)
    
    p <- lapply(samples, function(samp) {
        # Used viridis colour scheme. Change to normal scale_color_gradient if you like.
        plt <- FeaturePlot(dat_per_tp[[samp]], to_plot, order = TRUE) + 
            ggtitle(paste(samp, title)) +
            scale_color_viridis_c(limits = c(min_cnt, max_cnt))
            
        return(plt)
    })
    names(p) <- samples
    return(p)
})
names(plots) <- to_plots

# Flatten all plots into 1 giant list
plots_flat <- do.call(c, plots)

# Stitch all the plots into 1.
plt <- patchwork::wrap_plots(plots_flat, ncol=4)

# Save
ggsave("plot_split_by_tp.pdf", plot = plt, width = 24, height = 60, limitsize = FALSE)

This is a minor thing, but it will be nice for FeaturePlot to automatically support the colour schemes available in the viridis package as they are rather popular. Yes you can change the colour scheme post FeaturePlot, but this will be painful if the scales have all been set up beforehand. Alternatively, add a small section in the visualisation vignette about how to change colour scheme or just point to an existing tutorial like this one https://www.r-graph-gallery.com/ggplot2-color.html so it is easy for users to find.

[Elo-mars]

a simpler workaround for the issue with consistent color scaling when using FeaturePlot + split.by is to use combine=FALSE and patchwork, then add theme/scale elements at the end that will be applied to all plots:

patchwork::wrap_plots( FeaturePlot( seurat_obj, features=gene, split.by = "treatment", combine=FALSE)) & theme_minimal() & scale_color_gradient(low = "gray", high = "blue", limits = c(0,5))

exactly what I needed! --> in the umap coordinates with brewer palette

patchwork::wrap_plots(FeaturePlot(object = seurat_ex,slot="data", features = c(pathway),split.by ="Layer", reduction = "umap", pt.size = 0.5, keep.scale = "feature", combine=FALSE)) & theme_minimal() & labs(title = pathway) & scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))