MOFA model on aIg stimulated BJAB cells

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

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

MOFAfactors<- as.data.frame(mofa@expectations$Z) %>%
  mutate(sample = rownames(as.data.frame(mofa@expectations$Z)[,1:mofa@dimensions$K])) 

MOFAfactors <- left_join(as.data.frame(MOFAfactors), meta.allcells)

weights.prot <- get_weights(mofa, views = "PROT",as.data.frame = TRUE)

topnegprots.factor1 <- weights.prot %>%
  filter(factor == "Factor1" & value <= 0) %>%
  arrange(value)

topposprots.factor3 <- weights.prot %>%
  filter(factor == "Factor3" & value >= 0) %>%
  arrange(-value)


weights.RNA <- get_weights(mofa, views = "RNA",as.data.frame = TRUE)

## violin plots phospho-proteins highlighted in main

f.violin.fact <- function(data , protein){
  
  ggplot(subset(data)  , aes(x = as.factor(condition), y =get(noquote(protein)))) +
    annotate("rect",
          xmin = 5 - 0.45,
             xmax = 6 + 0.6,
           ymin = -4.5, ymax = 4.5, fill = "lightblue",
           alpha = .4
  )+
    geom_hline(yintercept=0, linetype='dotted', col = 'black')+ 
  geom_violin(alpha = 0.9,aes(fill = condition))+
  geom_jitter(width = 0.05,size = 0.1, fill = "black", shape = 21)+ 
  stat_summary(fun=median, geom="point", shape=95, size=2, inherit.aes = T, position = position_dodge(width = 0.9), color = "red")+
  theme_few()+
  ylab(paste0(protein)) +
  scale_x_discrete(labels = labels, expand = c(0.1,0.1), name = "Time aIg (minutes)") +
    scale_y_continuous(expand = c(0,0), name = strsplit(protein, split = "\\.")[[1]][2]) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg ",)+ 
  scale_fill_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg ",) +
  add.textsize +
  theme(#axis.title.x=element_blank(),
        #axis.text.x=element_blank(),
        #axis.ticks.x=element_blank(),
        legend.position="none") 
}

factors_toplot <- c(colnames(MOFAfactors)[c(1,3)])

for(i in factors_toplot) {
assign(paste0("plot.factor.", i), f.violin.fact(data = MOFAfactors,protein = i)) 
}

legend.violinfactor <- as_ggplot( get_legend( ggplot(subset(MOFAfactors)  , aes(x = as.factor(condition), y =get(noquote(factors_toplot[1])))) +
    annotate("rect",
          xmin = 5 - 0.45,
             xmax = 6 + 0.6,
           ymin = -4.5, ymax = 4.5, fill = "lightblue",
           alpha = .4
  )+
    geom_hline(yintercept=0, linetype='dotted', col = 'black')+ 
  geom_violin(alpha = 0.9,aes(fill = condition))+
  geom_jitter(width = 0.05,size = 0.1, aes(col = condition), fill = "black", shape = 21)+ 
  stat_summary(fun=median, geom="point", shape=23, size=2, inherit.aes = T, position = position_dodge(width = 0.9), color = "black")+
  theme_few()+
  ylab(paste0(factors_toplot[1])) +
  scale_x_discrete(labels = labels, expand = c(0.1,0.1), name = "Time aIg  (minutes)") +
    scale_y_continuous(expand = c(0,0), name = strsplit(factors_toplot[1], split = "\\.")[[1]][2]) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg \n(minutes)",)+ 
  scale_fill_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg \n(minutes)",) +
  add.textsize ) ) 

plot.violin.factor1 <- plot.factor.group1.Factor1 +
  scale_y_reverse(expand = c(0,0), name = "Factor 1 value")+
  annotate(geom="text", x=1.5, y=-4.3, label="minutes", size = 2,
              color="grey3") +
  annotate(geom="text", x=5.5, y=-4.3, label="hours", size = 2,
              color="grey3") +
  labs(title = "Factor 1 captures \nminute times-scale signal transduction")

