Anova + visualization

[france-hub]

Hello Belinda
Thanks for this package!
I need some help to understand if I am using your package properly and on how to visualize the data. I'll put here my attempt to visualize the results in case other users find it useful.

I have single cell RNA-seq data, 3 conditions labeled as group_id: Responders (Res), Non responders (NonRes) and healthy donors (HD); paired: for Res and NonRes I have the timepoint pre-treatment and post-treatment; 3 batches. My seurat object is called CD8.
First I used:

props <- getTransformedProps(CD8$clusters, CD8$sample_id, transform="logit")
plotCellTypePropsMeanVar(props$Counts)

That gave me this:

Then, I decided to use the function propeller.anova using the argument robust = TRUE

condition <- c(rep("HD",2), rep("NonRes",2), rep("Res", 4), rep("NonRes", 2), 
               rep("Res", 2))
pair <- rep(c(1,2,3,4,5,6), each = 2)
batch <- c(rep(1,2), rep(2,2), rep(3,2), rep(2,2), rep(1,2), rep(2,2)) #assign the batches to each sample_id
design <- model.matrix(~0 + condition + pair + batch)
df.anova <- propeller.anova(prop.list=props, design=design, coef = c(1,2,3), 
                            robust=TRUE, trend=FALSE, sort=TRUE) %>% as.data.frame() #coeff 1,2,3 conditions of interest

I got that the FDR of one of my subsets (StL) is significant. How to plot the results? I used a heatmap:

df.anova <- df.anova %>% set_colnames(c("HD", "NonRes", "Res", "Fstat", "FDR"))
fdr = 0.05
s <- factor(
  ifelse(df.anova$FDR < fdr, "yes", "no"), 
  levels = c("no", "yes"))
fdr_pal <- c("blue", "red3")
names(fdr_pal) <- levels(s)
F_pal <- c("blue",  "red3")
F_lims <- range(df.anova$Fstat)
F_brks <- c(0, f_lims[2])
Fstat = df.anova$Fstat
anno_cols <- list(Fstat = colorRamp2(F_brks, F_pal))
anno_cols$significant <- fdr_pal

right_anno <- rowAnnotation(
  Fstat = Fstat,
  significant = s,
  "foo" = row_anno_text(
    scientific(df.anova$Fstat, 2),
    gp = gpar(fontsize = 8)),
  col = anno_cols,
  gp = gpar(col = "white"),
  show_annotation_name = FALSE,
  simple_anno_size = unit(4, "mm"))

mat <- as.matrix(df.anova[,1:3])

Heatmap(mat, name = "Transformed proportions",
        right_annotation = right_anno,
        cluster_columns = FALSE,
        cluster_rows = FALSE,
        heatmap_legend_param = list(title_gp = gpar(
          fontsize = 10, fontface = "bold", lineheight = 0.8)))

and got this:

Does this work? Would you try other models for these data? Is the heatmap the way to go?

Thanks!
Francesco