genes_selected <- 
  c("MKI67",
    "CCR7", 
    "KLF2",
    "TCF7", 
    "TOX",
    "TOX2", 
    "ASCL2",
    "FOXP3", 
    "IKZF3",
    "GZMA", 
    "GZMK", 
    "CCL5", 
    "NKG7")

DefaultAssay(Combined_T) <- "integratedRNA"

perc_expr <- 
  FetchData(Combined_T, slot = "counts", vars = c("IdentI", paste0("rna_", genes_selected))) %>% 
  mutate(IdentI=as.factor(IdentI)) %>% 
  mutate_if(.predicate = is.numeric, .funs = ~ifelse(isZero(.), 1, 0)) %>% 
  pivot_longer(cols = 2:ncol(.), names_to = "Gene") %>% 
  group_by(IdentI, Gene) %>% 
  count(value) %>% 
  mutate(Prop=n/sum(n)) %>% 
  filter(value==0) %>% 
  select(-value, -n) %>% 
  mutate(Gene=substr(Gene, 5, nchar(.)))

DefaultAssay(Combined_T) <- "integratedRNA"

plot_genex <- 
  FetchData(Combined_T, slot = "data", vars = c("IdentI", paste0(genes_selected))) %>% 
  mutate(IdentI=factor(IdentI, levels = rev(cluster_order))) %>% 
  group_by(IdentI) %>% 
  summarise_all(mean) %>% 
  pivot_longer(cols = 2:ncol(.), names_to = "Gene") %>% 
  group_by(Gene) %>% 
  mutate(value=(value-min(value))/(max(value)-min(value))) %>% 
  left_join(., perc_expr) %>% 
  ggplot(aes(x=Gene, y=IdentI, size=100*Prop, fill=value))+
  geom_point(shape=21, stroke=0.1, color="grey45")+ 
  scale_size_continuous(range=c(0, 3), name="% pos. cells", limits=c(0, 100))+
  scale_fill_gradientn(name="Expression", colours = brewer.pal(5, "BuGn"), limits=c(0,1))+
  scale_y_discrete(limits=factor(rev(cluster_order)), labels=rev(unlist(labels_cl)))+
  scale_x_discrete(limits=genes_selected)+
  geom_hline(yintercept = c(1.5, 5.5, 9.5, 10.5, 13.5), linetype="solid", size=0.25, alpha=0.1)+
  ggtitle("RNA level")+
  coord_cartesian(clip = "off")+
  theme_bw()+
  mytheme_1+
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size=7),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(0.25,0.35,0,2), "cm"))

lines <- c(1, 5, 9, 10, 13, 14)

for(i in 1:length(cluster_order)) {
  
  plot_genex <- plot_genex+
   annotation_custom(grob = rectGrob(gp = gpar(fill=colors_umap_cl[as.character(rev(cluster_order)[i])], lex=1, col="white")), 
                    ymin = seq(0.5, length(cluster_order)-0.5, 1)[i], 
                    ymax = seq(1.5, length(cluster_order)+0.5, 1)[i],
                    xmin = 0, xmax = -1.5)+
   annotation_custom(grob = textGrob(label = paste0("C", c(14:1)[i]), gp = gpar(cex=0.6, col=ifelse(i %in% c(6,7,11,14), "white", "black"))),
                    ymin = seq(0.5, length(cluster_order)-0.5, 1)[i], 
                    ymax = seq(1.5, length(cluster_order)+0.5, 1)[i],
                    xmin = 0, xmax = -1.5)
}

for(i in 1:length(lines)) {
  
  plot_genex <- plot_genex+
    annotation_custom(grob = textGrob(label = rev(labels_celltypes_expr)[[i]], rot = 0, hjust = 1, gp = gpar(cex=0.6)), 
                      ymin = c(0,lines)[i]+0.5,
                      ymax = c(lines)[i]+0.5,
                      xmin = -1.65, xmax = -1.65)+
    annotation_custom(grob = linesGrob(gp = gpar(col="white", lex=3)),
                      ymin = c(0,lines)[i]+0.5,
                      ymax = c(0,lines)[i]+0.5,
                      xmin = -0.01, xmax = -1.5)
}

plot_genex <- plot_genex+labs(tag = "B")+
 theme(plot.tag.position = c(-0.25,1))

proteins_selected <- 
  c("CD4"="CD4", 
    "CD8a"="CD8a", 
    "CD45RA"="CD45RA", 
    "CD45RO"="CD45RO", 
    "CD95"="CD95", 
    "CD62L"="CD62L", 
    "CD127"="CD127", 
    "CD69"="CD69", 
    "CD38"="CD38", 
    "CD25"="CD25", 
    "ICOS"="CD278", 
    "CXCR5"="CD185", 
    "CD31"="CD31", 
    "KLRG1"="KLRG1", 
    "CD244"="CD244", 
    "PD1"="CD279", 
    "TIM3"="CD366"
    )

