--- a
+++ b/figures/Figure6.Rmd
@@ -0,0 +1,331 @@
+---
+title: "Figure 6"
+author: Tobias Roider
+date: "Last compiled on `r format(Sys.time(), '%d %B, %Y, %X')`"
+output: 
+  rmdformats::readthedown: 
+  
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r options, include=FALSE, warning = FALSE}
+
+library(knitr)
+opts_chunk$set(echo=TRUE, tidy=FALSE, include=TRUE, message=FALSE,
+               dpi = 100, cache = FALSE, warning = FALSE)
+opts_knit$set(root.dir = "../")
+options(bitmapType = "cairo")
+
+```
+
+# Read data, functions and packages
+```{r read data}
+
+source("R/ReadPackages.R")
+source("R/Functions.R")
+source("R/ReadData.R")
+source("R/ThemesColors.R")
+source("R/Helpers.R")
+
+```
+
+# Part 1
+## Frequencies
+```{r frequencies treg}
+ 
+df_freq_treg <- df_comb %>% 
+  add_prop(vars = c("PatientID", "Entity", "IdentI"), group.vars = 1) %>% 
+  fill_zeros(values_from = "Prop", names_from = "IdentI") %>% 
+  filter(IdentI %in% c(6,11))
+
+pvalues <- df_freq_treg %>% 
+  group_by(IdentI) %>% 
+  wilcox_test(data=., formula =  Prop ~ Entity, comparisons = list(c("FL", "rLN"), c("MZL", "rLN"))) %>% 
+  select(IdentI, Entity=group1, p) %>% 
+  mutate(p_s=format(p, scientific = TRUE, digits=1)) %>% 
+  mutate(p_f=case_when(p > 0.05 ~ "NA",
+                           p==0.05 ~ "0.05",
+                           p < 0.05 & p > 0.001 ~ as.character(round(p, 3)),
+                           p==0.001 ~ "0.001",
+                           p < 0.001 ~ p_s)) 
+
+p1 <- 
+  df_freq_treg %>% 
+  filter(IdentI %in% c(6)) %>% 
+  ggplot(aes(x=Entity, y=100*Prop))+
+  geom_boxplot(width=0.5, outlier.alpha = 0, size=0.25)+
+  ggbeeswarm::geom_beeswarm(size=0.75, shape=21, stroke=0.25, cex = 2.25, aes(fill=Entity))+
+  geom_text(inherit.aes = F, data = pvalues %>% filter(IdentI==6) %>% mutate(Y=c(38,45)),
+            aes(x=Entity, y=Y, label=p_f), size=2.5)+
+  scale_fill_brewer(palette = "Paired", limits=c("DLBCL", "MCL", "FL", "MZL", "rLN"))+
+  ggtitle(expression('T'[FH]))+
+  scale_y_continuous(limits = c(0,60), name="% of total T-cells (CITE-seq)")+
+  scale_x_discrete(limits=c("rLN", "DLBCL", "MCL", "FL", "MZL"))+
+  mytheme_1+
+  theme(legend.position = "none",
+        strip.background = element_rect(color=NA),
+        axis.title.x = element_blank(),
+        panel.border = element_rect(size=0.5),
+        plot.title = element_text(vjust = -1, color=colors_umap_cl[["6"]]),
+        axis.text.x = element_text(angle=45, hjust = 1),
+        panel.background = element_rect(fill=NA),
+        plot.margin = unit(c(0,0.1,0,0.25), "cm"))+
+  labs(tag = "A")
+
+p2 <- 
+  df_freq_treg %>% 
+  filter(IdentI %in% c(11)) %>% 
+  ggplot(aes(x=Entity, y=100*Prop))+
+  geom_boxplot(width=0.5, outlier.alpha = 0, size=0.25)+
+  ggbeeswarm::geom_beeswarm(size=0.75, shape=21, stroke=0.25, cex = 2.25, aes(fill=Entity))+
+  geom_text(inherit.aes = F, data = pvalues %>% filter(IdentI==11) %>% mutate(Y=c(16,17.5)),
+            aes(x=Entity, y=Y, label=p_f), size=2.5)+
+  scale_fill_brewer(palette = "Paired", limits=c("DLBCL", "MCL", "FL", "MZL", "rLN"))+
+  ggtitle(expression('T'[REG]~'EM'[2]))+
+  scale_y_continuous(limits = c(0,18.25), name="% of total T-cells (CITE-seq)")+
+  scale_x_discrete(limits=c("rLN", "DLBCL", "MCL", "FL", "MZL"))+
+  mytheme_1+
+  theme(strip.background = element_rect(color=NA),
+        plot.title = element_text(vjust = -1, color=colors_umap_cl[["11"]]),
+        axis.text.x = element_text(angle=45, hjust = 1),
+        panel.border = element_rect(size=0.5),
+        axis.title = element_blank(),
+        panel.background = element_rect(fill=NA),
+        plot.margin = unit(c(0,0.25,0,0.1), "cm"))
+
+```
+
+## Surface proteins
+```{r surface proteins}
+
+proteins_selected <- c("CD69"="CD69", 
+                       "CD25"="CD25", 
+                       "ICOS"="CD278", 
+                       "CD134"="CD134", 
+                       "CD161"="CD161", 
+                       "CCR4"="CD194", 
+                       "CCR5"="CD195", 
+                       "CXCR3"="CD183", 
+                       "CXCR5"="CD185", 
+                       "PD1"="CD279", 
+                       "CD38"="CD38", 
+                       "CD39"="CD39", 
+                       "TIGIT"="TIGIT")
+
+p3 <- 
+  left_join(percentageADT, meanADT) %>% 
+  filter(IdentI %in% c(6, 8, 11, 15, 13), Epitope %in% proteins_selected) %>% 
+  ggplot(aes(x=Epitope, y=IdentI, size=100*Prop, fill=Expression))+
+  geom_point(shape=21, stroke=0.1, color="grey45")+ 
+  scale_size_continuous(range=c(0, 4), 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=rev(as.character(c(6, 8, 13, 15, 11))), 
+                   labels=parse(text=rev(labels_cl_parsed[5:9]))
+                   )+
+  scale_x_discrete(limits=unname(proteins_selected), labels=names(proteins_selected))+
+  ggtitle("Protein level")+
+  coord_cartesian(clip = "off")+
+  theme_bw()+
+  theme(axis.title = element_blank(),
+        legend.position = "right",
+        axis.text.x = element_text(angle = 45, hjust = 1),
+        axis.text.y = element_blank(),
+        axis.text = element_text(size=7),
+        axis.ticks.y = element_blank(),
+        legend.text = element_text(size = 7, color="black"),
+        legend.title = element_text(size = 7, color="black", vjust = 0.5),
+        legend.box.margin=margin(-10,-10,-8,-8),
+        legend.key.height = unit(0.3, "cm"),
+        legend.key.width = unit(0.3, "cm"),
+        panel.border = element_rect(size=0.25),
+        plot.margin = unit(c(0,0.5,0,0.5), "cm"),
+        plot.tag = element_text(margin =  unit(c(0,1,0,0), "cm")),
+        plot.title = element_text(size = 7, color="black", vjust = -1, hjust = 0.5, face = "bold"),
+        panel.grid = element_blank())+
+  labs(tag = "B")
+
+for(i in 5:9) {
+  
+  p3 <- p3+
+    annotation_custom(grob = rectGrob(gp = gpar(fill=colors_umap_cl[as.character(cluster_order)[i]], col="white")), 
+                    ymin = rev(c(seq(0.5, 4.5, 1)+0.25))[i-4], 
+                    ymax = rev(c(seq(1.5, 5.5, 1)-0.25))[i-4],
+                    xmin = 0, xmax = -0.75)+
+    annotation_custom(grob = textGrob(label = labels_cl[[i]], rot = 0, hjust = 1, gp = gpar(cex=0.6)), 
+                      ymin = rev(c(seq(0.5, 4.5, 1)+0.25))[i-4], 
+                      ymax = rev(c(seq(1.5, 5.5, 1)-0.25))[i-4],
+                    xmin = -1, xmax = -1)
+}
+
+
+```
+
+## Assemble plot
+```{r assemble plot I, fig.height=2.5}
+
+p1+p2+p3+plot_layout(widths = c(1,1,1.45))
+ggsave(width = 18.5, height = 5.8, units = "cm", filename = "Figure6_p1.pdf")
+
+```
+
+# Part 2
+## Differentially expressed genes
+```{r differentially expressed genes}
+
+height.label <- 0.94
+position.label <- 1.1
+Idents(Combined_T) <- "IdentI"
+
+df_markers1115 <- FindMarkers(Combined_T, ident.1 = 11, ident.2 = c(15), test.use = "roc", assay = "integratedRNA") %>% 
+  rownames_to_column("Feature") %>% mutate(Assay="Gene")
+
+labels <- c("KLF2", "IKZF3", "IL21", "ASCL2", "IKZF2", "FOXP3", "CXCL13")
+
+df_tmp <- df_markers1115 %>% mutate(Label=ifelse(Feature %in% labels, Feature, NA)) %>% 
+  mutate(Label=gsub(Label, pattern = ".", fixed = T, replacement = ""))
+
+df_tmp1 <- df_tmp %>% filter(avg_log2FC<0)
+df_tmp2 <- df_tmp %>% filter(avg_log2FC>0)
+
+p4 <- 
+  ggplot()+
+  geom_point(data=df_tmp1, aes(x=avg_log2FC, y=power), alpha=ifelse(!is.na(df_tmp1$Label), 1, 0.25), stroke=0, size=1.25)+
+  geom_point(data=df_tmp2, aes(x=avg_log2FC, y=power), alpha=ifelse(!is.na(df_tmp2$Label), 1, 0.25), stroke=0, size=1.25)+
+  ggrepel::geom_text_repel(data=df_tmp1, aes(x=avg_log2FC, y=power, label=Label), show.legend = F, size=2.5, segment.size=0.25, xlim = c(-1.25, -1.75))+
+  ggrepel::geom_text_repel(data=df_tmp2, aes(x=avg_log2FC, y=power, label=Label), show.legend = F, size=2.5, segment.size=0.25, xlim = c(1.3, 1.5))+
+  geom_vline(xintercept = 0, linetype="dashed", size=0.25)+
+  scale_y_continuous(breaks = c(0.2, 0.4, 0.6, 0.8), limits=c(0.1, 0.85), name="2 x abs(AUC-0.5)")+
+  scale_x_continuous(name=expression('log'[2]~'fold change'), limits = c(-2.2, 2.2), expand = c(0,0))+
+  annotation_custom(grob = textGrob(label = expression('T'[REG]~'EM'[1]), hjust = 0.5, gp = gpar(cex=0.6, fontface="bold", col=colors_umap_cl["15"])), 
+                    xmin = -position.label, xmax = -position.label,
+                    ymin = height.label, ymax = height.label)+
+  annotation_custom(grob = textGrob(label = expression('T'[REG]~'EM'[2]), hjust = 0.5, gp = gpar(cex=0.6, fontface="bold", col=colors_umap_cl["11"])),
+                    xmin = position.label, xmax = position.label,
+                    ymin = height.label, ymax = height.label)+
+  mytheme_1+
+  coord_cartesian(clip = "off")+
+  theme(legend.position = c(0.17, 0.15),
+        legend.key.height = unit(units="cm", 0.3),
+        legend.box.spacing = unit(units="cm", 0.01),
+        legend.text = element_text(size=7),
+        panel.border = element_rect(size=0.25),
+        plot.margin = unit(c(0,0.25,0,0), "cm"),
+        legend.background = element_rect(fill = NA),
+        legend.box.margin=margin(-20,-20,-20,-20),
+        legend.key.width = unit(units="cm", 0.1))+
+    labs(tag = "C")
+
+```
+
+## Shared clonotypes
+```{r shared clonotypes}
+
+df_clonotypes_shared <- 
+  left_join(DFtotal_5prime %>% filter(!is.na(raw_clonotype_id)) %>% 
+            select(Barcode_fulla=Barcode_full, PatientID, refUMAP_1a=refUMAP_1, refUMAP_2a=refUMAP_2, IdentIa=IdentI, raw_clonotype_id) %>% distinct(),
+          DFtotal_5prime %>% filter(!is.na(raw_clonotype_id)) %>% 
+            select(Barcode_fullb=Barcode_full, PatientID, refUMAP_1b=refUMAP_1, refUMAP_2b=refUMAP_2, IdentIb=IdentI, raw_clonotype_id) %>% distinct()
+          ) %>% 
+  filter(Barcode_fulla!=Barcode_fullb) %>% 
+  filter(IdentIa!=IdentIb)
+
+df_subset <- 
+  df_clonotypes_shared %>% 
+  add_entity() %>% 
+  filter(IdentIb==11) %>% 
+  filter(refUMAP_1b<1, refUMAP_2b>4)
+
+label6 <- paste0(100*round(nrow(df_subset %>% filter(IdentIa==6))/nrow(df_subset), 3), " %")
+label14 <- paste0(100*round(nrow(df_subset %>% filter(IdentIa==14))/nrow(df_subset), 3), " %")
+label5 <- paste0(100*round(nrow(df_subset %>% filter(IdentIa==5))/nrow(df_subset), 3), " %")
+
+p5 <- ggplot()+
+  geom_point_rast(data=DFtotal_5prime,
+                  aes(x=refUMAP_1, y=refUMAP_2, fill=IdentI), size=0.25, 
+                  alpha=ifelse(DFtotal_5prime$IdentI==11, 0.75, 0.05), stroke=0, shape=21)+
+  geom_curve(data= df_subset, 
+             aes(x=refUMAP_1a, y=refUMAP_2a, xend=refUMAP_1b, yend=refUMAP_2b, color=IdentIa,
+                 group=paste(raw_clonotype_id, PatientID)), curvature = -0.4, size=0.1, alpha=0.4)+
+  geom_text(inherit.aes = F, aes(x=-0.75, y=8.5, label=label6), color=colors_umap_cl[["6"]], size=2.5)+
+  geom_text(inherit.aes = F, aes(x=4.5, y=6.75, label=label14), color=colors_umap_cl[["14"]], size=2.5)+
+  geom_text(inherit.aes = F, aes(x=8, y=1, label=label5), color=colors_umap_cl[["5"]], size=2.5)+
+  scale_fill_manual(values = colors_umap_cl, guide="none")+
+  scale_color_manual(values = colors_umap_cl, guide="none")+
+  coord_cartesian(clip = "off")+
+  labs(
+    x="refUMAP-1",
+    y="refUMAP-2",
+    title="Paired clonotypes of <span style='color:#08306B'>T<sub>REG</sub> EM<sub>2</sub></span>")+
+  mytheme_1+
+  theme(panel.border = element_rect(size=0.25),
+        plot.title = element_textbox_simple(size = 7, width = NULL, padding = margin(1.25, 0, 1, 0), 
+                                            lineheight = 1.25, halign=0.5, face = "plain"),
+        plot.margin = unit(c(0,0.25,0,0), units = "cm"))+
+  labs(tag = "D")
+
+
+```
+
+## Association with FL grade
+```{r tumor grading}
+
+p6 <- df_freq %>% 
+  left_join(., df_meta %>% select(PatientID, FL_Grade, Entity) %>% distinct, by="PatientID") %>% 
+  filter(Population==11, Entity=="FL") %>% 
+  ggplot(aes(x=FL_Grade, y=RNA))+
+  geom_boxplot(size=0.25, width=0.3, outlier.alpha = 0)+
+  ggbeeswarm::geom_beeswarm(size=0.8, shape=21, stroke=0.25, cex = 2.75, fill="grey65")+
+  stat_compare_means(comparisons = list(c("1/2", "3A")), vjust = -0.35, label.y = c(16.75),
+                     size=2.5, tip.length = 0.02, bracket.size = 0.25)+
+  ggtitle(expression('T'[REG]~'EM'[2]))+
+  xlab("Grade")+
+  ylim(0,18.25)+
+  ylab("% of total T-cells")+
+  mytheme_1+
+  theme(legend.position = "none",
+        strip.background = element_rect(color=NA),
+        panel.border = element_rect(size=0.5),
+        plot.title = element_text(vjust = -1, color=colors_umap_cl[["11"]]),
+        panel.background = element_rect(fill=NA),
+        plot.margin = unit(c(0,0.25,0,0.25), "cm"))+
+    labs(tag = "E")
+
+p7 <- df_freq %>% 
+  left_join(., df_meta %>% select(PatientID, FL_Grade, Entity) %>% distinct, by="PatientID") %>% 
+  filter(Population==6, Entity=="FL") %>% 
+  ggplot(aes(x=FL_Grade, y=RNA))+
+  geom_boxplot(size=0.25, width=0.3, outlier.size = 1, outlier.alpha = 0, outlier.shape = 21, outlier.fill = "grey70")+
+  ggbeeswarm::geom_beeswarm(size=0.8, shape=21, stroke=0.25, cex = 2.75, fill="grey65")+
+  stat_compare_means(comparisons = list(c("1/2", "3A")), vjust = -0.35, label.y = c(59.5),
+                     size=2.5, tip.length = 0.02, bracket.size = 0.25)+
+  ggtitle(expression('T'[FH]))+
+  ylim(0,65)+
+  xlab("Grade")+
+  ylab("% of total T-cells")+
+  mytheme_1+
+  theme(legend.position = "none",
+        strip.background = element_rect(color=NA),
+        panel.border = element_rect(size=0.5),
+        plot.title = element_text(vjust = -1, color=colors_umap_cl[["6"]]),
+        panel.background = element_rect(fill=NA),
+        plot.margin = unit(c(0,0.25,0,0.25), "cm"))+
+    labs(tag = "F")
+
+```
+
+## Assemble plot
+```{r compose plot II, fig.height=2.8}
+
+p4+p5+p6+p7+plot_layout(widths = c(1.05,1,0.3,0.3))
+
+#ggsave(width = 18.5, height = 6.5, units = "cm", filename = "Figure6_p2.pdf")
+
+```
+
+# Session info
+```{r session info}
+
+sessionInfo()
+
+```