---
title: "Supplementary Figures"
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, cache.lazy = 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")

```

# Supplementary Figure 1
## Patient characteristics
```{r patient characteristics SF1 part 1}

p1 <- 
  df_meta %>% 
  ggplot(aes(x=Order, y=Tcells, fill=Entity))+
  geom_bar(stat = "identity", color="white", width=0.5, size=0.25)+
  geom_text(aes(x=6, y=84, label=paste0("n = ", nrow(df_meta))), check_overlap = T, size=2.75)+
  scale_y_continuous(name="% T-cells", limits=c(0, 90), expand = c(0.025, 0.025))+
  scale_x_discrete(limits=as.character(1:101))+
  scale_fill_manual(values = colors_characteristics)+
  mytheme_1+
  coord_cartesian(clip = "off")+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=7),
        axis.ticks.x = element_blank(),
        plot.tag = element_text(margin = unit(c(0,0,0,-0.75), "cm")),
        plot.margin = unit(c(0.1,0.1,0,1), "lines"))

cex_ <- 0.6
pos_ <- -0.25
size_ <- 2.4

p_ann <- 
  ggplot()+
  geom_tile(data=df_meta, aes(y=0, x=Order, fill=Entity))+
  geom_tile(data=df_meta, aes(y=-1.75, x=Order, fill=Status))+
  geom_tile(data=df_meta, aes(y=-3.5, x=Order, fill=Sex))+
  geom_text(data=df_meta, aes(y=-5.25, x=Order, label=`CITEseq`), size=size_)+
  geom_text(data=df_meta, aes(y=-7, x=Order, label=`TCRseq`), size=size_)+
  geom_text(data=df_meta, aes(y=-8.75, x=Order, label=`MultiIF`), size=size_)+
  geom_text(data=df_meta, aes(y=-10.5, x=Order, label=`MultiFlow`), size=size_)+
  annotation_custom(grob=textGrob(label = "Entity", gp = gpar(cex=cex_), just="right"), xmin = pos_, xmax = pos_, ymin = 0, ymax = 0)+
  annotation_custom(grob=textGrob(label = "Collection", gp = gpar(cex=cex_), just="right"), xmin = pos_, xmax = pos_, ymin = -1.75, ymax = -1.75)+
  annotation_custom(grob=textGrob(label = "Sex", gp = gpar(cex=cex_), just="right"), xmin = pos_, xmax = pos_, ymin = -3.5, ymax = -3.5)+
  annotation_custom(grob=textGrob(label = "CITE-seq", gp = gpar(cex=cex_), just="right"), xmin = pos_, xmax = pos_, ymin = -5.25, ymax = -5.25)+
  annotation_custom(grob=textGrob(label = "TCR-seq", gp = gpar(cex=cex_), just="right"), xmin = pos_, xmax = pos_, ymin = -7, ymax = -7)+
  annotation_custom(grob=textGrob(label = "Multi-IF", gp = gpar(cex=cex_), just="right"), xmin = pos_, xmax = pos_, ymin = -8.75, ymax = -8.75)+
  annotation_custom(grob=textGrob(label = "Multi-flow", gp = gpar(cex=cex_), just="right"), xmin = pos_, xmax = pos_, ymin = -10.5, ymax = -10.5)+
   scale_x_discrete(limits=as.character(1:101))+
  scale_y_continuous(expand = c(0,0), name=NULL)+
  geom_vline(xintercept = seq(0.5, nrow(df_meta), 1), size=1, color="white")+
  scale_fill_manual(values = colors_characteristics)+
  coord_cartesian(clip = "off")+
  theme_void()+
  theme(legend.position = "none",
        plot.margin = unit(c(0,0.1,0,0.1), "lines"))

plot_legend1 <- 
  ggplot()+
  geom_tile(data=df_meta, aes(y=4, x=Order, fill=Entity))+
  scale_x_discrete(expand = c(0,0), name=NULL)+
  scale_y_discrete(expand = c(0,0), name=NULL)+
  geom_vline(xintercept = seq(0.5, 51.5, 1), color="white")+
  scale_fill_manual(values = colors_characteristics, limits=df_meta$Entity %>% unique())+
  guides(fill=guide_legend(direction = "horizontal", keywidth = 0.3, keyheight = 0.4))+
  theme_void()+
  mytheme_1+
  theme(legend.position = "bottom",
        legend.title = element_text(size=7, face = "bold"),
        legend.text = element_text(size=7))

legend1 <- get_legend(plot_legend1)

plot_legend2 <-
  ggplot()+
  geom_tile(data=df_meta, aes(y=4, x=Order, fill=Sex))+
  scale_x_discrete(expand = c(0,0), name=NULL)+
  scale_y_discrete(expand = c(0,0), name=NULL)+
  geom_vline(xintercept = seq(0.5, 51.5, 1), color="white")+
  scale_fill_manual(values = colors_characteristics, limits=df_meta$Sex %>% unique())+
  guides(fill=guide_legend(direction = "horizontal", keywidth = 0.3, keyheight = 0.4))+
  theme_void()+
  mytheme_1+
  theme(legend.position = "bottom",
        legend.title = element_text(size=7, face = "bold"),
        legend.text = element_text(size=7))

legend2 <- get_legend(plot_legend2)

plot_legend3 <- 
  ggplot()+
  geom_tile(data=df_meta, aes(y=4, x=Order, fill=Status))+
  scale_x_discrete(expand = c(0,0), name=NULL)+
  scale_y_discrete(expand = c(0,0), name=NULL)+
  geom_vline(xintercept = seq(0.5, 51.5, 1), color="white")+
  scale_fill_manual(values = colors_characteristics, name="Collection",
                    limits=filter(df_meta, Status!="NA")$Status %>% unique())+
  guides(fill=guide_legend(direction = "horizontal", keywidth = 0.3, keyheight = 0.4))+
  theme_void()+
  mytheme_1+
  theme(legend.position = "bottom",
        legend.title = element_text(size=7, face = "bold"),
        legend.text = element_text(size=7))

legend3 <- get_legend(plot_legend3)

```

## Associations with overall T-cell frequencies
```{r patient characteristics SF1 part 2}

# Sex
p2 <- df_meta %>% 
  filter(!Entity %in% c("CLL", "rLN")) %>%
  mutate(Tcells=ifelse(is.na(Tcells), TcellsDx, Tcells)) %>% 
  filter(Tcells>1) %>% 
  ggplot(aes(x=Sex, y=Tcells))+
  geom_boxplot(width=0.35, size=0.25)+
  #ggbeeswarm::geom_beeswarm(size=0.8, shape=21, stroke=0.25, cex = 2.75, aes(fill=Entity))+
  stat_compare_means(comparisons = list(c("F", "M")), label.y = 78, size=2.5, bracket.size = 0.25)+
  scale_fill_manual(values = hue_pal()(5))+
  scale_x_discrete(labels=c("Female", "Male"))+
  ylim(0,90)+
  ylab("T-cells in %")+
  mytheme_1+
  theme(axis.title.x = element_blank(),
        plot.tag = element_text(margin = unit(c(0,0,0,-0.75), "cm")))+
  labs(tag = "B")

# Status
p3 <- df_meta %>% 
  filter(!Entity %in% c("CLL", "rLN")) %>%
  mutate(Tcells=ifelse(is.na(Tcells), TcellsDx, Tcells)) %>% 
  filter(Tcells>1) %>% 
  ggplot(aes(x=Status, y=Tcells))+
  geom_boxplot(width=0.35, size=0.25)+
  #ggbeeswarm::geom_beeswarm(size=0.8, shape=21, stroke=0.25, cex = 2.75, aes(fill=Entity))+
  stat_compare_means(comparisons = list(c("Initial diagnosis", "Relapse")), label.y = 78, bracket.size = 0.25, size=2.5)+
  scale_x_discrete(labels=c("Initial \ndiagnosis", "Relapse"))+
  scale_fill_manual(values = hue_pal()(5))+
  ylim(0,90)+
  ylab("T-cells in %")+
  mytheme_1+
  theme(axis.title.x = element_blank())+
  labs(tag = "C")

# Age
p5 <- df_meta %>% 
  filter(!Entity %in% c("CLL", "rLN")) %>%
  mutate(Tcells=ifelse(is.na(Tcells), TcellsDx, Tcells)) %>% 
  filter(Tcells>1) %>% 
  ggplot(aes(x=Age, y=Tcells))+
  geom_point(size=1.25, shape=21, stroke=0.25, color="white", aes(fill=Entity))+
  stat_cor(aes(label=..r.label..),  size=2.5, label.y = 82)+
  scale_fill_manual(values = colors_characteristics, name=NULL, limits=c("DLBCL", "FL", "MCL", "MZL"))+
  ylim(0,90)+
  ylab("T-cells in %")+
  mytheme_1+
  theme(legend.position = "right",
        axis.title.x = element_text(vjust = 6.5),
        legend.box.margin = unit(c(0,0,0,-0.38), "cm"),
        legend.key.width = unit("cm", x = 0.15),
        legend.key.height = unit("cm", x = 0.35))+
  labs(tag = "D")

# Entity
df_tmp <- df_meta %>% 
  filter(!Entity %in% c("CLL", "rLN")) %>%
  mutate(Tcells=ifelse(is.na(Tcells), TcellsDx, Tcells)) %>% 
  filter(Tcells>1)
  
p6 <- ggplot(df_meta, aes(x=Entity, y=Tcells))+
  geom_boxplot(width=0.4, size=0.25)+
  ggbeeswarm::geom_beeswarm(size=0.75, shape=21, stroke=0.25, cex = 2.25, aes(fill=Entity))+
  stat_compare_means(data=df_tmp %>% filter(!Entity %in% c("rLN")), 
                     label.y = 83, label.x = 2.5, size=2.5, hjust=0.5, bracket.size = 0.25)+
  scale_fill_manual(values = colors_characteristics, name=NULL, 
                    limits=c("DLBCL", "FL", "MCL", "MZL", "rLN"))+
  geom_segment(data = NULL, aes(x=1, xend=4, y=80, yend=80), size=0.1)+
  ylim(0,90)+
  ylab("T-cells in %")+
  mytheme_1+
  theme(axis.title.x = element_blank())+
  labs(tag = "E")

```

## Assemble plot
```{r assemble SF1, fig.height=5.5}

# Plot
wrap_plots(p1+labs(tag = "A")+p_ann+plot_layout(nrow = 2, heights = c(1, 0.5)))/
  plot_spacer()/
  wrap_plots(p2+p3+p5+p6+plot_layout(nrow = 1, widths = c(0.55,0.55,1,1)))+
  plot_layout(heights = c(1.25,0.075,0.7))

#ggsave(width = 18.5, height = 12, units = "cm", filename = "SF1.pdf")

```

## Legend
```{r legend SF1, fig.height=0.5}

# Legend
as_ggplot(legend1)+as_ggplot(legend2)+as_ggplot(legend3)+plot_layout(nrow = 1, widths = c(1, 0.5, 0.75))
ggsave(width = 13, height = 1, units = "cm", filename = "S1_p1_legend.pdf")

```

# Supplementary Figure 2
## Confusion Matrix
```{r, confusion SF2, fig.height=4.5}

confmat <- confusionMatrix(GBclasses_surface$test, GBclasses_surface$predict)

conf_freq <- confmat$table %>% 
  data.frame() %>% 
  group_by(Reference) %>% 
  mutate(Prop=Freq/sum(Freq)) 
  
plot_surface <- ggplot(conf_freq, aes(x=Reference, y=Prediction, fill=Prop, label=round(Prop, 2)))+
  geom_tile()+
  geom_text(data=conf_freq %>% filter(Prop>0.4), color="white", size=1.75)+
    scale_fill_gradientn(colours = colorRampPalette(RColorBrewer::brewer.pal(9, "BuPu"))(100), name="Sensitivity",
                       limits=c(0,0.9), breaks=seq(0,0.8,0.2))+
  geom_hline(yintercept = seq(1.5, 13.5), color="black",linetype="solid", size=0.25, alpha=0.25)+
  geom_vline(xintercept = seq(1.5, 13.5), color="black",linetype="solid", size=0.25, alpha=0.25)+
  scale_y_discrete(expand = c(0,0), limits=factor(cluster_order), labels=unlist(labels_cl))+
  scale_x_discrete(expand = c(0,0), limits=factor(cluster_order), labels=unlist(labels_cl))+
  coord_fixed()+
  ggtitle("Surface marker only")+
  mytheme_1+
  theme(legend.position = "none",
        axis.text.x = element_text(angle=45, hjust = 1))+
  labs(tag = "A")