plot_protex <- 
  left_join(percentageADT, meanADT) %>% 
  filter(Epitope %in% proteins_selected) %>% 
  ggplot(aes(x=Epitope, y=IdentI, size=100*Prop, fill=Expression))+
  geom_point(shape=21, stroke=0.1, color="grey45")+ 
  geom_hline(yintercept = c(1.5, 5.5, 9.5, 10.5, 13.5), linetype="solid", size=0.25, alpha=0.1)+
  scale_size_continuous(range=c(0, 3), name="% pos. cells", limits=c(0, 100))+
  scale_fill_gradientn(name="Expression", colours = brewer.pal(5, "BuGn"), limits=c(0,1))+
  scale_y_discrete(limits=factor(rev(cluster_order)), labels=rev(unlist(labels_cl)))+
  scale_x_discrete(limits=proteins_selected, labels=names(proteins_selected))+
  ggtitle("Protein level")+
  coord_cartesian(clip = "off")+
  theme_bw()+
  mytheme_1+
  theme(axis.title = element_blank(),
        legend.position = "right",
        axis.text.x = element_text(angle = 45, hjust = 1, size=7),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.text = element_text(size = 7, color="black"),
        legend.title = element_text(size = 7, color="black", vjust = 0.8),
        legend.key.height = unit(0.3, "cm"),
        legend.key.width = unit(0.3, "cm"),
        legend.box.spacing = unit(0.1, "cm"),
        plot.margin = unit(c(0.25,0,0,0.15), "cm"),
        plot.tag.position = c(-0.025,1))+
        labs(tag = "C")

tfs_selected <- c("TCF7"="tfactivity_TCF7-E",  
                  "FOXP3"="tfactivity_FOXP3-E",
                  "ASCL2"="tfactivity_ASCL2-E", 
                  "KLF2"="tfactivity_KLF2-E")

df_tfs <- 
  FetchData(Combined_T, vars = c("Barcode_full", unname(tfs_selected))) %>% 
  left_join(df_comb %>% select(IdentI, Barcode_full, CellType), .) %>% 
  pivot_longer(cols =4:ncol(.)) %>% 
  mutate(name=gsub(name, pattern = "tfactivity_|-E", replacement = "")) %>% 
  mutate(name=factor(name, levels = names(tfs_selected))) %>% 
  group_by(name, IdentI) %>% 
  summarise(Mean=mean(value, na.rm=T)) %>% 
  group_by(name) %>% 
  mutate(Mean=2*((Mean-min(Mean))/(max(Mean)-min(Mean)))-1) 

plot_tfact <- 
  ggplot(df_tfs, aes(y=as.character(IdentI), x=name, fill=Mean))+
  geom_tile()+ 
  scale_fill_gradientn(name="TF activity", colours = colorRampPalette(colors = c("#762a83", "#f7f7f7", "#1b7837"))(100))+
  geom_vline(xintercept = seq(1.5, 4.5, 1), color="white", size=0.25)+
  geom_hline(yintercept = seq(1.5, 14.5, 1), color="white", size=0.25)+
  scale_y_discrete(limits=rev(factor(cluster_order)), expand = c(0,0))+
  scale_x_discrete(expand = c(0,0))+
  ggtitle("TF activity")+
  coord_fixed(clip = "off")+
  theme_bw()+
  mytheme_1+
  theme(axis.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, size=7),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.border = element_rect(size=0.25),
        plot.background = element_rect(fill = NA, color=NA),
        legend.position = "right",
        legend.text = element_text(size = 7, color="black"),
        legend.key.height = unit(0.3, "cm"),
        legend.key.width = unit(0.3, "cm"),
        legend.box.spacing = unit(0.1, "cm"),
        plot.margin = unit(c(0.25,0,0,0.65), "cm"),
        plot.tag.position = c(-0.2,1))+
        labs(tag = "D")

lines <- c(1, 5, 9, 10, 13, 14)

for(i in 1:length(cluster_order)) {
  
  plot_tfact <- plot_tfact+
   annotation_custom(grob = rectGrob(gp = gpar(fill=colors_umap_cl[as.character(rev(cluster_order)[i])], lex=1, col="white")), 
                    ymin = seq(0.5, length(cluster_order)-0.5, 1)[i], 
                    ymax = seq(1.5, length(cluster_order)+0.5, 1)[i],
                    xmin = 0, xmax = -1.5)+
   annotation_custom(grob = textGrob(label = paste0("C", c(14:1)[i]), gp = gpar(cex=0.6, col=ifelse(i %in% c(6,7,11,14), "white", "black"))),
                    ymin = seq(0.5, length(cluster_order)-0.5, 1)[i], 
                    ymax = seq(1.5, length(cluster_order)+0.5, 1)[i],
                    xmin = 0, xmax = -1.5)
}

for(i in 1:length(lines)) {
  
  plot_tfact <- plot_tfact+
    annotation_custom(grob = linesGrob(gp = gpar(col="white", lex=3)),
                      ymin = c(0,lines)[i]+0.5,
                      ymax = c(0,lines)[i]+0.5,
                      xmin = -0.01, xmax = -1.5)
}

plot_tfact

#ggsave(plot_tfact, filename = "Figure1_p3.pdf", width = 5, height = 7.35, units = "cm")