MOFA model on aIg stimulated BJAB cells

seu_combined_selectsamples <- FindVariableFeatures(seu_combined_selectsamples, selection.method = "mvp", assay = "SCT.RNA", num.bin =20, mean.cutoff = c(0, 5), dispersion.cutoff = c(0.5,10), nfeatures =500)

genes.variable <- seu_combined_selectsamples@assays$SCT.RNA@var.features[-grep("^MT", seu_combined_selectsamples@assays$SCT.RNA@var.features)] # without the mitochondrial genes 

proteins.all <- row.names(seu_combined_selectsamples[["PROT"]])

meta.allcells <- seu_combined_selectsamples@meta.data %>%
  mutate(sample = rownames(seu_combined_selectsamples@meta.data))

## Note, if mofa object was generated before, it will read the generated rds file. (this will speed-up the process of generating this html document if edits are required)

myfiles <- list.files(path="output/", pattern = ".rds$")

if("MOFA_aIg.rds" %in%  myfiles){mofa <- readRDS("output/MOFA_aIg.rds")} else { #If so, read object, else do:
     
      mofa <- create_mofa(list(
        "RNA" = as.matrix( seu_combined_selectsamples@assays$SCT.RNA@scale.data[genes.variable,] ),
        "PROT" = as.matrix( seu_combined_selectsamples@assays$PROT@scale.data[proteins.all,] ))
        )

      # Default settings used (try 15 factors, excludes all non-informative factors)
      data_opts <- get_default_data_options(mofa)
      model_opts <- get_default_model_options(mofa)
      train_opts <- get_default_training_options(mofa)
      train_opts$seed <- 42 # use same seed for reproducibility
      mofa <- prepare_mofa(
        object = mofa,
        data_options = data_opts,
        model_options = model_opts,
        training_options = train_opts
        )

    mofa <- run_mofa(mofa, outfile = "output/MOFA_aIg.hdf5")
    mofa <- run_umap(mofa, factors = c(1:7))
    samples_metadata(mofa) <- meta.allcells
    saveRDS(mofa, file= "output/MOFA_aIg.rds")
  
}


mofa

# Export Factor 1 and 3 weights for Cytoscape plot
data_weights_prot <- get_weights(object = mofa, views = "PROT", factors = c(1,3), as.data.frame = TRUE)

data_weights_prot <- as.data.frame(data_weights_prot)  %>%
  spread(factor, value) %>%
  mutate(Factor1_signflip = Factor1 *-1)

write.csv2(data_weights_prot, file = "output/data_weights_prot_fact1and3.csv")

## UMAP
plot.umap.data <-  plot_dimred(mofa, method="UMAP", color_by = "condition",stroke = 0.001, dot_size =1, alpha = 0.2, return_data = T)

plot.umap.all <- ggplot(plot.umap.data, aes(x=x, y = y, fill = color_by))+
  geom_point(size = 0.8, alpha = 0.6, shape = 21, stroke = 0) +
  theme_half_open() +
  scale_fill_manual(values = colorgradient6_manual, labels = c(labels), name = "Time aIg",)+
   theme(legend.position="none")+
  add.textsize +
  labs(title = "UMAP on MOFA+ factors shows minute and \nhour time-scale cell-state transitions", x = "UMAP 1", y = "UMAP 2")

## UMAP legend
legend.umap <- get_legend( ggplot(plot.umap.data, aes(x=x, y = y, fill = color_by))+
  geom_point(size = 2, alpha = 1, shape = 21, stroke = 0) +
  theme_half_open() +
  scale_fill_manual(values = colorgradient6_manual, labels = c(labels), name = "Time aIg \n(minutes)",)+
  add.textsize +
  labs(title = "UMAP on MOFA+ factors shows minute and \nhour time-scale cell-state transitions", x = "UMAP 1", "y = UMAP 2"))
legend.umap <- as_ggplot(legend.umap)

plot_fig1_umap <- plot_grid(plot.umap.all,legend.umap, labels = c(""), label_size = 10, ncol = 2, rel_widths = c(1, 0.2))

ggsave(plot_fig1_umap, filename = "output/paper_figures/Fig1_UMAP.pdf", width = 68, height = 62, units = "mm",  dpi = 300, useDingbats = FALSE)
ggsave(plot_fig1_umap, filename = "output/paper_figures/Fig1_UMAP.png", width = 68, height = 62, units = "mm",  dpi = 300)

plot_fig1_umap