confmat_plus <- confusionMatrix(GBclasses_surfaceplus$test, GBclasses_surfaceplus$predict)

conf_freq <- confmat_plus$table %>% 
  data.frame() %>% 
  group_by(Reference) %>% 
  mutate(Prop=Freq/sum(Freq)) 
  
plot_surfaceplus <- ggplot(conf_freq, aes(x=Reference, y=Prediction, fill=Prop, label=round(Prop, 2)))+
  geom_tile()+
  geom_text(data=conf_freq %>% filter(Prop>0.4), color="white", size=1.75)+
  scale_fill_gradientn(colours = colorRampPalette(RColorBrewer::brewer.pal(9, "BuPu"))(100), name="Sensitivity",
                       limits=c(0,0.9), breaks=seq(0,0.8,0.2))+
  geom_hline(yintercept = seq(1.5, 13.5), color="black",linetype="solid", size=0.25, alpha=0.25)+
  geom_vline(xintercept = seq(1.5, 13.5), color="black",linetype="solid", size=0.25, alpha=0.25)+
  scale_y_discrete(expand = c(0,0), limits=factor(cluster_order), labels=unlist(labels_cl))+
  scale_x_discrete(expand = c(0,0), limits=factor(cluster_order), labels=unlist(labels_cl))+
  coord_fixed()+
  ggtitle("Surface plus FoxP3, IKZF3, and Ki67")+
  mytheme_1+
  theme(legend.position = "right",
        legend.key.height = unit(0.35, "cm"),
        legend.key.width = unit(0.28, "cm"),
        legend.margin = margin(c(0,0.2,0,0), unit = "cm"),
        axis.text.x = element_text(angle=45, hjust = 1))+
  labs(tag = "B")

plot_surface+plot_surfaceplus+plot_layout(guides = "collect", widths = c(1,1))

#ggsave(width = 19, height = 9, units = "cm", filename = "SF2.pdf")

```

# Supplementary Figure 3
## Frequencies overview
```{r freq overview SF3}

df_pop

mat_complete <- rbind(
  df_facs %>% 
    left_join(., df_meta %>% select(PatientID, `CITEseq`)) %>% 
    filter(`CITEseq`=="-") %>% 
    select(PatientID, Population, Prop=FACS),
  df_freq %>% 
    mutate(RNA=ifelse(is.nan(RNA), 0, RNA)) %>% 
    select(PatientID, Population, Prop=RNA)) %>% 
  filter(Population %in% df_pop$Population) %>% 
  left_join(., df_pop) %>% 
  select(-Population) %>% 
  left_join(., df_meta %>% select(PatientID, Tcells)) %>% 
  mutate(Prop=Tcells*Prop/100) %>% 
  select(-Tcells) %>% 
  pivot_wider(names_from = "IdentI", values_from = "Prop", values_fill = 0) %>% 
  column_to_rownames("PatientID") 

cl_order <- c(6,1,2,9,8,13,15,11,12,16,3,5,14,19)
names(labels_cl_parsed) <- as.character(cluster_order)

df_freqPlot <- 
  mat_complete %>% 
  rownames_to_column("PatientID") %>% 
  pivot_longer(cols=2:ncol(.), names_to = "IdentI", values_to = "Prop") %>% 
  add_entity() %>% 
  mutate(Entity=factor(Entity, levels = c("rLN", "DLBCL", "MCL", "FL", "MZL"))) %>% 
  mutate(IdentI=factor(IdentI, levels=cl_order)) %>% 
  mutate(Label=factor(IdentI, levels=cl_order, labels = labels_cl_parsed[as.character(cl_order)])) %>% 
  group_by(Entity, IdentI) %>% 
  mutate(outlier = (Prop > quantile(Prop, 0.75) + IQR(Prop) * 1.5) | (Prop < quantile(Prop, 0.25) - IQR(Prop) * 1.5))

df_medianLines <- df_freqPlot %>% 
  filter(Entity=="rLN") %>% 
  group_by(IdentI, Label) %>% 
  summarise(MedianProp=median(Prop))

df_freqPlot_pvalues <- df_freqPlot %>% 
  group_by(IdentI) %>% 
  wilcox_test(data=., formula = Prop ~ Entity, detailed = T, ref.group = "rLN") %>% 
  adjust_pvalue(method = "BH") %>% 
  select(IdentI, Entity=group2, p.adj, estimate) %>% 
  mutate(Entity=factor(Entity, levels = c("rLN", "DLBCL", "MCL", "FL", "MZL"))) %>% 
  mutate(p.adj=ifelse(p.adj>0.05, NA, p.adj)) %>% 
  mutate(p.adj_s=format(p.adj, scientific = TRUE, digits=1)) %>% 
  mutate(p.adj_f=case_when(p.adj > 0.05 ~ "NA",
                           p.adj==0.05 ~ "0.05",
                           p.adj < 0.05 & p.adj > 0.001 ~ as.character(round(p.adj, 3)),
                           p.adj==0.001 ~ "0.001",
                           p.adj < 0.001 ~ p.adj_s)) %>% 
  filter(!is.na(p.adj)) %>% 
  left_join(., df_freqPlot %>% select(IdentI,  Label) %>% distinct) %>% 
  left_join(., data.frame(IdentI=factor(cl_order), height=c(28, 28, 21, 21, 10, 10, 12, 12, 14, 14, 50, 50, 12, 12)))

```

## Plot
```{r freq overview, fig.height=6, fig.width=5.2}

p <- list()

for(i in c(1:7)){
  
  y <- list(c(1,6),c(2,9),c(8,13),c(15,11),c(12,16),c(3,5),c(14,19))[[i]]
  ylim <- c(38,30,13,15.5,20,65,15)
  
  p[[i]] <- 
    ggplot(data=df_freqPlot %>% filter(IdentI %in% y) %>% 
             mutate(Label=factor(Label, levels = labels_cl_parsed[as.character(y)])), 
             aes(y=Prop, x=Entity, fill=IdentI))+
    geom_hline(data=df_medianLines %>%filter(IdentI %in% y), aes(yintercept=MedianProp),
               size=0.25, linetype="dashed", color="grey60")+
    geom_boxplot(width=0.4, outlier.shape = 21, outlier.size = 1, outlier.color = "white", 
                 outlier.alpha = 0, show.legend = F, size=0.25)+
    ggbeeswarm::geom_beeswarm(data = function(x) dplyr::filter_(x, ~ outlier), cex = 3, stroke=0.25, 
                              groupOnX = TRUE, shape = 21, size = 1, color = "white", alpha = 1)+
    geom_text(data=df_freqPlot_pvalues %>% filter(IdentI %in% y), 
              inherit.aes = F, aes(y=height, x=Entity, label=p.adj_f), hjust=0.1, size=2.3, angle=45)+
    scale_fill_manual(values = colors_umap_cl)+
    scale_y_continuous(name="% of total cells", limits=c(0,ylim[i]))+
    scale_x_discrete(expand = c(0.17,0.17))+
    facet_wrap(~Label, ncol = 2, labeller = label_parsed)+
    mytheme_1+
    theme(axis.title.x = element_blank(),
          strip.background = element_rect(color=NA),
          plot.margin = unit(c(0.1,0.25,0,0), units = "cm"),
          panel.border = element_rect(size = 0.5),
          axis.text.x = element_text(angle=45, hjust = 1))
  
  if(i!=7){
    p[[i]] <- p[[i]]+
      theme(axis.text.x = element_blank(),
            axis.ticks.x = element_blank())
  }
  
  if(!i %in% c(4,5)){
    p[[i]] <- p[[i]]+
      theme(axis.title.y = element_blank())
    }
  
}

plot_freq <- wrap_plots(p, ncol = 2)
plot_freq
ggsave(width = 18.5, height = 15, units = "cm", filename = "SF3.pdf")

```


# Supplementary Figure 4
## Lasso model (beta coefficients)
```{r lasso model SF4}

mat_complete <- rbind(
  df_facs %>% 
    left_join(., df_meta %>% select(PatientID, `CITEseq`)) %>% 
    filter(`CITEseq`=="-") %>% 
    select(PatientID, Population, Prop=FACS),
  df_freq %>% 
    mutate(RNA=ifelse(is.nan(RNA), 0, RNA)) %>% 
    select(PatientID, Population, Prop=RNA)) %>% 
  filter(Population %in% df_pop$Population) %>% 
  left_join(., df_pop) %>% 
  select(-Population) %>% 
  pivot_wider(names_from = "IdentI", values_from = "Prop", values_fill = 0) %>% 
  column_to_rownames("PatientID") 

total <- mat_complete %>% 
  rownames_to_column("PatientID") %>% 
  left_join(., df_meta %>% select(PatientID, Tcells)) %>% 
  add_entity() %>% 
  mutate_if(is.numeric, .funs = ~./100)

cell_types <- total %>% select(-Entity, -PatientID) %>% colnames()

gt <- my_glmnet(total)

limits_coefs=c(cluster_order, "Tcells", "(Intercept)")

plot_coefs <- gt$coefs %>% 
  mutate(Entity=factor(Entity, levels = c("rLN", "DLBCL", "MCL", "FL", "MZL"))) %>% 
  ggplot(aes(x=beta, y=cell_type,  fill=cell_type, color=cell_type))+
  geom_hline(yintercept = c(1.5, 5.5, 9.5, 10.5, 13.5)+2, linetype="solid", size=0.25, alpha=0.1)+
  geom_bar(stat = "identity", width = 0.5)+
  scale_y_discrete(limits=rev(limits_coefs), labels=c("Intercept", "T cells", rev(unlist(unname(labels_cl)))))+
  scale_color_manual(limits=rev(limits_coefs), values = c("black", "black", rev(unname(colors_umap_cl[as.character(cluster_order)]))))+
  scale_fill_manual(limits=rev(limits_coefs), values = c("black", "black", rev(unname(colors_umap_cl[as.character(cluster_order)]))))+
  facet_wrap(~Entity, ncol=5)+
  mytheme_1+
  theme(axis.title.y = element_blank(),
        panel.border = element_rect(size = 0.5))

```

## Multivariate Model
```{r multivariate model SF4, fig.height=7}

models <- list()
entities <- c("Overall", "FL", "MZL", "MCL")

for(i in colnames(mat_complete)){
  for(e in entities) {
  ent <- e
  if(e=="Overall"){ent <- c("FL", "MZL", "MCL", "DLBCL")}
    
models[[paste0(i, "_", e)]] <- mat_complete %>% 
  rownames_to_column("PatientID") %>% 
  pivot_longer(cols=2:ncol(.), names_to = "IdentI", values_to = "Prop") %>% 
  left_join(., df_meta %>% select(PatientID, Status, Pretreatment, Sex, Entity, Department, Age, Subtype, Tcells, `CITEseq`) %>% distinct, 
            by="PatientID") %>%
  filter(Entity %in% ent) %>% 
  filter(IdentI==i) %>% 
  glm(formula = Prop ~ Status + Age + Sex, 
      family = "gaussian") %>% summary()

models[[paste0(i, "_", e)]] <- models[[paste0(i, "_", e)]]$coefficients %>% 
  `colnames<-`(c("Estimate", "Std.error", "t.value", "p.value")) %>% 
  data.frame() %>% 
  dplyr::mutate(IdentI=as.character(i)) %>% 
  rownames_to_column("Parameter")

}
}

for(i in colnames(mat_complete)){
  ent <- "DLBCL"
  e <- "DLBCL"
  models[[paste0(i, "_", e)]] <- mat_complete %>% 
    rownames_to_column("PatientID") %>% 
    pivot_longer(cols=2:ncol(.), names_to = "IdentI", values_to = "Prop") %>% 
    left_join(., df_meta %>% select(PatientID, Status, Pretreatment, Sex, Entity, Department, Age, Subtype, Tcells, `CITEseq`) %>% distinct, 
              by="PatientID") %>%
    filter(Entity %in% ent) %>% 
    filter(IdentI==i) %>% 
    glm(formula = Prop ~ Status + Age + Sex + Subtype, 
        family = "gaussian") %>% summary()

models[[paste0(i, "_", e)]] <- models[[paste0(i, "_", e)]]$coefficients %>% 
  `colnames<-`(c("Estimate", "Std.error", "t.value", "p.value")) %>% 
  data.frame() %>% 
  dplyr::mutate(IdentI=as.character(i)) %>% 
  rownames_to_column("Parameter")

}