plot.violin.factor3 <- plot.factor.group1.Factor3 +
  annotate(geom="text", x=5.5, y=4.3, label="hours", size = 2,
              color="grey3")+
  annotate(geom="text", x=1.5, y=4.3, label="minutes", size = 2,
              color="grey3") +
  scale_y_continuous(expand = c(0,0), name = "Factor 3 value")+
  labs(title = "Factor 3 captures \nhour timescale signal transduction")

#### protein Loadings factor 1

list.bcellpathway.protein <- c("p-CD79a", "p-SYK", "p-SRC", "p-ERK1/2", "p-PLC-y2", "p-BLNK","p-PLC-y2Y759","p-PKC-b1", "p-p38", "p-AKT", "p-S6", "p-TOR", "CD79a") # , "p-PLC-y2", "p-PKC-b1", "p-IKKa/b","p-JNK", "p-p38", "p-p65","p-Akt", "p-S6", "p-TOR"
list.top20 <- topnegprots.factor1$feature[1:20]
## protein factor 1
plotdata.rank.PROT.1 <-plot_weights(mofa, 
  view = "PROT", 
  factors = c(1), 
  nfeatures = 15, 
  text_size = 1,
  manual = list(list.top20, list.bcellpathway.protein, "p-JAK1" ), 
  color_manual = list("grey","blue4","red"),
  return_data = TRUE
)

plotdata.rank.PROT.1 <- plotdata.rank.PROT.1 %>%
  mutate(Rank = 1:80,
         Weight = value, 
         colorvalue = ifelse(labelling_group == 2 & value <= -0.2,"grey", ifelse(labelling_group == 3 & value <= -0.2, "blue4", ifelse(labelling_group == 4 & value <= -0.2, "red", "grey"))),
         highlights = ifelse(labelling_group >= 1 & value <= -0.25, as.character(feature), "")
         ) 

plot.rank.PROT.1 <- ggplot(plotdata.rank.PROT.1, aes(x=Weight, y = Rank, label = highlights)) +
  labs(title = "(Phospho)proteins: <p><span style='color:blue4'>BCR </span>&<span style='color:red'> JAK1</span> signaling", #<span style='color:blue4'>BCR signaling</span> and <span style='color:red'>p-JAK1</span>
       x= "Factor 1 loading value",
       y= "Factor 1 loading rank") +
  geom_point(size=0.1, color =plotdata.rank.PROT.1$colorvalue, size =1) +
  geom_text_repel(size = 2, 
                  segment.size = 0.2, 
                  color =plotdata.rank.PROT.1$colorvalue,
                  nudge_x       = 1 - plotdata.rank.PROT.1$Weight,
                  direction     = "y",
                  hjust         = 1,
                  segment.color = "grey50")+
  theme_few()+
  add.textsize  +
  scale_y_continuous(trans = "reverse") +
  add.textsize +
  theme(
    plot.title = element_markdown()
  ) +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank()
        )+
  xlim(c(-1,1))


## violin plots phospho-proteins highlighted in main

f.violin.prot <- function(data , protein){
  
  ggplot(subset(data)  , aes(x = as.factor(condition), y =get(noquote(protein)))) +
    annotate("rect",
          xmin = 5 - 0.45,
             xmax = 6 + 0.6,
           ymin = -5.5, ymax = 5.5, fill = "lightblue",
           alpha = .4
  )+
  geom_violin(alpha = 0.9,aes(fill = condition))+
  geom_jitter(width = 0.05,size = 0.1, color = "black")+ 
  stat_summary(fun=median, geom="point", shape=95, size=2, inherit.aes = T, position = position_dodge(width = 0.9), color = "red")+
  theme_few()+
  ylab(paste0(protein)) +
  scale_x_discrete(labels = labels, expand = c(0.1,0.1), name = "Time aIg  (minutes)") +
  scale_y_continuous(expand = c(0,0), name = protein) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg ",)+ 
  scale_fill_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg ",) +
  add.textsize +
  theme(#axis.title.x=element_blank(),
        #axis.text.x=element_blank(),
        #axis.ticks.x=element_blank(),
        legend.position="none") 
}

