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

## Rename protein names
featurenamesmofa <- features_names(mofa)


## todo cleanup/more efficient
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="p-Src" ,"p-SRC",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="Syk" ,"SYK",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="p-c-Jun" ,"p-cJUN",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="Akt" ,"AKT",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="Btk" ,"BTK",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="p-Akt" ,"p-AKT",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="p-Syk" ,"p-SYK",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="Erk1/2" ,"ERK1/2",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="p-Erk1/2" ,"p-ERK1/2",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="p-Myc" ,"p-MYC",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="p-Btk" ,"p-BTK",x), how = "replace")
featurenamesmofa <- rapply(featurenamesmofa,function(x) ifelse(x=="Bcl-6" ,"BCL6",x), how = "replace")

features_names(mofa) <- featurenamesmofa

# 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")+
  scale_x_reverse()+
  scale_y_reverse()+
  add.textsize +
  labs(title = "UMAP on MOFA+ factors shows minute and \nhour timescale signal transductions", 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 timescale signal transductions", 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 timescale", 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 timescale", 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))
                               )

PCA.PROTPC1.data <- PCA.PROTPC1.data %>%
  mutate(highlights = as.character(highlights)) %>%
  mutate(highlights = case_when(as.character(highlights) == "p-Syk" ~ "p-SYK",
                                as.character(highlights) == "p-Src" ~ "p-SRC",
                           TRUE ~ highlights))

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:length(genes.variable), 
                               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("",(length(genes.variable)-4)),names(sort(seu_combined_selectsamples@reductions$pca.RNA@feature.loadings[,1]))[(length(genes.variable)-3):length(genes.variable)])
                               ) 
## 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")

seu_combined_selectsamples <- FindMultiModalNeighbors(
  seu_combined_selectsamples, reduction.list = list("pca.RNA", "pca.PROT"), 
  dims.list = list(c(1:7), c(1:7)), modality.weight.name = "RNA.weight"
)

seu_combined_selectsamples <- RunUMAP(seu_combined_selectsamples, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")

WNN.Dimplot <- DimPlot(seu_combined_selectsamples, reduction = 'wnn.umap', label = FALSE, repel = TRUE, label.size = 2.5) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg",)+
  labs(title = "Seurats WNN on PC's results in similar\nUMAP as MOFA+ based UMAP", x= "wnnUMAP 1", y = "wnnUMAP 2") +
  add.textsize +
  #  scale_x_reverse()+
  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,NULL,WNN.Dimplot,legend, labels = c(panellabels[1:2], "",panellabels[3:4],"",panellabels[5:6],"",panellabels[7]), 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,1))

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

Fig1.pca

## Variance per factor
plot.variance.perfactor.all <- plot_variance_explained(mofa, x="factor", y="view") +
    add.textsize +
    labs(title = "Variance explained by each factor per modality") 

## variance total
plot.variance.total <- plot_variance_explained(mofa, x="view", y="factor", plot_total = T)
plot.variance.total <- plot.variance.total[[2]] +
  add.textsize +
    labs(title = "Total \nvariance") +
  geom_text(aes(label=round(R2,1)), vjust=1.6, color="white", size=2.5)

## Significance correlation covariates
plot.heatmap.pval.covariates <- as.ggplot(correlate_factors_with_covariates(mofa, 
  covariates = c("time"),
  factors = 1:mofa@dimensions$K,
  fontsize = 7, 
  cluster_row = F,
  cluster_col = F
))+ 
  add.textsize +
  theme(axis.text.y=element_blank(),
        axis.text.x=element_blank(),
        plot.title = element_text(size=7, face = "plain"),
        )