p.value <- bind_rows(models, .id = "model") %>% 
  filter(Parameter!="(Intercept)") %>% 
  mutate(Entity=strsplit(model, split = "_") %>% sapply(., "[[", 2)) %>% 
  group_by(IdentI, Entity) %>% 
  mutate(p.adj=p.adjust(p.value, method = "BH")) %>% 
  mutate(sign=ifelse(Estimate<0 & p.adj<0.05, "-", NA)) %>% 
  mutate(sign=ifelse(Estimate>0 & p.adj<0.05, "+", sign)) %>% 
  mutate(Entity=factor(Entity, levels = c("Overall", "DLBCL", "MCL", "FL", "MZL")))


plot_multi <- ggplot()+
  geom_point(data=p.value %>% filter(Entity=="Overall"), aes(x=IdentI, y=Parameter, size=-log10(p.adj), fill=t.value), shape=21, stroke=0.25,
             position = position_nudge(y=0.16))+
  geom_point(data=p.value %>% filter(Entity!="Overall"), aes(x=IdentI, y=Parameter, size=-log10(p.adj), fill=t.value, group=Entity), shape=21, stroke=0.25, 
             position = ggpp::position_dodgenudge(y = -0.14, width = 0.78))+
  geom_text(data=p.value %>% filter(Entity=="Overall", p.adj<0.05), aes(x=IdentI, y=Parameter, label=round(p.adj, 3)), hjust=0.5, nudge_y = 0.4, size=2.5)+
  geom_text(data=p.value %>% filter(Entity=="Overall", p.adj<0.05), aes(x=IdentI, y=Parameter, label=sign), hjust=0.5, nudge_y = 0.18, size=2.5)+
  geom_text(data=p.value %>% filter(Entity!="Overall"), aes(x=IdentI, y=Parameter, label=sign, group=Entity), size=2.5, 
             position = ggpp::position_dodgenudge(y = -0.11, width = 0.78))+
  scale_x_discrete(limits=factor(cluster_order), labels=unlist(labels_cl))+
  scale_size_continuous(range=c(0.25,3.5), limits=c(0,3), name=expression('-log'[10]~'p'))+
  scale_fill_gradientn(colors=brewer.pal(9, name = "RdBu"), limits=c(-5, 5), breaks=c(-3,0,3), name="Statistic")+
  guides(fill=guide_colorbar(barwidth = 2.5))+
  guides(size=guide_legend(keywidth = 0.3))+
  geom_hline(yintercept = seq(1,6,1), color="grey90", size=0.25)+
  geom_vline(xintercept = seq(1.5,13.5,1), color="grey90", size=0.25)+
  scale_y_discrete(limits=rev(c("StatusRelapse",  "Age", "SexM", "Subtypenon-GCB")),
                   labels=rev(c("<b style='color:#B2182B'>Initial diagnosis </b><br>vs. <b style='color:#2166AC'>Relapse</b>",
                                "<b style='color:#B2182B'>Younger </b>vs.<br><b style='color:#2166AC'>Older</b>",
                                "<b style='color:#B2182B'>Female </b>vs.<br><b style='color:#2166AC'>Male</b>",
                                "<b style='color:#B2182B'>GCB</b> vs.<br><b style='color:#2166AC'>non-GCB</b><br>[only DLBCL]")))+
  mytheme_1+
  theme(axis.title = element_blank(),
        legend.position = "top",
        plot.title = element_text(margin = unit(c(0,0,0,0), units = "cm")),
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(size=0.5),
        legend.key.height = unit(x = 0.3, units = "cm"),
        legend.box.margin = unit(c(0,-8.5,-0.35,0), "cm"),
        plot.margin = unit(c(0,0.25,0,0.25), "cm"),
        legend.key.width = unit(x = 0.35, units = "cm"),
        axis.text.y = element_markdown(hjust = 1,  lineheight = 1.2))+
  labs(tag = "B")

plot_coefs+labs(tag = "A")+plot_multi+plot_layout(heights = c(1.2,1))

ggsave(width = 18, height = 15, units = "cm", filename = "SF4.pdf")

```

# Supplementary Figure 5
## Proportions of T-cell subsets in 5' scRNA versus CITE-seq
### Data handling
```{r data handling SF5}

df_freq_5prime <- DFtotal_5prime %>% 
  select(Barcode_full, IdentI, PatientID) %>% 
  distinct() %>% 
  add_prop(vars = c("PatientID", "IdentI"), group.vars = 1) %>% 
  mutate(Class="5prime") %>% 
  filter(PatientID %in% unique(df_comb$PatientID))

df_freq_3prime <- df_comb %>% 
  add_prop(vars = c( "PatientID", "IdentI"), group.vars = 1) %>% 
  mutate(Class="3prime") %>% 
  filter(PatientID %in% DFtotal_5prime$PatientID)

df_freq_prime <- rbind(df_freq_3prime, df_freq_5prime) %>% 
  mutate(Prop=100*Prop) %>% 
  pivot_wider(names_from = "Class", values_from = "Prop", values_fill = 0) 

```

### Correlation plots
```{r cor plots SF5}

this_theme <- theme(plot.margin = unit(c(0,0.25,0.1,0), units = "cm"),
  plot.title = element_text(vjust = -1),
  axis.title = element_blank(),
  axis.text = element_text(size=7, color="black"))

cor_plots_prime <- list()