proteins_toplot <- c("p-CD79a","p-Syk","p-JAK1", "IgM")

for(i in proteins_toplot) {
assign(paste0("plot.violin.", i), f.violin.prot(data = proteindata ,protein = i)) 

}

## Fig 2 bottom row panels
`plot.violin.p-CD79a` <- `plot.violin.p-CD79a` +
  labs(title = "BCR signaling pathway and \nJAK1 activation")+
  add.textsize+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),)

`plot.violin.p-Syk` <- `plot.violin.p-Syk` +
  add.textsize+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),) +
  scale_y_continuous(name = "p-SYK")

`plot.violin.p-JAK1` <-`plot.violin.p-JAK1` +
  add.textsize +
  labs(y = "<span style='color:red'>p-JAK1</span>") +
  theme(axis.title.y = element_markdown()) 

topposrna.factor3 <- weights.RNA %>%
  filter(factor == "Factor3" & value >= 0) %>%
  arrange(-value)
rownames(topposrna.factor3) <- topposrna.factor3$feature

### Convert Gene-names to gene IDs (using 'org.Hs.eg.db' library)

topposrna.factor3 <- topposrna.factor3 %>%
  mutate(ENTREZID = mapIds(org.Hs.eg.db, as.character(topposrna.factor3$feature), 'ENTREZID', 'SYMBOL'))

listgenes.factor3 <- topposrna.factor3$value

names(listgenes.factor3) <- topposrna.factor3$ENTREZID


go.pb.fct3.pos <- enrichGO(gene         = topposrna.factor3$feature,
                OrgDb         = org.Hs.eg.db,
                keyType       = 'SYMBOL',
                ont           = "BP",
                pAdjustMethod = "BH")

go.pb.fct3.pos <- simplify(go.pb.fct3.pos, cutoff=0.6, by="p.adjust", select_fun=min)


go.pb.fct3.pos.dplyr <- mutate(go.pb.fct3.pos, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio))) 

## plot main panel
plot.enriched.go.bp.factor3 <- ggplot(go.pb.fct3.pos.dplyr, showCategory = 10, 
  aes(-log10(p.adjust), fct_reorder(Description,  -log10(p.adjust)))) +   geom_segment(aes(xend=0, yend = Description)) +
  geom_point(aes(size = Count)) +
  scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
  scale_size_continuous(range=c(0.5, 3)) +
  theme_minimal() + 
  add.textsize +
  xlab("-log adj pval") +
  ylab(NULL) + 
  ggtitle("Enriched Biological Processes") +
  theme(legend.position="none")

legend.plot.enriched.go.bp.factor3 <- as.ggplot(get_legend(ggplot(go.pb.fct3.pos.dplyr, showCategory = 10, 
  aes(-log10(p.adjust), fct_reorder(Description,  -log10(p.adjust)))) +   geom_segment(aes(xend=0, yend = Description)) +
  geom_point(aes(size = Count)) +
  scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
  scale_size_continuous(range=c(0.5, 3)) +
  theme_minimal() + 
  add.textsize +
  xlab("-log adj pval") +
  ylab(NULL) + 
  ggtitle("Enriched Biological Processes") +
     theme(legend.position="right")))

## Plot supplementary
plot.enriched.go.bp.factor3_50 <- ggplot(go.pb.fct3.pos.dplyr, showCategory = 50, 
  aes(-log10(p.adjust), fct_reorder(Description,  -log10(p.adjust)))) +   geom_segment(aes(xend=0, yend = Description)) +
  geom_point(aes(size = Count)) +
  scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
  scale_size_continuous(range=c(0.5, 3)) +
  theme_minimal() + 
  add.textsize +
  xlab("-log adj pval") +
  ylab(NULL) + 
  ggtitle("Enriched Biological Processes") 

message(c("Significance protein folding GO term: \n",go.pb.fct3.pos.dplyr@result[which(go.pb.fct3.pos.dplyr@result$Description == "protein folding"),"p.adjust"]))