## Factor values over time
plot.violin.factorall <- plot_factor(mofa, 
  factor = c(2,4:mofa@dimensions$K),
  color_by = "condition",
  dot_size = 0.2,      # change dot size
  dodge = T,           # dodge points with different colors
  legend = F,          # remove legend
  add_violin = T,      # add violin plots,
  violin_alpha = 0.9  # transparency of violin plots
) +   
  add.textsize +
  scale_color_manual(values=c(colorgradient6_manual2, labels = c(labels), name = "Time aIg")) +
  scale_fill_manual(values=c(colorgradient6_manual2, labels = c(labels), name = "Time aIg"))+
  labs(title = "Factor values per time-point of additional factors not correlating with time." ) +
  theme(axis.text.x=element_blank())

## Loadings factors stable protein

plot.rank.PROT.2.4to7 <- plot_weights(mofa, 
  view = "PROT", 
  factors = c(1:mofa@dimensions$K), 
  nfeatures = 3, 
  text_size = 1.5
) +   
  add.textsize +
  labs(title = "Top 3 Protein loadings per factor") +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank()
)

## Loadings factors stable protein

plot.rank.RNA.2.4to7 <- plot_weights(mofa, 
  view = "RNA", 
  factors = c(1:mofa@dimensions$K), 
  nfeatures = 3, 
  text_size = 1.5
) +   
  add.textsize +
  labs(title = "Top 3 RNA loadings per factor") +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank()
)

## correlation time and factors
plot.correlation.covariates <- correlate_factors_with_covariates(mofa, 
  covariates = c("time"),
  factors = mofa@dimensions$K:1,
  plot = "r"
)