cor_plots_prime[["TFH"]] <- 
  df_freq_prime %>% 
  filter(IdentI==6) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["6"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["6"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,60), breaks = c(0, 20, 40, 60))+
  scale_y_continuous(limits = c(0,60), breaks = c(0, 20, 40, 60))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["6"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme

cor_plots_prime[["TPR"]] <- df_freq_prime %>% 
   filter(IdentI==14) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["14"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["14"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,12), breaks = c(0, 4, 8, 12))+
  scale_y_continuous(limits = c(0,12), breaks = c(0, 4, 8, 12))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["14"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme+
  theme(axis.title.y = element_text(size=7))

cor_plots_prime[["TDN"]] <- df_freq_prime %>% 
  filter(IdentI==19) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["19"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["19"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,8), breaks = c(0, 2, 4, 6, 8))+
  scale_y_continuous(limits = c(0,8), breaks = c(0, 2, 4, 6, 8))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["19"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  theme_axis_sub3

cor_plots_prime[["THCM2"]] <- df_freq_prime %>% 
  filter(IdentI==9) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["9"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["9"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,15), breaks = c(0, 5, 10, 15))+
  scale_y_continuous(limits = c(0,15), breaks = c(0, 5, 10, 15))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["9"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme

cor_plots_prime[["THCM1"]] <- df_freq_prime %>% 
  filter(IdentI==2) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["2"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["2"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,20), breaks = c(0, 6, 12, 18))+
  scale_y_continuous(limits = c(0,20), breaks = c(0, 6, 12, 18))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["2"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme


cor_plots_prime[["THNaive"]] <- df_freq_prime %>% 
  filter(IdentI==1) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["1"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["1"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,45), breaks = c(0, 15, 30, 45))+
  scale_y_continuous(limits = c(0,45), breaks = c(0, 15, 30, 45))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["1"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme

cor_plots_prime[["TREGCM1"]] <- df_freq_prime %>% 
  filter(IdentI==8) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["8"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["8"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["8"]])+
    scale_x_continuous(limits = c(0,25), breaks = c(0, 8, 16, 24))+
  scale_y_continuous(limits = c(0,25), breaks = c(0, 8, 16, 24))+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme

cor_plots_prime[["TREGEM2"]] <- df_freq_prime %>% 
  filter(IdentI==11) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["11"]],size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["11"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,16), breaks = c(0, 5, 10, 15))+
  scale_y_continuous(limits = c(0,16), breaks = c(0, 5, 10, 15))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["11"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  theme_axis_sub3

cor_plots_prime[["TREGEM1"]] <- df_freq_prime %>% 
  filter(IdentI==15) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["15"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["15"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,20), breaks = c(0, 6, 12, 18))+
  scale_y_continuous(limits = c(0,20), breaks = c(0, 6, 12, 18))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["15"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme+
  theme(axis.title.y = element_text(size=7),
        axis.title.x = element_text(size=7))

cor_plots_prime[["TREGCM2"]] <- df_freq_prime %>% 
  filter(IdentI==13) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["13"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["13"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,20), breaks = c(0, 6, 12, 18))+
  scale_y_continuous(limits = c(0,20), breaks = c(0, 6, 12, 18))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["13"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  this_theme

cor_plots_prime[["TTOXNaive"]] <-df_freq_prime %>% 
  filter(IdentI==12) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["12"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["12"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,27), breaks = c(0, 8, 16, 24))+
  scale_y_continuous(limits = c(0,27), breaks = c(0, 8, 16, 24))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["12"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  theme_axis_sub3
 
cor_plots_prime[["TTOXEM3"]] <- df_freq_prime %>% 
  filter(IdentI==5) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color="#b50923", size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color="#b50923", se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,80), breaks = c(0, 25, 50, 75))+
  scale_y_continuous(limits = c(0,80), breaks = c(0, 25, 50, 75))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["5"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  theme_axis_sub3

cor_plots_prime[["TTOXEM1"]] <- df_freq_prime %>% 
  filter(IdentI==3) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color="#fea044", size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x,na.rm = T,
              color="#fea044", se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,50), breaks = c(0, 15, 30, 45))+
  scale_y_continuous(limits = c(0,50), breaks = c(0, 15, 30, 45))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["3"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  theme_axis_sub3

cor_plots_prime[["TTOXEM2"]] <- df_freq_prime %>% 
  filter(IdentI==16) %>% 
  ggplot(aes(y=`5prime`, x=`3prime`))+
  geom_point(color=colors_umap_cl[["16"]], size=0.65, alpha=0.75)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_umap_cl[["16"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson",size=2.5, label.x.npc = c(0.45), label.y.npc = c(0.9))+
  scale_x_continuous(limits = c(0,15), breaks = c(0, 4, 8, 12))+
  scale_y_continuous(limits = c(0,15), breaks = c(0, 4, 8, 12))+
  labs(x="3' scRNA - %", y="5' scRNA - %")+
  ggtitle(labels_cl[["16"]])+
  theme_bw()+
  coord_fixed()+
  mytheme_1+
  theme_axis_sub3

cor_plots_prime_all <- cor_plots_prime$TPR+labs(tag ="A")+
  cor_plots_prime$THNaive+cor_plots_prime$THCM1+cor_plots_prime$THCM2+cor_plots_prime$TFH+
  cor_plots_prime$TREGCM1+cor_plots_prime$TREGCM2+cor_plots_prime$TREGEM1+cor_plots_prime$TREGEM2+
  cor_plots_prime$TTOXNaive+cor_plots_prime$TTOXEM1+cor_plots_prime$TTOXEM2+cor_plots_prime$TTOXEM3+
  cor_plots_prime$TDN+
  plot_layout(nrow = 2)

df_freq_prime %>% 
  group_by(IdentI) %>% 
  summarise(R=cor.test(`3prime`, `5prime`)$estimate) %>% pull(R) %>% median()

```

## CD4 and CD8 expression in clonal T-cells
```{r CD4 and CD8 expression in SF5}

lim_CD4 <- 0.3
lim_CD8 <- 0.6

df_clone_expr <- lapply(sobjs_T_5prime, function(sobj){
  FetchData(sobj, vars=c("CD4", "CD8A")) %>% 
              data.frame() %>% 
              rownames_to_column("Barcode_full")
  }) %>% bind_rows() %>% remove_rownames()

df_clone_CD4 <- DFtotal_5prime %>% 
  select(Barcode_full, PatientID, refUMAP_1, refUMAP_2, Entity, IdentI, raw_clonotype_id) %>% 
  distinct() %>% 
  add_count(IdentI, PatientID, raw_clonotype_id) %>% 
  left_join(., df_clone_expr) %>% 
  mutate(IdentI_new=factor(IdentI, levels = cluster_order, labels = labels_cl_parsed)) %>% 
  filter(IdentI %in% c(3,5)) %>% 
  mutate(Expanded=n>2) %>% 
  mutate(CD4CD8=case_when(CD4 > lim_CD4 & CD8A < lim_CD8 ~ "CD4pos",
                          CD4 < lim_CD4 & CD8A > lim_CD8 ~ "CD8pos",
                          CD4 > lim_CD4 & CD8A > lim_CD8 ~ "CD4CD8pos",
                          CD4 < lim_CD4 & CD8A < lim_CD8 ~ "CD4CD8neg"))
  
df_clone_CD4_freq <- df_clone_CD4 %>% 
  drop_na() %>% 
  add_prop(vars = c("IdentI", "Expanded", "CD4CD8"), group.vars = c(1,2)) %>% 
  fill_zeros(names_from = "IdentI", values_from = "Prop") %>% 
  mutate(Prop=round(Prop, 2)) %>% 
  left_join(., data.frame(CD4CD8=c("CD4CD8neg", "CD4pos", "CD4CD8pos",  "CD8pos"), 
                          CD8=c(-0.35, -0.35, 3.7, 3.7),
                          CD4=c(-0.35, 1.9, 1.9, -0.35))) %>% 
    mutate(IdentI_new=factor(IdentI, levels = cluster_order, labels = labels_cl_parsed)) 
  
plot_expr1 <- ggplot()+
  geom_point(data=df_clone_CD4 %>% filter(n<3), inherit.aes = F, aes(x=CD4, y=CD8A, color=IdentI), 
             position = position_jitter(width = 0.1, height = 0.1),  na.rm = T,
             size=0.25, alpha=0.25, stroke=0)+
  scale_color_manual(values = colors_umap_cl, guide="none")+
  scale_fill_manual(values = colors_umap_cl, guide="none")+
  geom_text(data=df_clone_CD4_freq %>% filter(Expanded==F), inherit.aes = F, aes(x=CD4, y=CD8, label=Prop), size=2.5, alpha=1)+
  geom_hline(yintercept = 0.6, linetype="dashed", size=0.25)+
  geom_vline(xintercept = 0.3, linetype="dashed", size=0.25)+
  scale_x_continuous(breaks = c(0,0.5,1,1.5,2), limits=c(-0.5, 2.4))+
  scale_y_continuous(breaks = c(0,1,2,3), limits=c(-0.5, 4))+
  labs(x="<i>CD4</i> - RNA expression",
       y="<i>CD8</i> - RNA expression",
       title="Clone size < 3",
       tag = "B")+
  facet_wrap(~IdentI_new, nrow = 1, labeller = label_parsed)+
  mytheme_1+
  theme(axis.title.x = element_textbox(size=7, halign = 0.5, margin = unit(units = "cm", c(0.1,0,0,0))),
        axis.title.y = element_textbox(size=7, orientation = "left-rotated", margin = unit(units = "cm", c(0,0,0.1,0))),
        panel.border = element_rect(size=0.4))

plot_expr2 <- ggplot()+
  geom_point(data=df_clone_CD4 %>% filter(n>2), inherit.aes = F, aes(x=CD4, y=CD8A, color=IdentI),
             position =  position_jitter(width = 0.1, height = 0.1), na.rm = T,
             size=0.25, alpha=0.25, stroke=0)+
  scale_color_manual(values = colors_umap_cl, guide="none")+
  scale_fill_manual(values = colors_umap_cl, guide="none")+
  scale_size_continuous(range=c(1, 5), limits=c(3, 50), breaks=c(3, 20, 35, 50),
                        labels=c("3", "20", "35", "> 50"), name = "Clonotype size")+
  geom_hline(yintercept = 0.6, linetype="dashed", size=0.25)+
  geom_vline(xintercept = 0.3, linetype="dashed", size=0.25)+
  geom_text(data=df_clone_CD4_freq %>% filter(Expanded==T), inherit.aes = F, aes(x=CD4, y=CD8, label=Prop), size=2.5, alpha=1)+
  facet_wrap(~IdentI_new, nrow = 1, labeller = label_parsed)+
  scale_x_continuous(breaks = c(0,0.5,1,1.5,2), limits=c(-0.5, 2.4))+
  scale_y_continuous(breaks = c(0,1,2,3), limits=c(-0.5, 4))+
  labs(x="<i>CD4</i> - RNA expression",
       y="<i>CD8</i> - RNA expression",
       title="Clone size > 2",
       tag = "C")+
  facet_wrap(~IdentI_new, nrow = 1, labeller = label_parsed)+
  mytheme_1+
  theme(axis.title.x = element_textbox(size=7, halign = 0.5, margin = unit(units = "cm", c(0.1,0,0,0))),
        axis.title.y = element_textbox(size=7, orientation = "left-rotated", margin = unit(units = "cm", c(0,0,0.1,0))),
        panel.border = element_rect(size=0.4))

```

## Assemble plot
```{r assemble SF5}

cor_plots_prime_all/wrap_plots(plot_expr1+plot_spacer()+plot_expr2+plot_layout(widths = c(1,0.1,1)))+
  plot_layout(heights = c(1.6,1))

#ggsave(width = 18.5, height = 12.5, units = "cm", filename = "SF5.pdf")

```

## TCR diversity
### Read
```{r TCR diversity read}

# Single cell T-cell receptor data read by immunarch package
# RNA–seq, epitope and TCR raw and processed data have been deposited in the Gene Expression Omnibus (GEO) under accession codes GSE252608 and GSE252455.
DF_immunarchTCR <- repLoad(list.files(path = "countMatrices", pattern = "TCRrep", full.names = T))
DF_immunarchTCR$meta$Sample <- strsplit(DF_immunarchTCR$meta$Sample, split = "_") %>% sapply("[[", 1)
names(DF_immunarchTCR$data) <- DF_immunarchTCR$meta$Sample

```

### Plot
```{r TCR diversity plot, fig.height=3.5}

plots <- list()
for(i in unique(DFtotal_5prime$PatientID)){
set.seed(substr(i, 4,7) %>% as.numeric()+5)
plots[[i]] <- 
  DFtotal_5prime %>% filter(PatientID==i) %>% 
  select(Barcode_full, PatientID, raw_clonotype_id) %>% 
  distinct() %>% 
  dplyr::count(raw_clonotype_id) %>% 
  drop_na() %>% 
  mutate(Prop=n/sum(n)) %>% 
  dplyr::arrange(-n) %>% 
  mutate(Cumsum=cumsum(n)) %>% 
  mutate(Max=sum(n)) %>% 
  filter(Cumsum < 0.1*Max) %>% 
  mutate(new_id=as.character(1:nrow(.))) %>% 
  mutate(PatientID=i) %>% 
  add_entity() %>% 
  mutate(PatientID_new=paste0(PatientID, " (", Entity, ")")) %>% 
  ggplot(aes(x=new_id, y=n))+
  geom_segment(aes(x=new_id, xend=new_id, y=0, yend=n), size=0.2)+
  facet_wrap(~PatientID_new)+
  xlab("Clonotype ID")+
  scale_y_continuous(name =  "Clonotype size")+
  scale_x_discrete(limits=sample(as.character(1:195), 195), name="Unique clonotype ID")+
  mytheme_1+
  theme(legend.position = "none",
        panel.border = element_rect(size=0.4),
        strip.background = element_rect(colour = NA),
        axis.title.x = element_text(vjust = 5),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

}

imm_raref <- repDiversity(DF_immunarchTCR$data, "raref", .verbose = F) %>% 
  rename(PatientID=Sample) %>% 
  add_entity() %>% 
  filter(!PatientID %in% c("LN0256", "LN0367"))

tops <- imm_raref %>% 
  group_by(PatientID) %>% 
  top_n(n=1, Size) %>% 
  #mutate(Mean=ifelse(PatientID=="LN0302", 5.0, Mean)) %>% 
  mutate(Size=ifelse(PatientID=="LN0132", 88, Size)) %>% 
  mutate(Size=ifelse(PatientID=="LN0110", 72, Size)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0110", 19, Mean)) %>% 
  mutate(Size=ifelse(PatientID=="LN0259", 90, Size)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0259", 20.75, Mean)) %>% 
  mutate(Size=ifelse(PatientID=="LN0264", 55, Size)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0264", 13.5, Mean)) %>% 
  mutate(Size=ifelse(PatientID=="LN0217", Size-10, Size)) %>% 
  mutate(Size=ifelse(PatientID=="LN0144", 145, Size)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0193", Mean+1, Mean)) %>% 
  mutate(Size=ifelse(PatientID=="LN0198", 115, Size)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0198", Mean+0.2, Mean)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0278", Mean+1, Mean)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0144", Mean+0.75, Mean)) %>% 
  mutate(Mean=ifelse(PatientID=="LN0302", Mean-1.5, Mean)) %>% 
  mutate(Size=ifelse(PatientID=="LN0302", 50, Size)) %>% 
  filter(!PatientID %in% c("LN0417", "LN0104", "LN0249")) # Manuel labelling in Powerpoint

plot_rare <- 
  imm_raref %>%
  ggplot(aes(x=Size, y=Mean, group=PatientID, color=Entity))+
  geom_line(linetype="solid", size=0.25)+
  geom_label(data=tops, aes(x=Size, y=Mean, label=PatientID),  size=2.25, 
                            show.legend = F, fill="white", color="white", 
                            label.padding = unit(units = "cm", 0.02), label.size = 0)+
  geom_text(data=tops, aes(x=Size, y=Mean, label=PatientID), size=2, show.legend = F)+
  guides(color=guide_legend(override.aes = list(size=0.35, linetype="solid")))+
  scale_color_brewer(palette = "Paired", limits=c("DLBCL", "MCL", "FL", "MZL", "rLN"))+
  ylab("Estimated diversity")+
  xlab("Clone size")+
  #xlim(0, 70)+
  mytheme_1+
  theme(#legend.position = "top",
        legend.title = element_blank(),
        panel.grid = element_blank(),
        legend.spacing.x = unit("cm", x = 0.05),
        legend.box.margin = unit(c(0,0,-0.95,0), "cm"),
        panel.border = element_rect(size=0.4),
        legend.key.height = unit("cm", x = 0.36),
        legend.key.width = unit("cm", x = 0.5))+
  labs(tag = "D")

plot_rare+wrap_plots((plots$LN0132+theme(axis.title.x = element_blank())+labs(tag = "E"))/plots$LN0217+labs(tag = "F"))+
  plot_layout(widths = c(1,1.7))
  
#ggsave(width = 18, height = 8, units = "cm", filename = "SF5.pdf")

```

# Supplementary Figure 6
## T-cell exhaustion UMAP
### Calculate exhaustion module
```{r calculate exhaustion module}

exhausted_cells <- ttox@meta.data %>%
  mutate(Exhausted=ifelse( Pseudotime>=24, "yes", "no")) %>% 
  filter(Exhausted=="yes") %>% rownames()

Combined_T@meta.data$Exhausted <- Combined_T@meta.data %>% 
  mutate(Exhausted=Barcode_full %in% exhausted_cells) %>% 
  pull(Exhausted)

Idents(Combined_T) <- "Exhausted"

module_exhausted <- FindMarkers(Combined_T, ident.1 = "TRUE", assay = "integratedRNA", test.use = "roc") %>% 
  rownames_to_column("Gene") %>% 
  mutate(Module=paste0(Gene, ifelse(myAUC>0.5, "+", "-")),
         Assay="RNA")

module_exhausted_prot <- FindMarkers(Combined_T, ident.1 = "TRUE", assay = "integratedADT", test.use = "roc") %>% 
  rownames_to_column("Gene") %>% 
  mutate(Module=paste0(Gene, ifelse(myAUC>0.5, "+", "-")),
         Assay="Protein")

#module_exhausted --> Supplementary Table 5
#WriteXLS::WriteXLS(rbind(module_exhausted, module_exhausted_prot), ExcelFileName = "SuppTable5.xlsx")

module_exhausted <- list(module_exhausted$Module)
names(module_exhausted) <- "exhausted"

Combined_T <- UCell::AddModuleScore_UCell(Combined_T, features = module_exhausted)

```

### Plot exhaustion module
```{r plot exhaustion module}

set.seed(1)
plot_exh <- FetchData(Combined_T, vars = c("wnnUMAP_1", "wnnUMAP_2", "exhausted_UCell")) %>% 
  sample_frac(0.2) %>% 
  ggplot(aes(x=wnnUMAP_1, y=wnnUMAP_2, fill= exhausted_UCell))+
  ggrastr::geom_point_rast(size=0.25, stroke=0, shape=21, raster.dpi = 600, alpha=0.75)+
  scale_fill_gradientn(colours = brewer.pal(n = 9, name = "YlOrRd")[2:9],
                       name="Score")+
  xlab("wnnUMAP-1")+
  ylab("wnnUMAP-2")+
  ggtitle("Exhaustion signature")+
  mytheme_1+
  theme(panel.border = element_rect(size = 0.2),
        axis.title.x = element_text(margin = unit(units = "cm", c(-0.75,0,0,0))),
        #legend.margin = margin(c(0,0,0,-0.35), unit = "cm"),
        legend.box.margin = unit(c(0,0,0,-0.35), units = "cm"),
        legend.title = element_text(size=6),
        legend.text = element_text(size=6),
        legend.position = "right",
        plot.title = element_text(face = "plain", vjust = -0.5),
        panel.background = element_rect(fill=NA),
        legend.key.height = unit(units="cm", 0.2),
        legend.key.width = unit(units="cm", 0.15),
        legend.box.background = element_rect(fill=NA, color=NA),
        legend.background = element_rect(fill=NA, color=NA)
        )+
  labs(tag="A")

```

## Association with cell-of-origin in DLBCL
```{r SF6 part 1}

### Schmitz
p1 <- left_join(df_surv_schmitz, df_ttoxcompl_schmitz) %>% 
  drop_na() %>% 
  ggplot(aes(x=Subtype, y=Exhausted/Absolute, group=Subtype))+
  geom_boxplot(outlier.alpha = 0, width=0.4, size=0.25)+
  stat_compare_means(size=2.25, vjust = 1, aes(label=paste0("p = ", ..p.format..)))+
  stat_compare_means(comparisons = list(c("ABC", "GCB")), size=2.25)+
  scale_y_continuous(limits=c(0, 0.65), name = "Exhausted T-cells")+
  #ggtitle("Schmitz et al. 2018")+
  xlab("Cell-of-origin")+
  mytheme_1+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        #axis.title.x = element_text(margin = unit(units = "cm", c(-0.1,0,0,0))),
        plot.title = element_text(face = "plain", vjust = -0.5))+
  labs(tag = "B")

### Chapuy
p2 <- left_join(df_surv_chapuy, df_ttoxcompl_chapuy) %>% 
  drop_na() %>% 
  ggplot(aes(x=Subtype, y=Exhausted/Absolute, group=Subtype))+
  geom_boxplot(outlier.alpha = 0, width=0.4, size=0.25)+
  stat_compare_means(size=2.25, vjust = 1, aes(label=paste0("p = ", ..p.format..)))+
  scale_y_continuous(limits=c(0, 0.65), name="Exhausted T-cells")+
  #ggtitle("Chapuy et al. 2018")+
  xlab("Cell-of-origin")+
  mytheme_1+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        #axis.title.x = element_text(margin = unit(units = "cm", c(-0.1,0,0,0))),
        plot.title = element_text(face = "plain", vjust = -0.5))+
  labs(tag = "C")

```

## Association with genetic Subtype
```{r SF6 part 2}

### Schmitz
p3 <- left_join(df_surv_schmitz, df_ttoxcompl_schmitz) %>% 
  drop_na() %>% 
  ggplot(aes(x=GenSubtype, y=Exhausted/Absolute, group=GenSubtype))+
  geom_boxplot(outlier.alpha = 0, width=0.4, size=0.25)+
  stat_compare_means(size=2.25, vjust = 1, aes(label=paste0("p = ", ..p.format..)))+
  scale_y_continuous(limits = c(0, 0.65), name="Exhausted T-cells")+
  #ggtitle("Schmitz et al. 2018")+
  xlab("Cluster")+
  mytheme_1+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(face = "plain", vjust = -0.5),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        #axis.title.x = element_text(margin = unit(units = "cm", c(-0.1,0,0,0)))
        )

### Chapuy
p4 <- left_join(df_surv_chapuy, df_ttoxcompl_chapuy) %>% 
  drop_na() %>% 
  mutate(Cluster=paste0("C", Cluster)) %>% 
  ggplot(aes(x=Cluster, y=Exhausted/Absolute, group=Cluster))+
  geom_boxplot(outlier.alpha = 0, width=0.4, size=0.25)+
  stat_compare_means(size=2.25, vjust = 1, aes(label=paste0("p = ", ..p.format..)))+
  scale_y_continuous(limits=c(0, 0.65), name="Exhausted T-cells")+
  #ggtitle("Chapuy et al. 2018")+
  xlab("Cluster")+
  mytheme_1+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(face = "plain", vjust = -0.5),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        #axis.title.x = element_text(margin = unit(units = "cm", c(-0.65,0,0,0)))
        )

```

## Assemble plot
### Part 1
```{r assemble plot SF6 part 1, fig.height=2.5}

plot_exh+p1+p3+p2+p4+plot_layout(nrow = 1, widths = c(1.4,0.75,1,0.75,1))

#ggsave(width = 19, height = 6, units = "cm", filename = "SF6_p1.pdf")

```


## Associations with genetic features
### Mutations
```{r mutations SF6}

pairs <- list(c(1:42), c(43:85))

df_mut <- df_snvs_chapuy %>% 
  filter(Description=="Mutation", value %in% c(0,2)) %>% 
  mutate(value=factor(value, levels = c("0", "2"), labels = c("wt", "mut"))) %>% 
  group_by(Name, value) %>%
  dplyr::summarise(MeanExhausted=mean(Exhausted/Absolute)) %>% 
  mutate(MeanExhausted=ifelse(MeanExhausted>0.15, 0.15, MeanExhausted),
         Group=case_when(Name %in% unique(.$Name)[pairs[[1]]] ~ "group1",
                         Name %in% unique(.$Name)[pairs[[2]]] ~ "group2")) 

p.values <- df_snvs_chapuy %>% 
  filter(Description=="Mutation", value %in% c(0,2)) %>% 
  mutate(Exhausted=Exhausted/Absolute) %>% 
  compare_means(formula = Exhausted ~ value, group.by = "Name", p.adjust.method = "BH") %>% 
  filter(p<0.05) %>% 
  left_join(., df_mut %>% select(Name, Group), by="Name")

p5 <- df_mut %>% 
  ggplot(aes(x=value, y=Name, fill=MeanExhausted))+
  geom_tile()+
  scale_fill_gradientn(name="Expression", colours = brewer.pal(5, "PRGn"), limits=c(0,0.17))+
  geom_vline(xintercept = 1.5, color="black", size=0.25)+
  geom_segment(data=data.frame(y=seq(1.5, 42.5,1)), inherit.aes = F, 
               aes(y = y, yend = y, x=0.05, xend=2.8), color="white", size=0.25)+
  geom_text(data = p.values, inherit.aes = F, aes(y=Name, x=2.8, label=round(p, 2)), size=2.5)+
  facet_wrap(~Group, ncol=3, scales = "free_y")+
  ggtitle("Mutations")+
  coord_cartesian(clip = "off")+
  mytheme_1+
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        axis.text.x = element_text(angle=45, hjust=1, size=7),
        axis.ticks = element_blank(),
        axis.text.y = element_text(size=7, margin = unit(c(0,-0.4,0,0), units = "cm")),
        plot.margin = unit(units = "cm", c(0,1.25,0,0)),
        axis.title = element_blank(),
        panel.border = element_blank())+
  labs(tag = "E")

```

### Copy number gain
```{r copy number SF6}

df_gain <- df_snvs_chapuy %>% 
  filter(Description=="CN gain") %>% 
  mutate(value=factor(value, levels = c("0", "1", "2"), labels = c("wt", "gain", "gain"))) %>% 
  group_by(Name, value) %>%
  dplyr::summarise(MeanExhausted=mean(Exhausted/Absolute)) 

p.values <- df_snvs_chapuy %>% 
  filter(Description=="CN gain") %>% 
  mutate(Exhausted=Exhausted/Absolute) %>% 
  mutate(value=factor(value, levels = c("0", "1", "2"), labels = c("wt", "gain", "gain"))) %>% 
  compare_means(formula = Exhausted ~ value, group.by = "Name", p.adjust.method = "BH") %>% 
  filter(p<0.05)

p6 <- df_gain %>% 
  ggplot(aes(x=value, y=Name, fill=MeanExhausted))+
  geom_tile()+
  scale_fill_gradientn(name="Expression", colours = brewer.pal(5, "PRGn"), limits=c(0,0.17))+
  geom_vline(xintercept = 1.5, color="black", size=0.25)+
  geom_segment(data=data.frame(y=seq(1.5, 31.5,1)), inherit.aes = F, 
               aes(y = y, yend = y, x=0.25, xend=2.8), color="white", size=0.25)+
  geom_text(data = p.values, inherit.aes = F, aes(y=Name, x=2.8, label=round(p, 2)), size=2.5)+
  coord_cartesian(clip = "off")+
  ggtitle("CN gain")+
  mytheme_1+
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        axis.text.x = element_text(angle=45, hjust=1, size=7),
        axis.ticks = element_blank(),
        axis.text.y = element_text(size=7, margin = unit(c(0,-0.15,0,0), units = "cm")),
        plot.margin = unit(units = "cm", c(0,1.25,0,0)),
        axis.title = element_blank(),
        panel.border = element_blank())+
  labs(tag = "F")

```

### Copy number loss
```{r copy number loss SF6}

df_loss <- df_snvs_chapuy %>% 
  filter(Description=="CN loss") %>% 
  mutate(value=factor(value, levels = c("0", "1", "2"), labels = c("wt", "loss", "loss"))) %>% 
  group_by(Name, value) %>%
  dplyr::summarise(MeanExhausted=mean(Exhausted/Absolute)) 

p.values <- df_snvs_chapuy %>% 
  filter(Description=="CN loss") %>% 
  mutate(Exhausted=Exhausted/Absolute) %>% 
  mutate(value=factor(value, levels = c("0", "1", "2"), labels = c("wt", "loss", "loss"))) %>% 
  compare_means(formula = Exhausted ~ value, group.by = "Name", p.adjust.method = "BH") %>% 
  filter(p<0.05) %>% 
  mutate(p=ifelse(p<0.005, 0.0051, p))

p7 <- df_loss %>% 
  ggplot(aes(x=value, y=Name, fill=MeanExhausted))+
  geom_tile()+
  scale_fill_gradientn(name="Expression", colours = brewer.pal(5, "PRGn"), limits=c(0,0.17))+
  geom_vline(xintercept = 1.5, color="black", size=0.25)+
   geom_segment(data=data.frame(y=seq(1.5, 32.5,1)), inherit.aes = F, 
               aes(y = y, yend = y, x=0.25, xend=2.8), color="white", size=0.25)+
  geom_text(data = p.values, inherit.aes = F, aes(y=Name, x=2.8, label=round(p, 2)), size=2.5)+
  ggtitle("CN loss")+
    coord_cartesian(clip = "off")+
  mytheme_1+
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        axis.text.x = element_text(angle=45, hjust=1, size=7),
        axis.ticks = element_blank(),
        axis.text.y = element_text(size=7, margin = unit(c(0,-0.15,0,0), units = "cm")),
        axis.title = element_blank(),
        panel.border = element_blank())+
  labs(tag = "G")

```

### Structural variants
```{r structural variants SF6}

df_struct <- df_snvs_chapuy %>% 
  filter(Description=="SV") %>% 
  mutate(value=factor(value, levels = c("0", "3"), labels = c("wt", "mut")))%>% 
   group_by(Name, value) %>%
  dplyr::summarise(MeanExhausted=mean(Exhausted/Absolute)) 

p.values <- df_snvs_chapuy %>% 
  filter(Description=="SV") %>% 
  mutate(Exhausted=Exhausted/Absolute) %>% 
   mutate(value=factor(value, levels = c("0", "3"), labels = c("wt", "mut"))) %>% 
  compare_means(formula = Exhausted ~ value, group.by = "Name", p.adjust.method = "BH") %>% 
  filter(p<0.05) 

p8 <- df_struct %>% 
  ggplot(aes(x=value, y=Name, fill=MeanExhausted))+
  geom_tile()+
  scale_fill_gradientn(name="Expression", colours = brewer.pal(5, "PRGn"), limits=c(0,0.17))+
  geom_vline(xintercept = 1.5, color="black", size=0.25)+
  geom_segment(data=data.frame(y=seq(1.5, 7.5,1)), inherit.aes = F, 
               aes(y = y, yend = y, x=0.25, xend=2.75), color="white", size=0.25)+
  ggtitle("Structural variants")+
  coord_cartesian(clip = "off")+
  mytheme_1+
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        axis.text.x = element_text(angle=45, hjust=1, size=7),
        axis.ticks = element_blank(),
        axis.text.y = element_text(size=7, margin = unit(c(0,-0.15,0,0), units = "cm")),
        plot.margin = unit(units = "cm", c(0,1.25,0,0)),
        axis.title = element_blank(),
        panel.border = element_blank())+
  labs(tag = "H")

```


### Part 2
```{r assemble plot SF6 part 2, fig.height=7}

p5+wrap_plots(p6/p8+plot_layout(heights = c(2,0.5)))+(p7/plot_spacer()+plot_layout(heights = c(2,0.5)))+
  plot_layout(nrow = 1, widths = c(2.9,1,1))

#ggsave(width = 18, height = 14, units = "cm", filename = "SF6_p2.pdf")

```

### Part 3 (Legend)
```{r assemble plot SF6 part 3, fig.height=1}

as_ggplot(get_legend(p8+guides(fill=guide_colorbar(nrow = 2, title = "Exhausted\nT-cells"))+
                       theme(legend.position = "right",
                             legend.key.height = unit(units = "cm", 0.3),
                             legend.key.width = unit(units = "cm", 0.3))))

ggsave(width = 2, height = 2.4, units = "cm", filename = "SF6_legend.pdf")

```

# Supplementary Figure 7
Panel A was generated using FlowJo. 

## Flow cytometry: IKZF3
```{r, fig.height=3}

med <- df_ikzf3 %>% filter(Entity=="rLN") %>% pull(`FoxP3+/IKZF3+`) %>% median()

pvalues <- df_ikzf3 %>% rename(IKZF3=`FoxP3+/IKZF3+`) %>% 
  data.frame() %>% 
  compare_means(data=., formula = IKZF3 ~ Entity, ref.group = "rLN") %>% 
  select(Entity=group2, p) %>% 
  filter(p<0.05) %>% 
  mutate(p=round(p,3))

nrow(df_ikzf3)

plot_aiolos <- df_ikzf3 %>% 
  ggplot(aes(x=Entity,y=`FoxP3+/IKZF3+`))+
  geom_hline(yintercept = med, size=0.25, linetype="dashed", color="grey60")+
  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 %>% mutate(Y=c(75, 75)),
            aes(x=Entity, y=Y, label=p), size=2.5, check_overlap = T)+
  scale_fill_brewer(palette = "Paired", limits=c("DLBCL", "MCL", "FL", "MZL", "rLN"))+
  scale_y_continuous(limits = c(0,80), name=expression('% IKZF3'^'+'~'of FoxP3'^'+'))+
  scale_x_discrete(limits=c("rLN", "DLBCL", "MCL", "FL", "MZL"))+
  ggtitle("Flow cytometry")+
  mytheme_1+
  theme(legend.position = "none",
        strip.background = element_rect(color=NA),
        axis.title.x = element_blank(),
        plot.title = element_text(face = "plain", size=7),
        panel.border = element_rect(size=0.5),
        axis.text.x = element_text(angle=45, hjust = 1, size=7),
        axis.text.y = element_text(size=7),
        axis.title.y = element_text(size=7),
        panel.background = element_rect(fill=NA),
        plot.margin = unit(c(0,0.25,0,0.25), "cm"))

plot_spacer()+plot_aiolos+plot_layout(widths = c(3,1))

ggsave(width = 19, height = 5.7, units = "cm", filename = "SF7.pdf")

```

## Treg clonotypes
```{r treg clonotypes SF7, fig.height=3}

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)

treg_shared <- list()
for(i in c(8,13,15)){

df_subset <- 
  df_clonotypes_shared %>% 
  add_entity() %>% 
  filter(IdentIb==i) 

treg_shared[[i]] <- ggplot()+
  geom_point_rast(data=DFtotal_5prime,
                  aes(x=refUMAP_1, y=refUMAP_2, fill=IdentI), size=0.25, 
                  alpha=ifelse(DFtotal_5prime$IdentI==i, 0.4, 0.04), 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.15, alpha=0.4)+
  scale_color_manual(values = colors_umap_cl, limits=factor(cluster_order),
                    labels=unlist(labels_cl), guide="none")+
  scale_fill_manual(values = colors_umap_cl, limits=factor(cluster_order), guide="none",
                    labels=unlist(labels_cl))+
  guides(fill=guide_legend(nrow = 7, byrow = F, override.aes = list(size=1.75, stroke=0, shape=21, alpha=1, color="white")))+
  coord_cartesian(clip = "off")+
  xlab("refUMAP-1")+
  ylab("refUMAP-2")+
  mytheme_1+
  theme(legend.position = "right",
        panel.border = element_rect(size=0.25),
        plot.title = element_textbox_simple(size=7, width = NULL, face = "plain",
                                            padding = margin(1.25, 0, 1, 0),
                                            lineheight = 1.25,
                                            halign=0.5),
        legend.text = element_text(size=7),
        legend.spacing.x = unit("cm", x = 0.13),
        axis.title.x = element_text(size=7),
        axis.title.y = element_text(size=7),
        axis.text = element_text(size=7),
        legend.spacing.y = unit("cm", x = 0.001),
        legend.key.width = unit("cm", x = 0.05),
        legend.key.height = unit("cm", x = 0.5),
        legend.box.margin = margin(unit = "cm",c(0,-0.35,0,-1)),
        legend.title = element_blank())

if(i==8)
  treg_shared[[i]] <- treg_shared[[i]]+
  labs(title="Paired clonotypes of <span style='color:#C6DBEF'>T<sub>REG</sub> CM<sub>1</sub></span>",
       tag = "C")
  
if(i==13)
 treg_shared[[i]] <- treg_shared[[i]]+
 labs(title="Paired clonotypes of <span style='color:#6BAED6'>T<sub>REG</sub> CM<sub>2</sub></span>")+
  theme(axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank())

if(i==15)
 treg_shared[[i]] <- treg_shared[[i]]+
 labs(title="Paired clonotypes of <span style='color:#2171B5'>T<sub>REG</sub> EM<sub>1</sub></span>")+
    theme(axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank())
  
}

plot_treg <- treg_shared[[8]]+treg_shared[[13]]+treg_shared[[15]]+plot_layout(guides = "collect")
#plot_treg
```

## Survival analysis
```{r survival plots}

kmplot_tfh <- readRDS("data/SurvPlot_Tfh_Gallium.rds")
kmplot_tfh$plot$plot_env$legend <- c(0.32, 0.2)
kmplot_tfh$plot$theme$legend.position <- c(0.32, 0.2)
kmplot_tfh$plot$theme$legend.background$fill <- NA
kmplot_tfh$plot$theme$legend.text$colour <- NA
kmplot_tfh$plot$theme$legend.text$size <- 7
kmplot_tfh$plot <- kmplot_tfh$plot+annotation_custom(grob = textGrob(label = expression('T'[FH]~'High'), gp = gpar(cex=0.5), x=0.36, y=0.215))+
  annotation_custom(grob = textGrob(label = expression('T'[FH]~'Low'), gp = gpar(cex=0.5), x=0.36, y=0.125))+
  labs(tag = "D")

kmplot_treg <- readRDS("data/SurvPlot_TregEff2_Gallium.rds")
kmplot_treg$plot$plot_env$legend <- c(0.32, 0.2)
kmplot_treg$plot$theme$legend.position <- c(0.32, 0.2)
kmplot_treg$plot$theme$legend.background$fill <- NA
kmplot_treg$plot$theme$legend.text$colour <- NA
kmplot_tfh$plot$theme$legend.text$size <- 7
kmplot_treg$plot <- kmplot_treg$plot+annotation_custom(grob = textGrob(label = expression('T'[REG]~'EM'[2]~'High'), gp = gpar(cex=0.5), x=0.36, y=0.215))+
  annotation_custom(grob = textGrob(label = expression('T'[REG]~'EM'[2]~'Low'), gp = gpar(cex=0.5), x=0.36, y=0.125))+
  labs(tag = "E")


```

## Assemble plot
```{r}

plot_treg/
  wrap_plots(kmplot_tfh$plot+kmplot_treg$plot+plot_spacer())+
  plot_layout(heights = c(1.2,1))

#ggsave(width = 18, height = 11, units = "cm", filename = "SF7.pdf")

```

# Supplementary Figure 8
## Dendrogram T-cells
```{r dendrogram}

# Create data frame
data <- data.frame(
  level1="_Tcells",
  level2=c("_'T'[Pr]",
           rep("_'T'[H]",2),  
           "_'T'[FH]", 
           rep("_'T'[REG]",1),  
           rep("_'T'[TOX]",3)),
  level3=c("_'T'[Pr]", 
           "TH_'CD4'^'+'*' Naive'",
           "TH_'non-Naive (CM'[1]*' + CM'[2]*')'", 
           "_'T'[FH]",  
           "_'T'[REG]",
           "TTOX_'CD8'^'+'*' Naive'",
           "TTOX_'non-Naive (EM'[1]*' + EM'[2]*')'", 
           "TTOX_'   Exhausted (EM'[3]*')'")
)

# Data handling
edges_level1_2 <- data %>% select(level1, level2) %>% unique %>% rename(from=level1, to=level2)
edges_level2_3 <- data %>% select(level2, level3) %>% unique %>% rename(from=level2, to=level3)
edge_list=rbind(edges_level1_2, edges_level2_3)
vert <- data.frame(
  name=unique(c(data$level1, data$level2, data$level3))) %>% 
  mutate(cluster=as.character(c(NA, 14, 'TH', 6, 'TREG', "TTOX", 1, 2, 12, 3, 5))) %>% 
  mutate(label=strsplit(name, split = "_") %>% sapply(., "[[", 2))

# Make ggraph object
mygraph_codex <- graph_from_data_frame( edge_list ,vertices = vert)

# Small codex dendrogramm
ggraph(mygraph_codex, layout = 'tree', circular = FALSE) + 
  geom_edge_diagonal(strength = 1.4, edge_width=0.25)+
  geom_node_label(aes(label=label), 
                  parse = T, nudge_y=0.11, label.padding =  unit(units = "cm", 0.12),
                  size=2.35, alpha=1, 
                  fill=c(rep("white", 4), "black", rep("white", 6)), vjust=1, color=NA,
                  label.size = 0, label.r = unit(units = "cm", 0))+
  geom_node_text(aes(label=label, color=cluster), 
                 vjust=1, nudge_y=0.05, 
                 parse = T, 
                 alpha=c(0,rep(1,5),rep(0,5)),
                 size=2.35)+
  scale_color_manual(values = colors_dendrogramm_codex)+
  coord_cartesian(clip = "off")+
  ggtitle("T-cell subsets \nidentified by mIF")+
  theme_void()+
  theme(legend.position = "none")+
  theme(plot.margin = unit(c(-0.5,0.35,0.25,0.25), units = "cm"),
        plot.title = element_text(hjust = 0.4, vjust=-7.5, size=7))

#ggsave(filename = "Figure8_p2.pdf", width = 12.6, height = 4.25, units = "cm")

```

## Handle data
```{r handle data SF8}
  
# Read CODEX expression data
# Available at BioStudies database (https://www.ebi.ac.uk/biostudies/) under accession number S-BIAD565
codex_expression <- data.table::fread("data/cells_expression.csv") %>% tibble() %>% 
  rename(unique_cell_id=V1)

proteins_selected <- c("PAX5", "CD20", "CD79a", "CD21", "PDPN", "CD38", "MCT",  "GRZB", "CD56", "CD163", "CD206", "CD11c",
                       "CD15", "CD34", "CD31", "CD90", "Ki67", "PD1", "CXCR5", "ICOS", "CD69", "CD45RO", "TIM3", "LAG3",
                       "CD57", "CD8", "CD45RA", "FOXP3", "CD4",  "CD3")

codex_meanExp <- codex_expression %>% 
  left_join(., codex_annotation %>% select(unique_cell_id, Merged_final)) %>% 
  filter(!is.na(Merged_final)) %>% 
  select(-unique_cell_id) %>% 
  group_by(Merged_final) %>% 
  summarise_all(mean) %>% 
  pivot_longer(cols = 2:ncol(.), names_to = "Protein", values_to = "Expression") %>% 
  group_by(Protein) %>% 
  mutate(Expression=(Expression-min(Expression))/(max(Expression)-min(Expression))) %>% 
  filter(Protein %in% proteins_selected)

```

## Makers and cell types
```{r plot markers SF8}

p2 <- codex_meanExp %>% 
  ggplot(aes(x=Protein, y=Merged_final, fill=Expression))+
  geom_tile()+
  scale_fill_gradientn(name="Expression", colours = brewer.pal(5, "GnBu"), limits=c(0,1), breaks=c(0,0.5,1))+
  geom_hline(yintercept = seq(1.5, 17.5, 1), size=0.25, color="white")+
  geom_vline(xintercept = seq(1.5, 32.5, 1), size=0.25, color="white")+
  scale_x_discrete( expand = c(0,0), limits=c(proteins_selected))+
  scale_y_discrete(expand = c(0,0), limits=c("B", "FDC", "PC", "MC", "NK", "NKT",  "Macro", "DC", "Granulo", "Stromal cells",
                                             "TPR", "TFH",  "TTOX_exh", "TTOX", "TTOXNaive",   "Treg",  "CD4T", "CD4TNaive"),
                   labels=unlist(list("B-cells", "FDC", "Plasma cells", "Mast cells", "NK cells", "NK T-cells", "Macrophages", "Dendritic cells", "Granulocytes",
                                      "Stromal cells", labels_codex$TPR, labels_codex$TFH, labels_codex$TTOX_exh, labels_codex$TTOX, labels_codex$TTOXNaive, 
                                      labels_codex$Treg, labels_codex$CD4T, labels_codex$CD4TNaive)))+
  guides(fill=guide_colorbar(ticks.colour = "black"))+
  theme_bw()+
  mytheme_1+
   theme(axis.title = element_blank(),
         legend.position = "top",
         axis.text.y = element_text(size=7),
         legend.text = element_text(size = 7, color="black"),
         legend.title = element_text(size = 7, color="black", vjust = 0.8, margin = unit(units = "cm", c(0,0.2,0,0))),
         legend.key.height = unit(0.25, "cm"),
         plot.margin = unit(c(0,0,0,0), units = "cm"),
         legend.key.width = unit(0.2, "cm"),
         legend.box.spacing = unit(0.1, "cm"),
         legend.box.margin = unit(c(0,0,0,0), units = "cm"),
         plot.title = element_text(face = "plain", vjust = -1),
         plot.tag = element_text(margin = unit(c(0,-0.5,-0.25,0), units = "cm")),
         axis.text.x = element_text(size=6.5, angle = 45, hjust = 1))+
  labs(tag = "C")

```

## T-cell numbers in codex
```{r T-cell numbers SF8}

df_codex_no <- 
  codex_annotation %>% 
  filter(Merged_final %in% c("TTOXNaive", "TTOX_exh", "TTOX", "Treg", "TPR", "TFH", "CD4T", "CDT4Naive")) %>% 
  count(PatientID, Merged_final, unique_region) %>% 
  group_by(PatientID) %>% 
  mutate(Sum=sum(n)) %>% 
  ungroup() %>% 
  mutate(No=dense_rank(desc(Sum)))

regions_random <- df_codex_no %>% 
  select(PatientID, unique_region) %>% 
  distinct() %>% 
  group_by(PatientID) %>% 
  sample_n(1)

p3 <- ggplot()+
  geom_hline(yintercept = 47280, size=0.25, linetype="dashed")+
  geom_bar(data=df_codex_no %>% filter(unique_region %in% regions_random$unique_region), 
           aes(x=No-0.15, y=n, fill=Merged_final), color="white",
           stat = "identity",  width=0.2, size=0.25, alpha=0.7)+
  geom_bar(data=df_codex_no %>% filter(!unique_region %in% regions_random$unique_region), 
           aes(x=No+0.15, y=n, fill=Merged_final), color="white",
           stat = "identity",  width=0.25, size=0.25, alpha=0.7)+
  scale_y_continuous(name = "Absolute number of cells - mIF", limits = c(0, 100000))+
 scale_fill_manual(values = colors_codex[c(2:9)],
                   limits=limits_codex[c(2:9)],
                   labels=labels_codex[c(2:9)],
                   name=NULL)+
  guides(fill=guide_legend(nrow = 2, default.unit = "cm", override.aes = list(color="white"), 
                           keywidth = 0.3, keyheight = 0.3, byrow = T))+
  scale_x_continuous(name="Patients", 
                     breaks=unique(df_codex_no$No), 
                     labels=unique(df_codex_no$PatientID), 
                     expand = c(0.02,0.02))+
  mytheme_1+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.margin = unit(c(0,0,0,0), units = "cm"),
        legend.box.margin = unit(c(0,0,-0.3,0), units = "cm"),
        axis.title.y = element_text(size=7, vjust = -15),
        axis.title.x = element_blank(),
        legend.spacing.y = unit(0.15, units = "cm"),
        plot.tag = element_text(margin = unit(c(0,-0.5,0,0), units = "cm")),
        legend.position = "top")+
  labs(tag = "D")

```

## Correlations
```{r correlations SF8}

freq_codex <- 
  codex_annotation %>% 
  filter(Merged_final %in% c("TTOXNaive", "TTOX_exh", "TTOX", "Treg", "TPR", "TFH", "CD4T", "CD4TNaive")) %>% 
  filter(!unique_region %in% c("191_1reg004", "191_4reg004", "191_4reg005", "191_1reg006")) %>% 
  add_prop(vars = c("Merged_final", "PatientID"), group.vars = 2) %>% 
  rename(Prop_codex=Prop, IdentI=Merged_final) 

freq_citeseq <- 
  Combined_T@meta.data %>% 
  add_prop(vars = c("IdentI", "PatientID"), group.vars = 2) %>% 
  mutate(IdentII=case_when(IdentI==1 ~ "CD4TNaive",
                           IdentI %in% c(2,9) ~ "CD4T",
                           IdentI==14 ~ "TPR",
                           IdentI==6 ~ "TFH",
                           IdentI %in% c(8,11,13,15) ~ "Treg",
                           IdentI==12 ~ "TTOXNaive",
                           IdentI %in% c(3,16) ~ "TTOX",
                           IdentI %in% c(5) ~ "TTOX_exh")) %>% 
  group_by(IdentII, PatientID) %>% 
  summarise(Prop=sum(Prop)) %>% 
  rename(Prop_citeseq=Prop, IdentI=IdentII) %>% 
  fill_zeros(names_from = "IdentI", values_from = "Prop_citeseq")

freq_joined <- left_join(freq_codex, freq_citeseq) %>% 
  mutate(Prop_codex=100*Prop_codex, Prop_citeseq=100*Prop_citeseq)

this_theme <- 
  theme_bw()+
  mytheme_1+
  theme(plot.margin = unit(c(0,0.2,0.2,0.35), units = "cm"),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        plot.title = element_textbox_simple(face = "plain", halign=0.5, width = 2, padding = margin(0, 0, 0, 0)),
        plot.tag = element_text(margin = unit(c(0,-0.5,0,0), units = "cm")),
        axis.text = element_text(size=7, color="black"))

cor_plots_codex <- list()

cor_plots_codex[["TFH"]] <- 
  freq_joined %>% 
  filter(IdentI=="TFH") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["TFH"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_codex[["TFH"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.5), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,40), breaks = c(0, 12, 24, 36), name = "mIF")+
  scale_y_continuous(limits = c(0,40), breaks = c(0, 12, 24, 36), name = "CITE-seq")+
  labs(title="T<sub>FH</sub>")+
  coord_fixed()+
  this_theme+
  theme(plot.margin = unit(c(0,0,0.2,0), units = "cm"))

cor_plots_codex[["TREG"]] <- 
  freq_joined %>% 
  filter(IdentI=="Treg") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["Treg"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color="black", se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.5), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,56), breaks = c(0, 16, 32, 48), name = "mIF")+
  scale_y_continuous(limits = c(0,56), breaks = c(0, 16, 32, 48), name = "CITE-seq")+
  labs(title="T<sub>REG</sub>")+
  coord_fixed()+
  this_theme+
  theme(axis.title.y = element_text(size=7, angle = 90, vjust = 2.5))