##For RNA loadings panels
topgeneset.fct3<- unlist(str_split(go.pb.fct3.pos.dplyr@result[1:10,"geneID"], pattern = "/"))

topgeneset.fct3 = bitr(topgeneset.fct3, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")

## RNA factor 1 loadings
plotdata.rank.RNA.3 <-plot_weights(mofa, 
  view = "RNA", 
  factors = c(3), 
  nfeatures = 5, 
  text_size = 1,
  manual = list(topposrna.factor3$feature[c(1:5)] ,topgeneset.fct3$SYMBOL), 
  color_manual = list("grey","blue4"),
  return_data = TRUE
)

plotdata.rank.RNA.3 <- plotdata.rank.RNA.3 %>%
  mutate(Rank = 1:nrow(plotdata.rank.RNA.3),
         Weight = value, 
         colorvalue = ifelse(labelling_group == 3 &value >= 0.25,"blue4", ifelse(labelling_group == 2&value >= 0.25, "grey2", "grey")),
         highlights = ifelse(labelling_group >= 1&value >= 0.25, as.character(feature), "")
         )

plot.rank.RNA.3 <- ggplot(plotdata.rank.RNA.3, aes(x=Weight, y = Rank, label = highlights)) +
  labs(title = "Contributing genes <p><span style='color:blue4'><span style='color:blue4'>GO-term vesicle related </span> ",  #<span style='color:blue4'>BCR signaling regulation (find GO term+list todo)</span>
       x= "Factor 3 loading value",
       y= "Factor 3 loading rank") +
  geom_point(size=0.1, color =plotdata.rank.RNA.3$colorvalue) +
  geom_text_repel(size = 2, 
                  segment.size = 0.2, 
                  color =plotdata.rank.RNA.3$colorvalue,
                  nudge_x       = -1 - plotdata.rank.RNA.3$Weight,
                  direction     = "y",
                  hjust         = 0,
                  segment.color = "grey50")+
  theme_few()+
  add.textsize +
  scale_y_continuous() +
  add.textsize +
  theme(
    plot.title = element_markdown()
  )  +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.title.y = element_blank())+
  xlim(c(-1,1))


#### Protein loadings factor 3

list.highlight.prot.fct3 <- c("p-p38", "p-AKT", "p-S6", "p-TOR", "XBP1_PROT", "p-STAT5")
list.top20 <- topnegprots.factor1$feature[1:20]

plotdata.rank.PROT.3 <-plot_weights(mofa, 
  view = "PROT", 
  factors = c(3), 
  nfeatures = 30, 
  text_size = 1,
  manual = list(topposprots.factor3$feature[c(1:8)] ,list.highlight.prot.fct3), 
  color_manual = list("grey","blue4"),
  return_data = TRUE
)

plotdata.rank.PROT.3 <- plotdata.rank.PROT.3 %>%
  mutate(Rank = 1:80,
         Weight = value, 
         colorvalue = ifelse(labelling_group == 3 &value >= 0.25,"blue4", ifelse(labelling_group == 2&value >= 0.4, "grey", "grey")),
         highlights = ifelse(labelling_group >= 1&value >= 0.25, as.character(feature), "")
         )  %>%
  mutate(highlights = case_when(as.character(highlights) == "XBP1_PROT" ~ "XBP1",
                           TRUE ~ highlights))

plot.rank.PROT.3 <- ggplot(plotdata.rank.PROT.3, aes(x=Weight, y = Rank, label = highlights)) +
  labs(title = "Signaling implicated in <p><span style='color:blue4'>Unfolded protein response</span>", 
       x= "Factor 3 loading value",
       y= "Factor 3 loading rank") +
  geom_point(size=0.1, color =plotdata.rank.PROT.3$colorvalue, size =1) +
  geom_text_repel(size = 2, 
                  segment.size = 0.2, 
                  color =plotdata.rank.PROT.3$colorvalue,
                  nudge_x       = -1 - plotdata.rank.PROT.3$Weight,
                  direction     = "y",
                  hjust         = 0,
                  segment.color = "grey50")+
  theme_few()+
  add.textsize  +
  scale_y_continuous() +
  add.textsize +
  theme(
    plot.title = element_markdown()
  ) +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank()
        )+
  xlim(c(-1,1))


