---

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