seu_combined_selectsamples <- RunPCA(seu_combined_selectsamples,assay = "SCT.RNA", features = genes.variable, verbose = FALSE, ndims.print = 0, reduction.name = "pca.RNA")
seu_combined_selectsamples <- RunPCA(seu_combined_selectsamples, assay = "PROT", features = proteins.all, verbose = FALSE, ndims.print = 0,reduction.name = "pca.PROT")

## PCA analysis

plot.PCA.RNA <- DimPlot(seu_combined_selectsamples, reduction = "pca.RNA", group.by = "condition", pt.size =0.2) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg",)+
  labs(title = "RNA PCA separates cells \nat hour time-scale", x= "RNA PC 1", y = " RNA PC 2") +
  add.textsize +
  theme(legend.position = "none")


plot.PCA.PROT <- DimPlot(seu_combined_selectsamples, reduction = "pca.PROT", group.by = "condition", pt.size = 0.2) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg",)+
  labs(title = "(phospho-)protein PCA separates cells\nat minutes time-scale", x= "Protein PC 1", y = "Protein PC 2") +
  add.textsize +
  #  scale_x_reverse()+
  theme(legend.position = "none")

legend <- get_legend(DimPlot(seu_combined_selectsamples, reduction = "pca.RNA", group.by = "condition", pt.size =0.1) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg",)+
  add.textsize)
legend <- as_ggplot(legend)


PCA.PROTPC1.data <- data.frame(rank = 1:80, 
                               protein = names(sort(seu_combined_selectsamples@reductions$pca.PROT@feature.loadings[,1])),
                               weight.PC1 = sort(seu_combined_selectsamples@reductions$pca.PROT@feature.loadings[,1]),
                               highlights = c(names(sort(seu_combined_selectsamples@reductions$pca.PROT@feature.loadings[,1]))[1:4],rep("",76))
                               ) 

plot.PCA.PROTweightsPC1 <- ggplot(PCA.PROTPC1.data, aes(x=weight.PC1, y = rank, label = highlights)) +
  geom_point(size=0.1) +
  labs(title = "Top 4 protein loadings \n", x= "Weight Protein PC1") +
  geom_point(color = ifelse(PCA.PROTPC1.data$highlights == "", "grey50", "red")) +
  geom_text_repel(size = 2, segment.size = 0.25)+
  theme_half_open()+
  add.textsize +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()
        ) +
  scale_x_reverse()+
  scale_y_reverse()


PCA.RNAPC1.data <- data.frame(rank = 1:2371, 
                               RNA = names(sort(seu_combined_selectsamples@reductions$pca.RNA@feature.loadings[,1])),
                               weight.PC1 = sort(seu_combined_selectsamples@reductions$pca.RNA@feature.loadings[,1]),
                               highlights = c(rep("",2366),names(sort(seu_combined_selectsamples@reductions$pca.RNA@feature.loadings[,1]))[2367:2371])
                               ) 
## Loadings RNA
plot.PCA.RNAweightsPC1 <- ggplot(PCA.RNAPC1.data, aes(x=weight.PC1, y = rank, label = highlights)) +
  geom_point(size=0.1) +
  labs(title = "Top 4 RNA loadings \n", x= "Weight RNA PC1") +
  geom_point(color = ifelse(PCA.RNAPC1.data$highlights == "", "grey50", "red")) +
  geom_text_repel(size = 2, segment.size = 0.25)+
  theme_half_open()+
  add.textsize +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()
        ) 

## ridgeplots
PC1.data <- data.frame(sample = rownames(seu_combined_selectsamples@reductions$pca.PROT@cell.embeddings),
                       PC1_PROT = seu_combined_selectsamples@reductions$pca.PROT@cell.embeddings[,1],
                       PC1_RNA = seu_combined_selectsamples@reductions$pca.RNA@cell.embeddings[,1]) %>%
  left_join(meta.allcells)

plot_ridge_PC1Prot <- ggplot(PC1.data, aes(x = PC1_PROT, y = condition, fill = condition)) + 
  scale_fill_manual(values = colorgradient6_manual, labels = c(labels), name = "Time aIg \n(minutes)")+
  geom_density_ridges2() +
  scale_y_discrete(labels = labels,  name = "Time aIg (minutes)")+
  scale_x_continuous(name = "PC1 Proteins") +
  theme_half_open() +
  add.textsize+
  theme(legend.position = "none")

plot_ridge_PC1RNA <- ggplot(PC1.data, aes(x = PC1_RNA, y = condition, fill = condition)) + 
  scale_fill_manual(values = colorgradient6_manual, labels = c(labels), name = "Time aIg \n(minutes)")+
  geom_density_ridges2() +
  scale_y_discrete(labels = labels,  name = "Time aIg (minutes)")+
  scale_x_continuous(name = "PC1 RNA") +
  theme_half_open() +
  add.textsize+
  theme(legend.position = "none")

