---

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()

```