#### RNA loadings factor 1

plotdata.rank.RNA.1 <-plot_weights(mofa, 
  view = "RNA", 
  factors = c(1), 
  nfeatures = 2, 
  text_size = 1,
  return_data = TRUE
)

plotdata.rank.RNA.1 <- plotdata.rank.RNA.1 %>%
  mutate(Rank = 1:nrow(plotdata.rank.RNA.1),
         Weight = value, 
         colorvalue = ifelse(labelling_group == 1,"blue4", ifelse(labelling_group == 1, "blue4", "grey")),
         highlights = ifelse(labelling_group >= 1, as.character(feature), "")
         )

plot.rank.RNA.1 <- ggplot(plotdata.rank.RNA.1, aes(x=Weight, y = Rank, label = highlights)) +
  labs(title = " Contributing genes<p><span style='color:blue4'>BCR activation</span>",  
       x= "Factor 1 loading value",
       y= "Factor 1 loading rank") +
  geom_point(size=0.1, color =plotdata.rank.RNA.1$colorvalue) +
  geom_text_repel(size = 2, 
                  segment.size = 0.2, 
                  color =plotdata.rank.RNA.1$colorvalue,
                  nudge_x       = 1 - plotdata.rank.RNA.1$Weight,
                  direction     = "y",
                  hjust         = 1,
                  segment.color = "grey50") +
  theme_few()+
  add.textsize +
  scale_y_continuous(trans = "reverse") +
  add.textsize +
  theme(
    plot.title = element_markdown()
  )  +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank(), 
        axis.title.y=element_blank()
        )+
  xlim(c(-1,1))


#### Violins genes


f.violin.rna <- function(data , protein){
  
  ggplot(subset(data)  , aes(x = as.factor(condition), y =get(noquote(protein)))) +
    annotate("rect",
          xmin = 5 - 0.45,
             xmax = 6 + 0.6,
           ymin = -3, ymax = 15, fill = "lightblue",
           alpha = .4
  )+
  geom_violin(alpha = 0.9,aes(fill = condition))+
  geom_jitter(width = 0.05,size = 0.1, color = "black")+ 
  stat_summary(fun=median, geom="point", shape=95, size=2, inherit.aes = T, position = position_dodge(width = 0.9), color = "red")+
  theme_few()+
  ylab(paste0(protein)) +
  scale_x_discrete(labels = labels, expand = c(0.1,0.1), name = "Time aIg  (minutes)") +
  scale_y_continuous(expand = c(0,0), name = protein) +
  scale_color_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg ",)+ 
  scale_fill_manual(values = colorgradient6_manual2, labels = c(labels), name = "Time aIg ",) +
  add.textsize +
  theme(#axis.title.x=element_blank(),
        #axis.text.x=element_blank(),
        #axis.ticks.x=element_blank(),
        legend.position="none") 
}
rnadata <- as.data.frame(t(seu_combined_selectsamples@assays$SCT.RNA@scale.data)) %>%
  mutate(sample = rownames(t(seu_combined_selectsamples@assays$SCT.RNA@scale.data))) %>%
  left_join(meta.allcells)

enriched.geneset.posregBcell.fact1 <- c("NPM1", "NEAT1", "BTF3", "IGHM", "IGKC")

##Violin genes
for(i in enriched.geneset.posregBcell.fact1) {
assign(paste0("plot.violin.rna.", i), f.violin.rna(data = rnadata ,protein = i)) 
}

plot.violin.rna.NEAT1 <- plot.violin.rna.NEAT1 +
  labs(title = "Expression of upregulated genes\n")+
  add.textsize+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),)

plot.violin.rna.NPM1 <- plot.violin.rna.NPM1 +
  add.textsize+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),)