Fig1.pca <- plot_grid(plot.PCA.PROT,plot.PCA.RNA,legend, plot_ridge_PC1Prot, plot_ridge_PC1RNA,NULL,  plot.PCA.PROTweightsPC1, plot.PCA.RNAweightsPC1, labels = c(panellabels[1:2], "","","","",panellabels[3:4]), label_size = 10, ncol = 3, rel_widths = c(0.9,0.9,0.2,0.9,0.9), rel_heights = c(1,0.8,1.3))

ggsave(filename = "output/paper_figures/Suppl_PCA_aIg.pdf", width = 143, height = 170, units = "mm",  dpi = 300, useDingbats = FALSE)
ggsave(filename = "output/paper_figures/Suppl_PCA_aIg.png", width = 143, height = 170, units = "mm",  dpi = 300)

Fig1.pca

Fig.1.suppl.mofa.row1 <- plot_grid(plot.variance.perfactor.all, plot.variance.total,NULL, plot.heatmap.pval.covariates, labels = c(panellabels[1:3]), label_size = 10, ncol = 4, rel_widths = c(1.35, 0.27,0.2,0.38))

Fig.1.suppl.mofa.row2 <- plot_grid(plot.violin.factorall,legend.umap, labels = panellabels[4], label_size = 10, ncol = 2, rel_widths = c(1,0.1))

Suppl_mofa <- plot_grid(Fig.1.suppl.mofa.row1, Fig.1.suppl.mofa.row2, plot.rank.PROT.2.4to7, plot.rank.RNA.2.4to7, labels = c("","", panellabels[5:6]),label_size = 10, ncol = 1, rel_heights = c(1.45,1,1.1,1.1))


ggsave(Suppl_mofa, filename = "output/paper_figures/Suppl_MOFAaIg.pdf", width = 183, height = 220, units = "mm",  dpi = 300, useDingbats = FALSE)
ggsave(Suppl_mofa, filename = "output/paper_figures/Suppl_MOFAaIg.png", width = 183, height = 220, units = "mm",  dpi = 300)

Suppl_mofa

Fig2.row1.violinsprot <- plot_grid(`plot.violin.p-CD79a`, `plot.violin.p-Syk`, `plot.violin.p-JAK1`,labels = c(panellabels[3]), label_size = 10, ncol = 1, rel_heights = c(1.2,0.95,1.2))
Fig2.row1 <- plot_grid(plot.violin.factor1,plot.rank.PROT.1, plot.rank.RNA.1, Fig2.row1.violinsprot, labels = c(panellabels[1:2], ""),  label_size = 10, ncol =4, rel_widths = c(0.8,0.55,0.5,0.8))

Fig2.row2.violinsrna <- plot_grid(plot.violin.rna.NEAT1,plot.violin.rna.NPM1, plot.violin.rna.BTF3,labels = "", label_size = 10, ncol = 1, rel_heights = c(1.2,0.95,1.2))

Fig2.row2 <- plot_grid(plot.violin.factor3,plot.rank.PROT.3,plot.rank.RNA.3,Fig2.row2.violinsrna, labels = c(panellabels[4:5], "", panellabels[7]), label_size = 10, ncol =4, rel_widths = c(0.8,0.55,0.5,0.8))

Fig2.row3 <- plot_grid(plot.enriched.go.bp.factor3,legend.plot.enriched.go.bp.factor3,NULL, labels = panellabels[6], label_size = 10, ncol =3, rel_widths = c(1.8,0.2,0.6))

plot_fig2 <- plot_grid(Fig2.row1, Fig2.row2,Fig2.row3, labels = c("", "", ""), label_size = 10, ncol = 1, rel_heights = c(1,1,0.55))

ggsave(plot_fig2, filename = "output/paper_figures/Fig2.pdf", width = 183, height = 183, units = "mm",  dpi = 300, useDingbats = FALSE)
ggsave(plot_fig2, filename = "output/paper_figures/Fig2.png", width = 183, height = 183, units = "mm",  dpi = 300)

plot_fig2

ggsave(plot.enriched.go.bp.factor3_50, filename = "output/paper_figures/Suppl_Enrichm_Fct3.pdf", width = 183, height = 183, units = "mm",  dpi = 300, useDingbats = FALSE)
ggsave(plot.enriched.go.bp.factor3_50, filename = "output/paper_figures/Suppl_Enrichm_Fct3.png", width = 183, height = 183, units = "mm",  dpi = 300)
plot.enriched.go.bp.factor3_50