cor_plots_codex[["TTOXNaive"]] <- 
  freq_joined %>% 
  filter(IdentI=="TTOXNaive") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["TTOXNaive"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_codex[["TTOXNaive"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.5), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,13), breaks = c(0, 4, 8, 12), name = "mIF")+
  scale_y_continuous(limits = c(0,13), breaks = c(0, 4, 8, 12), name = "CITE-seq")+
  labs(title="CD8<sup>+</sup> Naive")+
  coord_fixed()+
  this_theme+
  theme(plot.margin = unit(c(0,0,0.2,0), units = "cm"))

cor_plots_codex[["THNaive"]] <- 
  freq_joined %>% 
  filter(IdentI=="CD4TNaive") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["CD4TNaive"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_codex[["CD4TNaive"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.53), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,30), breaks = c(0, 10, 20, 30), name = "mIF")+
  scale_y_continuous(limits = c(0,30), breaks = c(0, 10, 20, 30), name = "CITE-seq")+
  labs(title="CD4<sup>+</sup> Naive")+
  coord_fixed()+
  this_theme+
  theme(plot.title = element_textbox_simple(face = "plain", halign=0.5, margin = unit(units = "cm", c(0,0,-1.75,0)), 
                                            width = 2, padding = margin(0, 0, 0, 0)),
        plot.margin = unit(c(0,0,0.2,0), units = "cm"))