plot.correlation.covariates <- ggcorrplot(plot.correlation.covariates, tl.col = "black", method = "square", lab = TRUE, ggtheme = theme_void, colors = c("#11304C", "white", "#11304C"), lab_size = 2.5) +
  add.textsize +
  labs(title = "Correlation of \nfactors with \ntime of treatment", y = "") +
  scale_y_discrete(labels = "") +
  coord_flip() +
  theme(axis.text.x=element_text(angle =0,hjust = 0.5), 
        axis.text.y=element_text(size = 5),
        legend.position="none",
        plot.title = element_text(hjust = 0.5))

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.3,0.25,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/Fig2.Suppl_MOFAaIg.pdf", width = 183, height = 220, units = "mm",  dpi = 300, useDingbats = FALSE)
ggsave(Suppl_mofa, filename = "output/paper_figures/Fig2.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[4]), 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, NULL, Fig2.row1.violinsprot, labels = c(panellabels[1:3], ""),  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[5:6], "", panellabels[8]), 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[7], 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

proteindata_counts <- as.data.frame(t(seu_combined_selectsamples@assays$PROT@scale.data)) %>%
  mutate(sample = rownames(t(seu_combined_selectsamples@assays$PROT@scale.data)))


metadata.all <- as.data.frame(seu_combined_selectsamples@meta.data) %>%
  mutate(sample = rownames((seu_combined_selectsamples@meta.data)))

proteindata_counts <- left_join(proteindata_counts, metadata.all)



toppJAK1 <- proteindata_counts %>%
  group_by(time) %>%
  top_frac(wt = `p-JAK1`, n = 0.05)

bottompJAK1 <- proteindata_counts %>%
  group_by(time) %>%
  top_frac(wt = `p-JAK1`, n = -0.05)

proteindata_counts <- proteindata_counts %>%
  mutate(highlowpJAK1 = ifelse(sample %in% toppJAK1$sample, "p-JAK1 high", ifelse(sample %in% bottompJAK1$sample, "p-JAK1 low","middle")))

addmeta <- proteindata_counts[,c("highlowpJAK1", "sample")]
rownames(addmeta) <- proteindata_counts$sample

seu_combined_selectsamples <- AddMetaData(seu_combined_selectsamples, addmeta)

seu.JAK1.180 <- subset(seu_combined_selectsamples, condition == "180.aIg.contr" & highlowpJAK1 != "middle")

seu.JAK1.180 <- SetIdent(seu.JAK1.180, value = "highlowpJAK1")

# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
markers.180 <- FindMarkers(seu.JAK1.180,ident.1 = "p-JAK1 high", ident.2 = "p-JAK1 low", assay = "PROT", slot = "scale.data", logfc.threshold = 0, return.thresh = 1, only.pos = F)
# view results
#markers.180 <- filter(markers.180, cluster == "p-JAK1 high")
#markers.180


# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
markers.180.RNA <-  FindAllMarkers(seu.JAK1.180,assay = "RNA", slot = "data", logfc.threshold = 0.3, return.thresh = 0.01, only.pos = T,min.pct = 0.1)
markers.180.RNA <-   FindMarkers(seu.JAK1.180,ident.1 = "p-JAK1 high", ident.2 = "p-JAK1 low", assay = "SCT.RNA", slot = "scale.data", logfc.threshold = 0, return.thresh = 1, only.pos = F)
# view results
#markers.180.RNA

markers.180$protein <-rownames(markers.180)
# add a column of NAs
markers.180$diffexpressed <- "NO"
# if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
markers.180$diffexpressed[markers.180$avg_diff > 0.25 & markers.180$p_val_adj < 0.01] <- "UP"
# if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
markers.180$diffexpressed[markers.180$avg_diff  < -0.25 & markers.180$p_val_adj < 0.01] <- "DOWN"

mycolors <- c("blue", "red", "black")
names(mycolors) <- c("DOWN", "UP", "NO")


markers.180$delabel <- NA
markers.180$delabel[markers.180$diffexpressed != "NO"] <- markers.180$protein[markers.180$diffexpressed != "NO"]

# Finally, we can organize the labels nicely using the "ggrepel" package and the geom_text_repel() function
# load library

# plot adding up all layers we have seen so far
plot.vulcano.180min <- ggplot(data=markers.180, aes(x=avg_diff , y=-log10(p_val_adj), col=diffexpressed, label=delabel)) +
        geom_point(size=0.5) + 
        theme_minimal() +
        geom_text_repel(size=2.2) +
        scale_color_manual(values=c("blue", "red", "black")) +
        geom_vline(xintercept=c(-0.25, 0.25), col="red") +
        geom_hline(yintercept=-log10(0.01), col="red") +
  labs(x = expression("Log"[2]*" Fold Change"), y = expression("-log"[10]*" adjusted p-value"), title = "p-JAK1 high vs p-JAK1 low (t = 180 min)") &
  add.textsize

sign.markers180 <- markers.180$protein[markers.180$avg_diff > 0.25 & markers.180$p_val_adj < 0.01 | markers.180$avg_diff  < -0.25 & markers.180$p_val_adj < 0.01] 

plot.vln.180min <- VlnPlot(seu.JAK1.180,assay = "PROT",slot = "scale.data", features =  sign.markers180, group.by = "highlowpJAK1",ncol = 6, pt.size = 0.5) &
  add.textsize

## RNA
markers.180.RNA$protein <-rownames(markers.180.RNA)
# add a column of NAs
markers.180.RNA$diffexpressed <- "NO"
# if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
markers.180.RNA$diffexpressed[markers.180.RNA$avg_diff > 0.25 & markers.180.RNA$p_val < 0.05] <- "UP"
# if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
markers.180.RNA$diffexpressed[markers.180.RNA$avg_diff < -0.25 & markers.180.RNA$p_val < 0.05] <- "DOWN"

mycolors <- c("blue", "red", "black")
names(mycolors) <- c("DOWN", "UP", "NO")


markers.180.RNA$delabel <- NA
markers.180.RNA$delabel[markers.180.RNA$diffexpressed != "NO"] <- markers.180.RNA$protein[markers.180.RNA$diffexpressed != "NO"]

# plot adding up all layers we have seen so far
plot.vulcano.180min.RNA <- ggplot(data=markers.180.RNA, aes(x=avg_diff, y=-log10(p_val), col=diffexpressed, label=delabel)) +
        geom_point() + 
        theme_minimal() +
        geom_text_repel(size=2.2) +
        scale_color_manual(values=c("blue", "red", "black")) +
        geom_vline(xintercept=c(-0.25, 0.25), col="red") +
        geom_hline(yintercept=-log10(0.05), col="red") +
  labs(x = expression("Log"[2]*" Fold Change"), y = expression("-log"[10]*" p-value"), title = "p-JAK1 high vs p-JAK1 low (t = 180 min)") &
  add.textsize

sign.markers180.RNA <- markers.180.RNA$protein[markers.180.RNA$avg_diff > 0.25 & markers.180.RNA$p_val < 0.05] 

plot.vln.180min.RNA <- VlnPlot(seu.JAK1.180,assay = "RNA", features =  sign.markers180.RNA[1:20], group.by = "highlowpJAK1",ncol = 10) &
  add.textsize

plot_180min <- plot_grid(plot.vulcano.180min, plot.vln.180min, labels = panellabels[c(5,6)], label_size = 10, ncol = 2, rel_widths = c(1,2))

## 6 minutes
seu.JAK1.006 <- subset(seu_combined_selectsamples, condition == "006.aIg.contr" & highlowpJAK1 != "middle")

seu.JAK1.006 <- SetIdent(seu.JAK1.006, value = "highlowpJAK1")

# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
markers.006 <- FindMarkers(seu.JAK1.006,ident.1 = "p-JAK1 high", ident.2 = "p-JAK1 low", assay = "PROT", slot = "scale.data", logfc.threshold = 0, return.thresh = 1, only.pos = F)
# view results
#markers.006 <- filter(markers.006, cluster == "p-JAK1 high")
#markers.006


# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
markers.006.RNA <-  FindAllMarkers(seu.JAK1.006,assay = "RNA", slot = "data", logfc.threshold = 0.3, return.thresh = 0.01, only.pos = T,min.pct = 0.1)
markers.006.RNA <-   FindMarkers(seu.JAK1.006,ident.1 = "p-JAK1 high", ident.2 = "p-JAK1 low", assay = "SCT.RNA", slot = "scale.data", logfc.threshold = 0, return.thresh = 1, only.pos = F)
# view results
#markers.006.RNA

markers.006$protein <-rownames(markers.006)
# add a column of NAs
markers.006$diffexpressed <- "NO"
# if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
markers.006$diffexpressed[markers.006$avg_diff > 0.25 & markers.006$p_val_adj < 0.01] <- "UP"
# if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
markers.006$diffexpressed[markers.006$avg_diff  < -0.25 & markers.006$p_val_adj < 0.01] <- "DOWN"

mycolors <- c("blue", "red", "black")
names(mycolors) <- c("DOWN", "UP", "NO")


markers.006$delabel <- NA
markers.006$delabel[markers.006$diffexpressed != "NO"] <- markers.006$protein[markers.006$diffexpressed != "NO"]

# Finally, we can organize the labels nicely using the "ggrepel" package and the geom_text_repel() function
# load library
# plot adding up all layers we have seen so far
plot.vulcano.006min <- ggplot(data=markers.006, aes(x=avg_diff , y=-log10(p_val_adj), col=diffexpressed, label=delabel)) +
        geom_point(size=0.5) + 
        theme_minimal() +
        geom_text_repel(size=2.2) +
        scale_color_manual(values=c("blue", "red", "black")) +
        geom_vline(xintercept=c(-0.25, 0.25), col="red") +
        geom_hline(yintercept=-log10(0.01), col="red") +
  labs(x = expression("Log"[2]*" Fold Change"), y = expression("-log"[10]*" adjusted p-value"), title = "p-JAK1 high vs p-JAK1 low (t = 006 min)") &
  add.textsize

sign.markers006 <- markers.006$protein[markers.006$avg_diff > 0.25 & markers.006$p_val_adj < 0.01 | markers.006$avg_diff  < -0.25 & markers.006$p_val_adj < 0.01] 

plot.vln.006min <- VlnPlot(seu.JAK1.006,assay = "PROT",slot = "scale.data", features =  sign.markers006, group.by = "highlowpJAK1",ncol = 6, pt.size = 0.5) &
  add.textsize

markers.006.RNA$protein <-rownames(markers.006.RNA)
# add a column of NAs
markers.006.RNA$diffexpressed <- "NO"
# if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
markers.006.RNA$diffexpressed[markers.006.RNA$avg_diff > 0.25 & markers.006.RNA$p_val < 0.05] <- "UP"
# if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
markers.006.RNA$diffexpressed[markers.006.RNA$avg_diff < -0.25 & markers.006.RNA$p_val < 0.05] <- "DOWN"

mycolors <- c("blue", "red", "black")
names(mycolors) <- c("DOWN", "UP", "NO")


markers.006.RNA$delabel <- NA
markers.006.RNA$delabel[markers.006.RNA$diffexpressed != "NO"] <- markers.006.RNA$protein[markers.006.RNA$diffexpressed != "NO"]

# Finally, we can organize the labels nicely using the "ggrepel" package and the geom_text_repel() function
# load library
# plot adding up all layers we have seen so far
plot.vulcano.006min.RNA <- ggplot(data=markers.006.RNA, aes(x=avg_diff, y=-log10(p_val), col=diffexpressed, label=delabel)) +
        geom_point() + 
        theme_minimal() +
        geom_text_repel(size=2.2) +
        scale_color_manual(values=c("blue", "red", "black")) +
        geom_vline(xintercept=c(-0.25, 0.25), col="red") +
        geom_hline(yintercept=-log10(0.05), col="red") +
  labs(x = expression("Log"[2]*" Fold Change"), y = expression("-log"[10]*" p-value"), title = "p-JAK1 high vs p-JAK1 low (t = 006 min)") &
  add.textsize

sign.markers006.RNA <- markers.006.RNA$protein[markers.006.RNA$avg_diff > 0.25 & markers.006.RNA$p_val < 0.05] 

plot.vln.006min.RNA <- VlnPlot(seu.JAK1.006,assay = "RNA", features =  sign.markers006.RNA[1:20], group.by = "highlowpJAK1",ncol = 10) &
  add.textsize

plot_006min <- plot_grid(plot.vulcano.006min, plot.vln.006min, labels = panellabels[c(3,4)], label_size = 10, ncol = 2, rel_widths = c(1,2))

## 2 minutes
seu.JAK1.002 <- subset(seu_combined_selectsamples, condition == "002.aIg.contr" & highlowpJAK1 != "middle")

seu.JAK1.002 <- SetIdent(seu.JAK1.002, value = "highlowpJAK1")

# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
markers.002 <- FindMarkers(seu.JAK1.002,ident.1 = "p-JAK1 high", ident.2 = "p-JAK1 low", assay = "PROT", slot = "scale.data", logfc.threshold = 0, return.thresh = 1, only.pos = F)
# view results
#markers.002 <- filter(markers.002, cluster == "p-JAK1 high")
#markers.002


# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
markers.002.RNA <-  FindAllMarkers(seu.JAK1.002,assay = "RNA", slot = "data", logfc.threshold = 0.3, return.thresh = 0.01, only.pos = T,min.pct = 0.1)
markers.002.RNA <-   FindMarkers(seu.JAK1.002,ident.1 = "p-JAK1 high", ident.2 = "p-JAK1 low", assay = "SCT.RNA", slot = "scale.data", logfc.threshold = 0, return.thresh = 1, only.pos = F)
# view results
#markers.002.RNA

markers.002$protein <-rownames(markers.002)
# add a column of NAs
markers.002$diffexpressed <- "NO"
# if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
markers.002$diffexpressed[markers.002$avg_diff > 0.25 & markers.002$p_val_adj < 0.01] <- "UP"
# if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
markers.002$diffexpressed[markers.002$avg_diff  < -0.25 & markers.002$p_val_adj < 0.01] <- "DOWN"

mycolors <- c("blue", "red", "black")
names(mycolors) <- c("DOWN", "UP", "NO")


markers.002$delabel <- NA
markers.002$delabel[markers.002$diffexpressed != "NO"] <- markers.002$protein[markers.002$diffexpressed != "NO"]

plot.vulcano.002min <- ggplot(data=markers.002, aes(x=avg_diff , y=-log10(p_val_adj), col=diffexpressed, label=delabel)) +
        geom_point(size=0.5) + 
        theme_minimal() +
        geom_text_repel(size=2.2) +
        scale_color_manual(values=c("blue", "red", "black")) +
        geom_vline(xintercept=c(-0.25, 0.25), col="red") +
        geom_hline(yintercept=-log10(0.01), col="red") +
  labs(x = expression("Log"[2]*" Fold Change"), y = expression("-log"[10]*" adjusted p-value"), title = "p-JAK1 high vs p-JAK1 low (t = 002 min)") &
  add.textsize

sign.markers002 <- markers.002$protein[markers.002$avg_diff > 0.25 & markers.002$p_val_adj < 0.01 | markers.002$avg_diff  < -0.25 & markers.002$p_val_adj < 0.01] 

plot.vln.002min <- VlnPlot(seu.JAK1.002,assay = "PROT",slot = "scale.data", features =  sign.markers002, group.by = "highlowpJAK1",ncol = 6, pt.size = 0.5) &
  add.textsize

markers.002.RNA$protein <-rownames(markers.002.RNA)
# add a column of NAs
markers.002.RNA$diffexpressed <- "NO"
# if log2Foldchange > 0.6 and pvalue < 0.05, set as "UP" 
markers.002.RNA$diffexpressed[markers.002.RNA$avg_diff > 0.25 & markers.002.RNA$p_val < 0.05] <- "UP"
# if log2Foldchange < -0.6 and pvalue < 0.05, set as "DOWN"
markers.002.RNA$diffexpressed[markers.002.RNA$avg_diff < -0.25 & markers.002.RNA$p_val < 0.05] <- "DOWN"

mycolors <- c("blue", "red", "black")
names(mycolors) <- c("DOWN", "UP", "NO")


markers.002.RNA$delabel <- NA
markers.002.RNA$delabel[markers.002.RNA$diffexpressed != "NO"] <- markers.002.RNA$protein[markers.002.RNA$diffexpressed != "NO"]

plot.vulcano.002min.RNA <- ggplot(data=markers.002.RNA, aes(x=avg_diff, y=-log10(p_val), col=diffexpressed, label=delabel)) +
        geom_point() + 
        theme_minimal() +
        geom_text_repel(size=2.2) +
        scale_color_manual(values=c("blue", "red", "black")) +
        geom_vline(xintercept=c(-0.25, 0.25), col="red") +
        geom_hline(yintercept=-log10(0.05), col="red") +
  labs(x = expression("Log"[2]*" Fold Change"), y = expression("-log"[10]*" p-value"), title = "p-JAK1 high vs p-JAK1 low (t = 002 min)") &
  add.textsize

sign.markers002.RNA <- markers.002.RNA$protein[markers.002.RNA$avg_diff > 0.25 & markers.002.RNA$p_val < 0.05] 

plot.vln.002min.RNA <- VlnPlot(seu.JAK1.002,assay = "RNA", features =  sign.markers002.RNA[1:20], group.by = "highlowpJAK1",ncol = 10) &
  add.textsize

plot_002min <- plot_grid(plot.vulcano.002min, plot.vln.002min, labels = panellabels[c(1,2)], label_size = 10, ncol = 2, rel_widths = c(1,2))

plot_all <- plot_grid(plot_002min,plot_006min,plot_180min , labels =c("", "", ""), ncol = 1, rel_heights = c(1,1,1))

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

plot_all