cor_plots_codex[["TPR"]] <- 
  freq_joined %>% 
  filter(IdentI=="TPR") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["TPR"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_codex[["TPR"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.5), label.y.npc = c(0.93))+
  scale_x_continuous(limits = c(0,13), breaks = c(0, 4, 8, 12), name = "mIF")+
  scale_y_continuous(limits = c(0,13), breaks = c(0, 4, 8, 12), name = "CITE-seq")+
  labs(title="T<sub>Pr</sub>")+
  coord_fixed()+
  this_theme+
  theme(axis.title.y = element_text(size=7, angle = 90, vjust = 2.5),
        plot.title = element_textbox_simple(face = "plain", halign=0.5, margin = unit(units = "cm", c(0,0,-1.75,0)), 
                                            width = 2, padding = margin(0, 0, 0, 0)),
        plot.tag = element_text(margin = unit(c(0,0,-0.25,0), units = "cm")))

cor_plots_codex[["TH"]] <- 
  freq_joined %>% 
  filter(IdentI=="CD4T") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["CD4T"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_codex[["CD4T"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.5), label.y.npc = c(0.93))+
  scale_x_continuous(limits = c(0,60), breaks = c(0, 20, 40, 60), name = "mIF")+
  scale_y_continuous(limits = c(0,60), breaks = c(0, 20, 40, 60), name = "CITE-seq")+
  labs(title="Memory T<sub>H</sub>")+
  coord_fixed()+
  this_theme+
  theme(axis.title.y = element_text(size=7, angle = 90, vjust = 2.5))

cor_plots_codex[["TTOX"]] <- 
  freq_joined %>% 
  filter(IdentI=="TTOX") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["TTOX"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_codex[["TTOX"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.5), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,45), breaks = c(0, 15, 30, 45), name = "mIF")+
  scale_y_continuous(limits = c(0,45), breaks = c(0, 15, 30, 45), name = "CITE-seq")+
  labs(title="Memory T<sub>TOX</sub>")+
  coord_fixed()+
  this_theme+
  theme(axis.title.y = element_text(size=7, angle = 90, vjust = 2.5),
        axis.title.x = element_text(size=7))

cor_plots_codex[["TTOX_exh"]] <- 
  freq_joined %>% 
  filter(IdentI=="TTOX_exh") %>% 
  ggplot(aes(x=Prop_codex, y=Prop_citeseq))+
  geom_point(fill=colors_codex[["TTOX_exh"]], size=1, alpha=0.75, shape=21, stroke=0.1)+
  geom_smooth(method = "lm", linetype="dashed", size=0.25, formula = y ~ x, na.rm = T,
              color=colors_codex[["TTOX_exh"]], se=F, fullrange=T)+
  stat_cor(aes(label=..r.label..), method = "pearson", size=2.5, label.x.npc = c(0.5), label.y.npc = c(0.1))+
  scale_x_continuous(limits = c(0,65), breaks = c(0, 20, 40, 60), name = "mIF")+
  scale_y_continuous(limits = c(0,65), breaks = c(0, 20, 40, 60), name = "CITE-seq")+
  labs(title="PD1<sup>+</sup> TIM3<sup>+</sup> T<sub>TOX</sub>")+
  coord_fixed()+
  this_theme+
  theme(axis.title.x = element_text(size=7),
        plot.margin = unit(c(0,0,0.2,0), units = "cm"))

freq_joined %>% 
  group_by(IdentI) %>% 
  summarise(R=cor.test(Prop_codex, Prop_citeseq)$estimate) %>% pull(R) %>% median()

```

## Assemble plot
```{r assemble SF9, fig.height=5.5}

p_full <- wrap_plots(p2/p3+plot_layout(heights = c(1.25,1)))+wrap_plots(cor_plots_codex$TPR+labs(tag = "E")+cor_plots_codex$THNaive+cor_plots_codex$TH+cor_plots_codex$TFH+
  cor_plots_codex$TREG+cor_plots_codex$TTOXNaive+cor_plots_codex$TTOX+cor_plots_codex$TTOX_exh+
  plot_layout(ncol = 2))+
  plot_layout(widths = c(2.25,1.15))
p_full
#ggsave(p_full, width = 18, height = 14, units = "cm", filename = "SF8.pdf")


```

# Supplementary Figure 9
## Immunofluorescence images
```{r images SF9, fig.height=3}

plots_codex <- list()

for(r in c("191_3reg008", "191_4reg004", "191_2reg007", "191_5reg002", "191_1reg003", "empty")) {
  
  df_tmp <- codex_annotation %>% filter(unique_region== r) %>% 
    mutate(Merged_all_simple=ifelse(Merged_final %in% c("Granulo", "Macro", "DC"), "Myeloid", Merged_final)) %>% 
    mutate(Merged_all_simple=ifelse(Merged_all_simple %in% c("MC", "NKT", "PC", "NK"), "Other", Merged_all_simple)) %>% 
    filter(((x-mean(.$x))^2+(y-mean(.$y))^2)<2500^2) 
  
  plots_codex[[r]] <- ggplot()+
    geom_point_rast(data=df_tmp %>% filter(Merged_all_simple=="B"), aes(x=x,y=y), shape=21, size=0.25, stroke=0, alpha=1, raster.dpi =300,
                    color=colors_codex[["B"]], fill=colors_codex[["B"]])+
    geom_point_rast(data=df_tmp %>% filter(Merged_all_simple!="B"), aes(x=x,y=y, fill=Merged_all_simple, color=Merged_all_simple), 
                    shape=21, size=0.25, stroke=0, alpha=1, raster.dpi=300)+
    scale_color_manual(values = colors_codex, limits=limits_codex, labels=labels_codex, name="Cell type")+
    scale_fill_manual(values = colors_codex, limits=limits_codex, labels=labels_codex, name="Cell type")+
    ggtitle(unique(df_tmp$Entity))+
    coord_fixed()+
    theme_void()+
    theme(legend.position = "none",
          plot.title = element_text(color="white", hjust=0.1, size=8, 
                                    margin = unit(units = "cm", c(0,0,-0.6,0)), face = "bold"),
          plot.margin = unit(units = "cm", c(0.1, 0.1, 0.1, 0.1)),
          panel.background = element_rect(fill = "black", color="black"))
}

plots_codex

#ggsave(wrap_plots(plots_codex), width = 19, height = 12.5, units = "cm", filename = "SFigure9_p2.pdf")

legend_plot_codex <- ggplot()+
  geom_point_rast(data=df_tmp %>% filter(Merged_all_simple!="B"), aes(x=x,y=y, fill=Merged_all_simple, color=Merged_all_simple), 
                  shape=21, size=0.25, stroke=0, alpha=0, raster.dpi=300)+
  scale_color_manual(values = colors_codex, limits=limits_codex, labels=labels_codex, name="Cell type")+
  scale_fill_manual(values = colors_codex, limits=limits_codex, labels=labels_codex, name="Cell type")+
  guides(fill=guide_legend(ncol = 2, override.aes = list(size=1.75, stroke=0, shape=21, alpha=1, color=NA)))+
  guides(color=guide_legend(ncol = 2))+
  ggtitle("")+
  coord_fixed()+
  mytheme_codex+
  theme(panel.background = element_rect(fill = "black", color="black"),
        legend.position = "right",
        legend.text = element_text(size=6, color="white"),
        legend.box.background = element_rect(fill = "black", color="black"),
        legend.spacing.x = unit("cm", x = 0.13),
        legend.spacing.y = unit("cm", x = 0.001),
        legend.key.width = unit("cm", x = 0.05),
        legend.key.height = unit("cm", x = 0.5),
        legend.title = element_blank())

as_ggplot(get_legend(legend_plot_codex))
#ggsave(width = 4, height = 4, units = "cm", filename = "SFigure9_p2_legend.pdf")

```

# Supplementary Figure 10
## Load analysis
```{r run analysis SF10}

# Read results from neighborhood analysis
# Please run file: analysis/NeighborhoodAnalysis.Rmd
load("output/Neighborhood_results.RData")

# Add codex annotation
codex_annotation <- left_join(codex_annotation, nn_classes, by="unique_cell_id")
codex_annotation

```

## Tissue cores
### Images
```{r plots SF10, fig.height=6}

regions <- codex_annotation %>% pull(unique_region) %>% unique()
plots <- list()
df <- list()
for(r in regions){
  
df[[r]] <- codex_annotation %>% 
  filter(!is.na(Region), unique_region %in% r) %>% 
  filter(x>500, x<7500) %>% 
  filter(y>500, y<7500) %>% 
  filter(((x-mean(.$x))^2+(y-mean(.$y))^2)<2500^2)
  
plots[[r]] <- df[[r]] %>% 
  ggplot()+
  ggrastr::geom_point_rast(aes(x=x,y=y,color=Region, fill=Region), shape=21, size=0.25, stroke=0, alpha=1, raster.dpi =400)+
  scale_color_manual(values = colors_nn)+
  scale_fill_manual(values = colors_nn)+
  guides(color=guide_legend(override.aes = list(size=3)))+
  ggtitle(unique(df[[r]]$PatientID))+
  coord_fixed(clip = "off")+
  theme_void()+
  theme(legend.position = "none",
        plot.margin = unit(units = "cm", c(0.1,0.1,0.1,0.1)),
        plot.subtitle = element_text(size=7, face = "bold", hjust=0.5, margin = unit(units = "cm", c(0,0,0,0))),
        plot.title = element_text(size=6.5, face = "plain", margin = unit(units = "cm", c(-0.1,0,-0.6,-0.1))),
        panel.background = element_rect(fill = NA, color = NA),
        plot.background = element_rect(fill = NA, color = NA))

}

empty <- codex_annotation %>% 
  filter(!is.na(Region), unique_region %in% "empty") %>% 
  ggplot()+
  ggrastr::geom_point_rast(aes(x=x,y=y,color=Region, fill=Region), alpha=0, raster.dpi =400)+
  guides(color=guide_legend(override.aes = list(size=3, alpha=1)))+
  ggtitle("")+
  coord_fixed(clip = "off")+
  theme_void()+
  theme(legend.position = "right",
        plot.margin = unit(units = "cm", c(0.1,0.1,0.1,0.1)),
        plot.title = element_text(size=8, face = "plain", margin = unit(units = "cm", c(0,0,-0.75,0))),
        panel.background = element_rect(fill = NA, color = NA),
        plot.background = element_rect(fill = NA, color = NA))


p_full <- 
  wrap_plots(plots$`191_1reg006`+labs(tag = "A", subtitle = "rLN")+
               plots$`191_3reg007`+
               plots$`191_5reg005`+plot_layout(ncol = 1))+
  
  wrap_plots(plots$`191_4reg004`+labs(tag = "B",  subtitle = "DLBCL")+#
               plots$`191_3reg001`+
               plots$`191_4reg006`+plot_layout(ncol = 1))+
  
  wrap_plots(plots$`191_2reg007`+labs(tag = "C",  subtitle = "MCL")+
               plots$`191_2reg002`+
               plots$`191_3reg006`+plot_layout(ncol = 1))+
  
  wrap_plots(plots$`191_5reg002`+labs(tag = "D",  subtitle = "FL")+
               plots$`191_3reg002`+
               plots$`191_5reg001`+plot_layout(ncol = 1))+
  
  wrap_plots(plots$`191_1reg003`+labs(tag = "E",  subtitle = "MZL")+
               plots$`191_1reg004`+
               empty+plot_layout(ncol = 1))+
  plot_layout(ncol = 5)

p_full

#ggsave(p_full, width = 18, height = 11.5, units = "cm", filename = "SF11.pdf")

```

### Legend
```{r legend SF10, fig.height=1}

labels_nn <- c(
  "N1: B-cells / FDC" ,
  'N2: B-cells / FDC / T'[FH]~'',
  'N3: T'[Pr]~'/ T'[REG]~'',
  'N4: Macrophages / B-cells / Exh. T'[TOX]~'',
  "N5: B-cells",
  "N6: T-cell area I" ,
  "N7: T-cell area II" ,
  'N8: PC / NK / Memory T'[TOX]~'' ,
  "N9: T-cell area III" ,
  "N10: Stromal cells / Macrophages"
   )

p_legend <- codex_annotation %>% 
  filter(!is.na(Region), unique_region %in% r) %>% 
  ggplot()+
  ggrastr::geom_point_rast(aes(x=x,y=y,color=Region, fill=Region), shape=21, size=0.25, 
                           stroke=0, alpha=1, raster.dpi =400)+
  scale_color_manual(values = colors_nn, limits=factor(1:10), labels=labels_nn)+
  scale_fill_manual(values = colors_nn, limits=factor(1:10), labels=labels_nn)+
  guides(color=guide_legend(nrow = 2, override.aes = list(size=2, color="black", stroke=0.25)))+
  ggtitle(unique(df[[r]]$PatientID))+
  coord_fixed(clip = "off")+
  theme_void()+
  theme(legend.position = "right",
        plot.margin = unit(units = "cm", c(0.1,0.1,0.1,0.1)),
        legend.text = element_text(size=6.5),
        legend.title = element_text(size=7, face = "bold"),
        plot.subtitle = element_text(size=7, face = "bold", hjust=0.5, margin = unit(units = "cm", c(0,0,0,0))),
        plot.title = element_text(size=6.5, face = "plain", margin = unit(units = "cm", c(-0.1,0,-0.6,-0.1))),
        panel.background = element_rect(fill = NA, color = NA),
        plot.background = element_rect(fill = NA, color = NA),
        legend.box.margin = unit(c(0,0,0,-0.38), "cm"),
        legend.key.width = unit("cm", x = 0.1),
        legend.spacing.x = unit("cm", x = 0.1),
        legend.key.height = unit("cm", x = 0.35))

as_ggplot(get_legend(p_legend))

#ggsave(width = 19, height = 1.5, units = "cm", filename = "SF11_legend.pdf")

```

## Composition of neighborhoods
```{r neighborhood composition, fig.height=3}

df_freq_nh <- codex_annotation %>% 
  add_prop(vars = c("Entity", "Region", "unique_region"), group.vars = 3) %>% 
  fill_zeros(names_from = "Region", values_from = "Prop") %>%
  group_by(Entity, Region) %>% 
  mutate(Max=0.06+max(Prop),
         Region_label=paste0("N", Region)) %>% 
  mutate(Region_label=factor(Region_label, levels = paste0("N", 1:10))) %>% 
  mutate(Max=ifelse(Region_label=="N8" & Entity=="MCL", 0.18, Max)) %>% 
  ungroup()

pvalues <- df_freq_nh %>% 
  compare_means(data=., formula = Prop ~ Entity, ref.group = "rLN", 
                group.by = "Region_label", p.adjust.method = "BH") %>% 
  filter(p.adj<0.05) %>% 
  mutate(p.adj_s=format(p.adj, scientific = TRUE, digits=1)) %>% 
  mutate(p.adj_f=case_when(p.adj > 0.01 ~ as.character(round(p.adj, 2)),
                           p.adj==0.01 ~ "0.01",
                           p.adj < 0.01 ~ p.adj_s),
         Entity=group2) %>% 
  mutate(Entity=factor(Entity, levels = c("rLN", "DLBCL", "MCL", "FL", "MZL"))) %>% 
  left_join(., df_freq_nh %>% select(Region_label, Max, Entity) %>% distinct, by = c("Region_label", "Entity"))

df_medianLines <- df_freq_nh %>% 
  filter(Entity=="rLN") %>% 
  group_by(Region_label) %>% 
  dplyr::summarise(MedianProp=median(Prop)) 

df_freq_nh %>% 
  ggplot(aes(x=Entity, y=Prop)) +
  geom_hline(data=df_medianLines, aes(yintercept=MedianProp),
             size=0.25, linetype="dashed", color="grey60")+
  geom_boxplot(width=0.5, outlier.alpha = 0, size=0.25)+
  ggbeeswarm::geom_beeswarm(size=1, shape=21, stroke=0.1, cex = 2, aes(fill=Region))+
  geom_text(data=pvalues, inherit.aes = F, aes(y=Max, x=Entity, label=p.adj_f), size=2.5)+
  scale_fill_manual(values = colors_nn)+
  scale_x_discrete(limits=c("rLN", "DLBCL", "MCL", "FL", "MZL"))+
  facet_wrap(~Region_label, strip.position = "top", scales = "free_y", nrow = 2)+
  scale_y_continuous(name = "% of total area", expand = c(0,0.075))+
  mytheme_1+
  theme(strip.text.y = element_text(angle = 0, size=6),
        axis.text.x = element_text(angle=45, hjust=1),
        axis.title.x = element_blank(),
        strip.background = element_blank(),
        plot.margin = unit(c(0,0.1,0,0.1), "cm"))+
  labs(tag = "F")

ggsave(width = 18, height = 7, units = "cm", filename = "SF10.pdf")

```

# Session info
```{r session}

sessionInfo()

```
