library(Seurat)
library(openxlsx)
library(readxl)
library(tidyr)
library("RColorBrewer")
library(dplyr)
library(scRepertoire)
library(reshape2)
library(fgsea)
library(ggalluvial)
library("gridExtra")
library(matrixStats)
library(fmsb)

cluster_dir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A019/Analysis_MM"

wdir <- paste0(cluster_dir, "/Result3_MNOT")
ddir <- paste0(cluster_dir, "/results3/Liver_Tumor_TnK/")
pdir <- paste0(cluster_dir, "/Result3_MNOT/Plots")

# Custom functions
source(file = "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A018/Analysis_MM/last_scripts_MNOT/TAD018/plot_stellare_function.R")

# Load RDS
obs1 <- readRDS(paste0(ddir, "Liver_Tumor_TnK_final.rds"))
#obs1 <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A019/Analysis_MM/results3/Liver_TandNK_NO_TCR/Liver_TandNK_final.rds")

########################
#####Annotations#######
########################

obs1@meta.data$Tissue <- obs1@meta.data$PROP.Source
obs1@meta.data$RNA_Group <- obs1@meta.data$PROP.Condition

obs1@meta.data$RNA_Group[obs1@meta.data$RNA_Group == "TA_only"] <- "TA33"
obs1@meta.data$RNA_Group[obs1@meta.data$RNA_Group == "TA_combo"] <- "TA33.Combo"

obs1@meta.data$RNA_Group <- factor(x = obs1@meta.data$RNA_Group, levels = c("TA33","TA33.Combo"))

# Tissue Sample assign
obs1$Tube <- paste0(obs1$hash.ID,obs1$RNA_Group)
obs1@meta.data$Orig<-paste0(obs1@meta.data$Tissue, ".", obs1@meta.data$RNA_Group)
table(obs1@meta.data$Tissue,obs1@meta.data$RNA_Group)
table(obs1@meta.data$Orig)

#################create labels for all CD8 and CD4 T cells#############
obs1@meta.data$NewlablesG <- obs1@meta.data$NewlablesM
obs1@meta.data$NewlablesG[grep("CD4",obs1@meta.data$NewlablesM) ] <- "CD4 T cells"
obs1@meta.data$NewlablesG[grep("CD8",obs1@meta.data$NewlablesM)] <- "CD8 T cells"
unique(obs1@meta.data$NewlablesM)

################create labels CD4 and CD8 clusters without naive####################
obs1@meta.data$NewlablesS <- obs1@meta.data$NewlablesM
indCd4<- grepl("CD4",obs1@meta.data$NewlablesM) & !grepl("Naive",obs1@meta.data$NewlablesM)
indCd8<- grepl("CD8",obs1@meta.data$NewlablesM) & !grepl("Naive",obs1@meta.data$NewlablesM)
obs1@meta.data$NewlablesS[indCd4] <- "CD4 T cells"
obs1@meta.data$NewlablesS[indCd8] <- "CD8 T cells"
obs1@meta.data$NewlablesST <- paste0(obs1@meta.data$Tissue,".",obs1@meta.data$NewlablesS)


################################################
###Calculating hyperexpanded T cells#############
################################################
obs1$aa3.Tube <-  paste0(obs1$cdr3_aa2,".",obs1$Tube)
obs1$aa3.Tube[is.na(obs1$cdr3_aa2)] <- NA
obs1@meta.data$ClonoFreq <- obs1@meta.data  %>% group_by(aa3.Tube)  %>% mutate(count = n()) %>% ungroup() %>%
  as.data.frame() %>% select (count)
obs1@meta.data$ClonoFreq[is.na(obs1$cdr3_aa2),] <- NA
obs1@meta.data$Clonetype  <- NA
obs1@meta.data$Clonetype[obs1@meta.data$ClonoFreq==1] <- "unique"
obs1@meta.data$Clonetype[obs1@meta.data$ClonoFreq > 1 & obs1@meta.data$ClonoFreq <= 5] <- "small"
obs1@meta.data$Clonetype[obs1@meta.data$ClonoFreq > 5 & obs1@meta.data$ClonoFreq <= 30] <- "large"
obs1@meta.data$Clonetype[obs1@meta.data$ClonoFreq > 30] <- "hyperexpanded"
obs1$Clonetype <- factor(x = obs1$Clonetype, levels = c("unique", "small", "large", "hyperexpanded")) #order clonotype variable

obs1$Expansion <- "unique"
obs1$Expansion[obs1@meta.data$ClonoFreq>5] <- "expanded"

table(obs1@meta.data$Clonetype)

######TISSUE SHARED CLONOTYPES########
# Identify cells with the same Cdr3aa sequence in both tissues
df1<-data.frame(obs1@meta.data$cdr3_aa2,obs1@meta.data$hash.ID,obs1@meta.data$RNA_Group, obs1@meta.data$Tissue)
df2<-na.omit(df1)
df2$cdr3_mouse<-paste0(df2$obs1.meta.data.cdr3_aa2,".",df2$obs1.meta.data.hash.ID,".",df2$obs1.meta.data.RNA_Group)
shared_cdr3aa <- intersect(df2$cdr3_mouse[df2$obs1.meta.data.Tissue == "Liver"], df2$cdr3_mouse[df2$obs1.meta.data.Tissue == "Tumor"])
obs1$cdr3_mouse<-paste0(obs1$cdr3_aa2,".",obs1$hash.ID,".",obs1$RNA_Group)
obs1$Shared <- ifelse(obs1$cdr3_mouse %in% shared_cdr3aa, "Shared", "Unique")

# Create variables share tissue 
obs1@meta.data$Sharedtissue<-paste0(obs1@meta.data$Tissue,".",obs1@meta.data$Shared)
obs1@meta.data$Shared.NewlablesS<-paste0(obs1@meta.data$Shared,".",obs1@meta.data$NewlablesS)
obs1@meta.data$Shared.NewlablesG<-paste0(obs1@meta.data$Shared,".",obs1@meta.data$NewlablesG)
obs1@meta.data$Expansion.NewlablesG<-paste0(obs1@meta.data$Shared,".",obs1@meta.data$NewlablesG,".",obs1@meta.data$Expansion)

################################################


Tobs <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A019/Analysis_MM/results3/Tumor_TandNK/Tumor_TandNK_final.rds")


########################
#####Annotations#######
########################

Tobs@meta.data$Tissue <- Tobs@meta.data$PROP.Source
Tobs@meta.data$RNA_Group <- Tobs@meta.data$PROP.Condition

Tobs@meta.data$RNA_Group[Tobs@meta.data$RNA_Group == "TA_only"] <- "TA33"
Tobs@meta.data$RNA_Group[Tobs@meta.data$RNA_Group == "TA_combo"] <- "TA33.Combo"

Tobs@meta.data$RNA_Group <- factor(x = Tobs@meta.data$RNA_Group, levels = c("TA33","TA33.Combo"))

# Tissue Sample assign
Tobs$Tube <- paste0(Tobs$hash.ID,Tobs$RNA_Group)
Tobs@meta.data$Orig<-paste0(Tobs@meta.data$Tissue, ".", Tobs@meta.data$RNA_Group)
table(Tobs@meta.data$Tissue,Tobs@meta.data$RNA_Group)
table(Tobs@meta.data$Orig)

#################create labels for all CD8 and CD4 T cells#############
Tobs@meta.data$NewlablesG <- Tobs@meta.data$NewlablesM
Tobs@meta.data$NewlablesG[grep("CD4",Tobs@meta.data$NewlablesM) ] <- "CD4 T cells"
Tobs@meta.data$NewlablesG[grep("CD8",Tobs@meta.data$NewlablesM)] <- "CD8 T cells"
unique(Tobs@meta.data$NewlablesM)

################create labels CD4 and CD8 clusters without naive####################
Tobs@meta.data$NewlablesS <- Tobs@meta.data$NewlablesM
indCd4<- grepl("CD4",Tobs@meta.data$NewlablesM) & !grepl("Naive",Tobs@meta.data$NewlablesM)
indCd8<- grepl("CD8",Tobs@meta.data$NewlablesM) & !grepl("Naive",Tobs@meta.data$NewlablesM)
Tobs@meta.data$NewlablesS[indCd4] <- "CD4 T cells"
Tobs@meta.data$NewlablesS[indCd8] <- "CD8 T cells"
Tobs@meta.data$NewlablesST <- paste0(Tobs@meta.data$Tissue,".",Tobs@meta.data$NewlablesS)


################################################
###Calculating hyperexpanded T cells#############
################################################
Tobs$aa3.Tube <-  paste0(Tobs$cdr3_aa2,".",Tobs$Tube)
Tobs$aa3.Tube[is.na(Tobs$cdr3_aa2)] <- NA
Tobs@meta.data$ClonoFreq <- Tobs@meta.data  %>% group_by(aa3.Tube)  %>% mutate(count = n()) %>% ungroup() %>%
  as.data.frame() %>% select (count)
Tobs@meta.data$ClonoFreq[is.na(Tobs$cdr3_aa2),] <- NA
Tobs@meta.data$Clonetype  <- NA
Tobs@meta.data$Clonetype[Tobs@meta.data$ClonoFreq==1] <- "unique"
Tobs@meta.data$Clonetype[Tobs@meta.data$ClonoFreq > 1 & Tobs@meta.data$ClonoFreq <= 5] <- "small"
Tobs@meta.data$Clonetype[Tobs@meta.data$ClonoFreq > 5 & Tobs@meta.data$ClonoFreq <= 30] <- "large"
Tobs@meta.data$Clonetype[Tobs@meta.data$ClonoFreq > 30] <- "hyperexpanded"
Tobs$Clonetype <- factor(x = Tobs$Clonetype, levels = c("unique", "small", "large", "hyperexpanded")) #order clonotype variable

Tobs$Expansion <- "unique"
Tobs$Expansion[Tobs@meta.data$ClonoFreq>5] <- "expanded"

table(Tobs@meta.data$Clonetype)



######TISSUE SHARED CLONOTYPES with tumors########
# clonotypes in Liver
df1<-data.frame(obs1@meta.data$cdr3_aa2,obs1@meta.data$hash.ID,obs1@meta.data$RNA_Group, obs1@meta.data$Tissue)
df2<-na.omit(df1)
df2$cdr3_mouse<-paste0(df2$obs1.meta.data.cdr3_aa2,".",df2$obs1.meta.data.hash.ID,".",df2$obs1.meta.data.RNA_Group)

# clonotypes in Tumor
df3<-data.frame(Tobs@meta.data$cdr3_aa2,Tobs@meta.data$hash.ID,Tobs@meta.data$RNA_Group, Tobs@meta.data$Tissue)
df4<-na.omit(df3)
df4$cdr3_mouse<-paste0(df4$Tobs.meta.data.cdr3_aa2,".",df4$Tobs.meta.data.hash.ID,".",df4$Tobs.meta.data.RNA_Group)

# List of Shared clonotypes
shared_cdr3aa <- intersect(df2$cdr3_mouse, df4$cdr3_mouse)

# Add variable shared in liver
obs1$cdr3_mouse<-paste0(obs1$cdr3_aa2,".",obs1$hash.ID,".",obs1$RNA_Group)
obs1$Shared <- ifelse(obs1$cdr3_mouse %in% shared_cdr3aa, "Shared", "Unique")

# Add variable shared in tumor
Tobs$cdr3_mouse<-paste0(Tobs$cdr3_aa2,".",Tobs$hash.ID,".",Tobs$RNA_Group)
Tobs$Shared <- ifelse(Tobs$cdr3_mouse %in% shared_cdr3aa, "Shared", "Unique")


# Add new variables in liver 
obs1@meta.data$New_division <- obs1@meta.data$NewlablesM_clu
ind<- grepl("CD4",obs1@meta.data$NewlablesM_clu) & grepl("Shared",obs1@meta.data$Shared)
obs1@meta.data$New_division[ind] <- "Shared CD4"
ind<- grepl("CD8",obs1@meta.data$NewlablesM_clu) & grepl("Shared",obs1@meta.data$Shared)
obs1@meta.data$New_division[ind] <- "Shared CD8"
obs1@meta.data$New_division2 <- "Other"
ind<- grepl("CD4",obs1@meta.data$NewlablesM_clu) & grepl("Shared",obs1@meta.data$Shared)
obs1@meta.data$New_division2[ind] <- "Shared CD4"

unique(obs1@meta.data$New_division)

obs1@meta.data$New_division_CD8 <- "Other"
obs1@meta.data$New_division_CD8[obs1@meta.data$New_division=="Shared CD8"] <- "Shared CD8"

# Variable with shared CD4 only in selected groups
obs1$New_division3<-paste0(obs1$New_division2,".",obs1$RNA_Group)
obs1@meta.data$New_division3[obs1@meta.data$New_division3 %in% "Shared CD4.OVA.IFNa"] <- "Other.OVA.IFNa"
obs1@meta.data$New_division3[obs1@meta.data$New_division3 %in% "Shared CD4.OVA"] <- "Other.OVA"

# Add new variables in tumor
Tobs@meta.data$New_division <- Tobs@meta.data$NewlablesM_clu
ind<- grepl("CD4",Tobs@meta.data$NewlablesM_clu) & grepl("Shared",Tobs@meta.data$Shared)
Tobs@meta.data$New_division[ind] <- "Shared CD4"
ind<- grepl("CD8",Tobs@meta.data$NewlablesM_clu) & grepl("Shared",Tobs@meta.data$Shared)
Tobs@meta.data$New_division[ind] <- "Shared CD8"
Tobs@meta.data$New_division2 <- "Other"
ind<- grepl("CD4",Tobs@meta.data$NewlablesM_clu) & grepl("Shared",Tobs@meta.data$Shared)
Tobs@meta.data$New_division2[ind] <- "Shared CD4"


#######################Plot UMAP################

# Save plot
pdf(paste0(pdir, "/AKTPF_TandNK_umap_Fig6B-1.pdf"), width=7, height=6) # Figure 6B-1
DimPlot(obs1, reduction = "umap.harmony.orig.ident", group.by = "NewlablesM_clu", pt.size = 0.5, label = T, label.size = 3, repel = T)
dev.off()

# DimPlot(obs1, reduction = "umap.harmony.orig.ident",group.by = "Expansion", pt.size = 0.5, label = F, order = c("expanded", "unique"), cols = c("lightgrey", "blue"))

pdf(paste0(pdir, "/AKTPF_tetramer_umap_Fig6B-2.pdf"), width=7, height=6) # Figure 6B-2
DimPlot(obs1, reduction = "umap.harmony.orig.ident", group.by = "New_division_CD8", pt.size = 0.5, label = F, order = c("Shared CD8", "Other"), cols = c("lightgrey", "blue"))
dev.off()

pdf(paste0(pdir, "/AKTPF_shared_cd4_umap_FigS12F.pdf"), width=7, height=6) # Figure S12F
DimPlot(obs1, reduction = "umap.harmony.orig.ident", group.by = "New_division2", pt.size = 0.5, label = F, order = c("Shared CD4", "Other"), cols = c("lightgrey", "blue"))
dev.off()

# Export DimPlot as excel
var_list <- c("UMAPh_1", "UMAPh_2","Tissue","orig.ident","Sample","RNA_Group","Tube","RNA_snn_h.orig.ident_res.1.2","NewlablesM_clu","Clonetype","Expansion","New_division","New_division_CD8","New_division2")
TandNK_df <- FetchData(obs1, vars = var_list)
TandNK_df$cell_ID <- rownames(TandNK_df)
write.xlsx(TandNK_df, file=paste0(pdir, "/AKTPF_TandNK_umap_Fig6B_FigS12F.xlsx"), overwrite=T)

#######################Plot UMAP################

# NO_TCR

# # Save plot
# pdf(paste0(pdir, "/AKTPF_TandNK_umap_NO_TCR_FigS13D-1.pdf"), width=7, height=6) # Figure S13D-1
# DimPlot(obs1, reduction = "umap.harmony.orig.ident", group.by = "NewlablesM_clu", pt.size = 0.5, label = T, label.size = 3, repel = T)
# dev.off()
# 
# pdf(paste0(pdir, "/AKTPF_tetramer_umap_NO_TCR_FigS13D-2.pdf"), width=7, height=6) # Figure FigS13D-2
# DimPlot(obs1, reduction = "umap.harmony.orig.ident", group.by = "New_division_CD8", pt.size = 0.5, label = F, order = c("Shared CD8", "Other"), cols = c("lightgrey", "blue"))
# dev.off()
# 
# # Export DimPlot as excel
# var_list <- c("UMAPh_1", "UMAPh_2","Tissue","orig.ident","Sample","RNA_Group","Tube","RNA_snn_h.orig.ident_res.1.2","NewlablesM_clu","Clonetype","Expansion","New_division","New_division_CD8","New_division2")
# TandNK_df <- FetchData(obs1, vars = var_list)
# TandNK_df$cell_ID <- rownames(TandNK_df)
# write.xlsx(TandNK_df, file=paste0(pdir, "/AKTPF_TandNK_umap_NO_TCR_Fig6B_FigS13D.xlsx"), overwrite=T)


#######################DOTPLOT stellare ########################################################
featuresCD4<-c("Ly6c2","Oas1a","Irf8","Ifitm3", # Ifna
               "Il18rap","Il18r1","Il12rb1","Gbp2", # Ifng genes 
               "Ifng","Tbx21","Gzmk","Stat1","Tnf","Cxcr3","Il2", # Th1/effector genes
               "Gata3","Il4","Il5","Lima1","Auh","Runx2","Smad2", # Th2               
               "Rorc","Il17a","Palld","Basp1","Tnfrsf8","Apoe") # Th17 

featuresCD8<-c("Stat1","Ly6c2","Oas1a","Irf8","Ifitm3",   # IFNa genes
               "Il18rap","Il18r1","Ifngas1","Il12rb1","Gbp2", # IL12/IFNg genes
               "Tbx21","Klrg1","Cx3cr1","Tcf7","Klf2","Sell", # pex
                "Gzmb","Gzma","Prf1","Ifng","Tnf","Fasl", # effector
               "Pdcd1","Lag3","Ctla4","Havcr2","Tigit","Il10ra","Tox") # exhaustion

# PlotStellari
pdf(file=paste0(pdir, "/CD8T_expanded_shared_FigS6D.pdf"), width=14, height=3) # Figure S6D
plot.b<-PlotStellare(obs1=obs1, features1=featuresCD8, valueG="Expansion.NewlablesG", i="Shared.CD8 T cells.expanded", group.by0="Orig", controGroup="Liver.TA33", noCtrlvsAll=T)
print(plot.b)
dev.off()

pdf(file=paste0(pdir, "/CD4_T_cells_Liver_FigS6F.pdf"), width=14, height=3) # Figure S6F
plot.c<-PlotStellare(obs1=obs1, features1=featuresCD4, valueG="NewlablesG", i="CD4 T cells", group.by0="Orig", controGroup="Liver.TA33", noCtrlvsAll=T)
print(plot.c)
dev.off()


# Subset putative tumor specific CD8 in liver and tumor
T.liver<-subset(obs1, subset = Tissue=="Liver")
T.liverCD8.shexp<-subset(T.liver, subset = Expansion.NewlablesG =="Shared.CD8 T cells.expanded")

T.tumor<-subset(obs1, subset = Tissue=="Tumor")
T.tumorCD8.shexp<-subset(T.tumor, subset = Expansion.NewlablesG =="Shared.CD8 T cells.expanded")


#########################################################
###############FIGURE 6D############################
##########################################################

##############RADAR PLOT LIVER EXPANDED########
Effector <- list(c("Gzmb","Gzma","Prf1","Ifng","Tnf","Fasl"))
Pex <- list(c("Tbx21","Klrg1","Cx3cr1","Tcf7","Klf2","Sell"))
Exhaustion <- list(c("Pdcd1","Lag3","Ctla4","Havcr2","Tigit","Il10ra","Tox"))
IFNa_genes <- list(c("Stat1","Ly6c2","Oas1a","Irf8","Ifitm3") )
IFNg_genes <-list(c("Il18rap","Il18r1","Ifngas1","Il12rb1","Gbp2"))

T.liverCD8.shexp<- AddModuleScore(T.liverCD8.shexp, features = Effector, name = "Effector")
T.liverCD8.shexp<- AddModuleScore(T.liverCD8.shexp, features = Pex, name = "Pex")
T.liverCD8.shexp<- AddModuleScore(T.liverCD8.shexp, features = Exhaustion, name = "Exhaustion")
T.liverCD8.shexp<- AddModuleScore(T.liverCD8.shexp, features = IFNa_genes, name = "IFNa_genes")
T.liverCD8.shexp<- AddModuleScore(T.liverCD8.shexp, features = IFNg_genes, name = "IFNg_genes")

#Normalize data
#we needed to find a way to normalize the module scores. 
#we used a quantile normalization method.

#M1<-T.liverCD8.shexp@meta.data[,c(16, 37, 54:58)]
M1<-T.liverCD8.shexp@meta.data[,c("NewlablesM", "RNA_Group", "Effector1","Pex1","Exhaustion1","IFNa_genes1","IFNg_genes1")]
M1[,3:7]<-limma::normalizeQuantiles(M1[,3:7])
Groups <-unique(T.liverCD8.shexp$RNA_Group)
CellType<-unique(M1$NewlablesM)
M2<-M1
M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
M4<-rbind(M3a,M3b)
rownames(M4)<-Groups
colnames(M4)<-colnames(M2[,3:7])
x<-as.data.frame(M4)
max1 <- rep(do.call("max", x),5)
min1 <- rep(do.call("min", x),5)
ploTData<- rbind(max1,min1,x)
colors_border=c(rgb(0.1,0.6,0.2,0.9),rgb(0,0,0,0.9))
colors_in=c(rgb(0.1,0.6,0.2,0.1),rgb(0,0,0,0.1))

#radarchart(ploTData, title = "Tumor reactive liver T cells", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty=1, cglwd=0.8)

pdf(paste0(pdir, "/radar.plot.liver.expandedandshared.CD8_newgenes_Fig6D-1.pdf"), width=5, height=5)
radarchart(ploTData, title = "CD8 T cells", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty=1, cglwd=0.8)
dev.off()


##############RADAR PLOT tumor EXPANDED########
Effector <- list(c("Gzmb","Gzma","Prf1","Ifng","Tnf","Fasl"))
Pex <- list(c("Tbx21","Klrg1","Cx3cr1","Tcf7","Klf2","Sell"))
Exhaustion <- list(c("Pdcd1","Lag3","Ctla4","Havcr2","Tigit","Il10ra","Tox"))
IFNa_genes <- list(c("Stat1","Ly6c2","Oas1a","Irf8","Ifitm3") )
IFNg_genes <-list(c("Il18rap","Il18r1","Ifngas1","Il12rb1","Gbp2"))

T.tumorCD8.shexp<- AddModuleScore(T.tumorCD8.shexp, features = Effector, name = "Effector")
T.tumorCD8.shexp<- AddModuleScore(T.tumorCD8.shexp, features = Pex, name = "Pex")
T.tumorCD8.shexp<- AddModuleScore(T.tumorCD8.shexp, features = Exhaustion, name = "Exhaustion")
T.tumorCD8.shexp<- AddModuleScore(T.tumorCD8.shexp, features = IFNa_genes, name = "IFNa_genes")
T.tumorCD8.shexp<- AddModuleScore(T.tumorCD8.shexp, features = IFNg_genes, name = "IFNg_genes")

#Normalize data
#we needed to find a way to normalize the module scores. 
#we used a quantile normalization method.

#M1<-T.tumorCD8.shexp@meta.data[,c(16, 37, 54:58)]
M1<-T.tumorCD8.shexp@meta.data[,c("NewlablesM", "RNA_Group", "Effector1","Pex1","Exhaustion1","IFNa_genes1","IFNg_genes1")]
M1[,3:7]<-limma::normalizeQuantiles(M1[,3:7])
Groups <-unique(T.tumorCD8.shexp$RNA_Group)
CellType<-unique(M1$NewlablesM)
M2<-M1
M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
M4<-rbind(M3a,M3b)
rownames(M4)<-Groups
colnames(M4)<-colnames(M2[,3:7])
x<-as.data.frame(M4)
max1 <- rep(do.call("max", x),5)
min1 <- rep(do.call("min", x),5)
ploTData<- rbind(max1,min1,x)
colors_border=c(rgb(0.1,0.6,0.2,0.9),rgb(0,0,0,0.9))
colors_in=c(rgb(0.1,0.6,0.2,0.1),rgb(0,0,0,0.1))


#radarchart(ploTData, title = "Tumor reactive tumor CD8 T cells", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty=1, cglwd=0.8)

pdf(paste0(pdir, "/radar.plot.tumor.expandedandshared.CD8_Fig6D-2.pdf"), width=5, height=5)
radarchart(ploTData, title = "CD8 T cells", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty=1, cglwd=0.8)
dev.off()



#################FIGURE S6G###########################
############Spider plot liver CD4######################
CD4<-subset(obs1, subset = Tissue =="Liver" &
              NewlablesG =="CD4 T cells")

IFNa_genes<- list(c("Ly6c2","Oas1a","Irf8","Ifitm3") )
IFNg_genes<-list(c("Il18rap","Il18r1","Il12rb1","Gbp2"))
Th1<-list(c("Ifng","Tbx21","Gzmk","Stat1","Tnf","Cxcr3","Il2"))
Th2<-list(c("Gata3","Il4","Il5","Lima1","Auh","Runx2","Smad2"))
Th17<-list(c("Rorc","Il17a","Palld","Basp1","Tnfrsf8","Apoe"))

CD4 <- AddModuleScore(CD4, features = Th1, name = "Th1")
CD4 <- AddModuleScore(CD4, features = Th2, name = "Th2")
CD4 <- AddModuleScore(CD4, features = Th17, name = "Th17")
CD4 <- AddModuleScore(CD4, features = IFNa_genes, name = "IFNa_genes")
CD4 <- AddModuleScore(CD4, features = IFNg_genes, name = "IFNg_genes")

#Normalize data
#we needed to find a way to normalize the module scores. 
#we used a quantile normalization method.

#M1<-CD4@meta.data[,c(16, 37, 54:58)]
M1<-CD4@meta.data[,c("NewlablesM", "RNA_Group", "Th11","Th21","Th171","IFNa_genes1","IFNg_genes1")]
M1[,3:7]<-limma::normalizeQuantiles(M1[,3:7])
Groups <-unique(CD4$RNA_Group)
CellType<-unique(M1$NewlablesM)
M2<-M1
M3a<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[1],3:7]),na.rm=T)
M3b<-colMedians(as.matrix(M2[M2$RNA_Group%in%Groups[2],3:7]),na.rm=T)
M4<-rbind(M3a,M3b)
rownames(M4)<-Groups
colnames(M4)<-colnames(M2[,3:7])
x<-as.data.frame(M4)
max1 <- rep(do.call("max", x),5)
min1 <- rep(do.call("min", x),5)
ploTData<- rbind(max1,min1,x)
colors_border=c(rgb(0.1,0.6,0.2,0.9),rgb(0,0,0,0.9))
colors_in=c(rgb(0.1,0.6,0.2,0.1),rgb(0,0,0,0.1))

pdf(paste0(pdir, "/radar.plot.liver_CD4_FigS6G.pdf"), width=5, height=5)
radarchart(ploTData, title = "CD4 T liver cells", pcol=colors_border, pfcol=colors_in, plwd=2, plty=1, cglcol="grey", cglty = 1,cglwd=0.8)
dev.off()



#########################################
###################GSEA##################
#########################################
# Load signature
#miDB_sig4.MLS <- readRDS(paste0(username, "Dropbox (HSR Global)/90-756142658_RNAseq_Squadrito/05-DGE-Tumor/miDB_sig4.MLS.rds"))
miDB_sig4.MLS <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_NotaroTAD_A018/reference/miDB_sig4.MLS.rds")


###############Figure S6#############################
CD8<-subset(obs1, subset= NewlablesG =="CD8 T cells")

CD8<-SetIdent(CD8, value = "Orig")
obs1.dge<-FindMarkers(CD8, ident.1 = "Liver.TA33.Combo", ident.2 = "Liver.TA33",min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1)
obs1.ranks <- obs1.dge$avg_log2FC
names(obs1.ranks) <- rownames(obs1.dge)

obs2.dge<-FindMarkers(CD8, ident.1 = "Tumor.TA33.Combo", ident.2 = "Tumor.TA33",min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1)
obs2.ranks <- obs2.dge$avg_log2FC
names(obs2.ranks) <- rownames(obs2.dge)


res.combo.liver<-fgsea(pathways = miDB_sig4.MLS, stats = obs1.ranks,minSize=0,maxSize=500, eps=0)
res.combo.tumor<-fgsea(pathways = miDB_sig4.MLS, stats = obs2.ranks,minSize=0,maxSize=500, eps=0)

# Select the columns of interest
res.combo.liver2<-res.combo.liver[,c(1,3,6,8)]
res.combo.tumor2<-res.combo.tumor[,c(1,3,6,8)]

# Merge in a single file and save xlsx
merge1 <- merge(x=res.combo.liver2,y=res.combo.tumor2,by="pathway", suffixes = c("-liver", "-tumor"))
write.xlsx(merge1, paste0(pdir, "GSEA_CD8.xlsx"), asTable=T, rowNames=T)


##################GSEA Heatmap of selected#######################################

# Rearrange data to plot heatmap
res.combo.liver$origin <- "Liver"
res.combo.tumor$origin <- "Tumor"


# res.OVA$origin<- factor(x = res.OVA$origin, levels = c("OVA.Combo","OVA.Il12","OVA.IFNa","liOVA")) 
merge.results<-rbind(res.combo.liver,res.combo.tumor)

#create variable with star for statistic
merge.results$stats <- ifelse(merge.results$padj < 0.0005, "***",
                              ifelse(merge.results$padj < 0.005, "**",
                                     ifelse(merge.results$padj < 0.05, "*", "")))


#pathway manually selected
a1<-c("HALLMARK_G2M_CHECKPOINT",
      "GOBP_DNA_REPLICATION",
      "HALLMARK_ANGIOGENESIS",
      "GOBP_RESPONSE_TO_TYPE_I_INTERFERON",
      "GOBP_RESPONSE_TO_INTERFERON_GAMMA",
      "GOBP_LEUKOCYTE_MEDIATED_CYTOTOXICITY",
      "GOBP_CELL_CYCLE_G2_M_PHASE_TRANSITION",
      "GOBP_DEFENSE_RESPONSE_TO_VIRUS",
      "GOBP_CELL_KILLING",
      "GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_MIGRATION",
      "GOBP_VASCULAR_ENDOTHELIAL_GROWTH_FACTOR_SIGNALING_PATHWAY",
      "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION")

#subset selcting pathways
merge.results1 <-merge.results[merge.results$pathway %in% a1,]

# Rename pathways
merge.results1[merge.results1 == "GOBP_RESPONSE_TO_TYPE_I_INTERFERON"] <- "Response to IFNa"
merge.results1[merge.results1 == "GOBP_RESPONSE_TO_INTERFERON_GAMMA"] <- "Response to IFNg"
merge.results1[merge.results1 == "GOBP_DEFENSE_RESPONSE_TO_VIRUS"] <- "Response to virus"
merge.results1[merge.results1 == "GOBP_LEUKOCYTE_MEDIATED_CYTOTOXICITY"] <- "Leukocyte cytotoxicity"
merge.results1[merge.results1 == "GOBP_CELL_CYCLE_G2_M_PHASE_TRANSITION"] <- "G2-M phase transition"
merge.results1[merge.results1 == "GOBP_CELL_KILLING"] <- "Cell Killing"
merge.results1[merge.results1 == "GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_MIGRATION"] <- "Positive regulation of endothelial cell migration"
merge.results1[merge.results1 == "GOBP_VASCULAR_ENDOTHELIAL_GROWTH_FACTOR_SIGNALING_PATHWAY"] <- "VEGF signaling"
merge.results1[merge.results1 == "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"] <- "Hallmark EMT"
merge.results1[merge.results1 == "HALLMARK_G2M_CHECKPOINT"] <- "Hallmark G2M checkpoint"
merge.results1[merge.results1 == "GOBP_DNA_REPLICATION"] <- "DNA replication"
merge.results1[merge.results1 == "HALLMARK_ANGIOGENESIS"] <- "Hallmark angiogenesis"

# Order in the way we want to see
merge.results1$pathway<- factor(merge.results1$pathway, levels = rev(c("Response to virus",
                                                                   "Response to IFNa",
                                                                   "Response to IFNg",
                                                                   "Cell Killing",
                                                                   "Leukocyte cytotoxicity",
                                                                   "G2-M phase transition",
                                                                   "Hallmark G2M checkpoint",
                                                                   "DNA replication",
                                                                   "Hallmark EMT",
                                                                   "Positive regulation of endothelial cell migration",
                                                                   "VEGF signaling",
                                                                   "PGE2 signature Cilenti et al",
                                                                   "Hallmark angiogenesis")))


# Heatmap
y<-ggplot(merge.results1, aes (x=origin, y=pathway, fill=NES)) + 
  geom_tile() +
  scale_fill_distiller(palette = "RdBu") +
  geom_text(aes(label = stats), vjust = 0.7,angle = 0) +
  labs(y='',x='')+
  theme_classic()+
  theme(axis.text.y = element_text(color = "Black"),
        axis.text.x = element_text(angle = 45, hjust=1, color = "Black"),
        axis.line.x=element_line(color="Black"),
        axis.line.y=element_line(color="Black")) +
  ggtitle("GSEA_CD8")
ggsave(filename = paste0(pdir, "/GSEA_CD8_FigS6C.pdf"), plot=y, width=5, height=2.5)


# Unbias heatmap
merge.results <- merge.results[order(merge.results$padj, decreasing = F)]
merge.results_temp <-merge.results[merge.results$padj < 0.0001,]
filter_in <- head(unique(merge.results_temp$pathway), 100)
merge.results2 <-merge.results[merge.results$pathway %in% filter_in,]
merge.results2$origin<- factor(x = merge.results2$origin, levels = c("Liver","Tumor"))

y<-ggplot(merge.results2, aes (x=origin, y=pathway, fill=NES)) + 
  geom_tile() +
  scale_fill_distiller(palette = "RdBu") +
  geom_text(aes(label = stats), vjust = 0.7,angle = 0) +
  labs(y='',x='')+
  theme_classic()+
  theme(axis.text.y = element_text(color = "Black"), text = element_text(size = 15), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, color = "Black"), 
        axis.line.x=element_line(color="Black"),
        axis.line.y=element_line(color="Black"))
ggtitle("GSEA_CD8")
ggsave(filename = paste0(pdir, "/GSEA_CD8_unbias.pdf"), plot=y, width=15, height=30)



###################Figure S6E################
CD4<-subset(obs1, subset= NewlablesG =="CD4 T cells")

CD4<-SetIdent(CD4, value = "Orig")
obs1.dge<-FindMarkers(CD4, ident.1="Liver.TA33.Combo", ident.2="Liver.TA33", min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1)
obs2.dge<-FindMarkers(CD4, ident.1="Tumor.TA33.Combo", ident.2="Tumor.TA33", min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1)

obs1.ranks <- obs1.dge$avg_log2FC
obs2.ranks <- obs2.dge$avg_log2FC

names(obs1.ranks) <- rownames(obs1.dge)
names(obs2.ranks) <- rownames(obs2.dge)


res.combo.liver <- fgsea(pathways = miDB_sig4.MLS, stats=obs1.ranks, minSize=0, maxSize=500, eps=0)
res.combo.tumor <- fgsea(pathways = miDB_sig4.MLS, stats=obs2.ranks, minSize=0, maxSize=500, eps=0)


# Select columns of interest
res.combo.liver2<-res.combo.liver[,c(1,3,6,8)]
res.combo.tumor2<-res.combo.tumor[,c(1,3,6,8)]

#merge in a single file and save xlsx
merge1 <- merge(x=res.combo.liver2,y=res.combo.tumor2,by="pathway", suffixes = c("-liver", "-tumor"))
write.xlsx(merge1, file = paste0(pdir, "/GSEA_CD4.xlsx"), asTable=T, rowNames=T)


##################GSEA Heatmap of selected#######################################

#rearrange data to plot heatmap
res.combo.liver$origin <- "Liver"
res.combo.tumor$origin <- "Tumor"


# res.OVA$origin<- factor(x = res.OVA$origin, levels = c("OVA.Combo","OVA.Il12","OVA.IFNa","liOVA")) 
merge.results<-rbind(res.combo.liver,res.combo.tumor)

#create variable with star for statistic
merge.results$stats <- ifelse(merge.results$padj < 0.0005, "***",
                              ifelse(merge.results$padj < 0.005, "**",
                                     ifelse(merge.results$padj < 0.05, "*", "")))


#pathway manually selected
a1<-c("HALLMARK_INTERFERON_ALPHA_RESPONSE",
      "HALLMARK_INTERFERON_GAMMA_RESPONSE",
      "LPS_RO",
      "GOBP_INTERLEUKIN_18_PRODUCTION",
      "GOBP_INTERFERON_GAMMA_PRODUCTION",
      "GOBP_T_CELL_DIFFERENTIATION_INVOLVED_IN_IMMUNE_RESPONSE",
      "HALLMARK_E2F_TARGETS",
      "HALLMARK_G2M_CHECKPOINT",
      "GOBP_NEGATIVE_REGULATION_OF_MACROPHAGE_ACTIVATION",
      "IL4_RO",
      "GOBP_RESPONSE_TO_INTERLEUKIN_17")

# #subset selcting pathways
merge.results1 <-merge.results[merge.results$pathway %in% a1,]

# ###rename pathways
merge.results1[merge.results1 == "HALLMARK_INTERFERON_ALPHA_RESPONSE"] <- "Response to IFNa"
merge.results1[merge.results1 == "HALLMARK_INTERFERON_GAMMA_RESPONSE"] <- "Response to IFNg"
merge.results1[merge.results1 == "LPS_RO"] <- "LPS signature (Cilenti et al)"
merge.results1[merge.results1 == "GOBP_INTERLEUKIN_18_PRODUCTION"] <- "IL-18 production"
merge.results1[merge.results1 == "GOBP_INTERFERON_GAMMA_PRODUCTION"] <- "Ifng production"
merge.results1[merge.results1 == "GOBP_T_CELL_DIFFERENTIATION_INVOLVED_IN_IMMUNE_RESPONSE"] <- "T cell differentiation"
merge.results1[merge.results1 == "HALLMARK_E2F_TARGETS"] <- "Hallmark E2F targets"
merge.results1[merge.results1 == "HALLMARK_G2M_CHECKPOINT"] <- "Hallmark G2M checkpoint"
merge.results1[merge.results1 == "GOBP_NEGATIVE_REGULATION_OF_MACROPHAGE_ACTIVATION"] <- "Negative regulation of Macrophages activation"
merge.results1[merge.results1 == "IL4_RO"] <- "IL-4 Signature (Cilenti et al)"
merge.results1[merge.results1 == "GOBP_RESPONSE_TO_INTERLEUKIN_17"] <- "Response to IL-17"



# Order in the way we want to see
merge.results1$pathway<- factor(merge.results1$pathway, levels = rev(c("Response to IFNa",
                                                                   "Response to IFNg",
                                                                   "LPS signature (Cilenti et al)",
                                                                   "IL-18 production",
                                                                   "Ifng production",
                                                                   "T cell differentiation",
                                                                   "Hallmark E2F targets",
                                                                   "Hallmark G2M checkpoint",
                                                                   "Negative regulation of Macrophages activation",
                                                                   "IL-4 Signature (Cilenti et al)",
                                                                   "Response to IL-17")))


# Heatmap
y<-ggplot(merge.results1, aes(x=origin, y=pathway, fill=NES)) + 
  geom_tile() +
  scale_fill_distiller(palette = "RdBu") +
  geom_text(aes(label = stats), vjust=0.7, angle=0) +
  labs(y='',x='')+
  theme_classic()+
  theme(axis.text.y = element_text(color = "Black"),
        axis.text.x = element_text(angle=45, hjust=1, color="Black"),
        axis.line.x=element_line(color="Black"),
        axis.line.y=element_line(color="Black")) +
  ggtitle("GSEA_CD4")
ggsave(filename = paste0(pdir, "/GSEA_CD4_T_cells_FigS6E.pdf"), plot=y, width=5, height=3)

# Unbias heatmap
merge.results <- merge.results[order(merge.results$padj, decreasing = F)]
merge.results_temp <-merge.results[merge.results$padj < 0.0001,]
filter_in <- head(unique(merge.results_temp$pathway), 100)
merge.results2 <-merge.results[merge.results$pathway %in% filter_in,]
merge.results2$origin<- factor(x = merge.results2$origin, levels = c("Liver","Tumor"))

y<-ggplot(merge.results2, aes (x=origin, y=pathway, fill=NES)) + 
  geom_tile() +
  scale_fill_distiller(palette = "RdBu") +
  geom_text(aes(label=stats), vjust=0.7, angle=0) +
  labs(y='',x='')+
  theme_classic()+
  theme(axis.text.y = element_text(color = "Black"), text = element_text(size = 15), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, color = "Black"), 
        axis.line.x=element_line(color="Black"),
        axis.line.y=element_line(color="Black")) +
ggtitle("GSEA_CD4")
ggsave(filename = paste0(pdir, "/GSEA_CD4_T_cells_unbias.pdf"), plot=y, width=15, height=30)



####################Figure 6E#####################################
sub.obj<-subset(obs1, subset= Expansion.NewlablesG =="Shared.CD8 T cells.expanded")

sub.obj<-SetIdent(sub.obj, value = "Orig")
obs1.dge<-FindMarkers(sub.obj, ident.1="Liver.TA33.Combo", ident.2="Liver.TA33", min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1)
obs2.dge<-FindMarkers(sub.obj, ident.1="Tumor.TA33.Combo", ident.2="Tumor.TA33", min.pct=0.01, logfc.threshold=0, min.cells.feature=1, min.cells.group=1)

obs1.ranks <- obs1.dge$avg_log2FC
obs2.ranks <- obs2.dge$avg_log2FC

names(obs1.ranks) <- rownames(obs1.dge)
names(obs2.ranks) <- rownames(obs2.dge)


res.combo.liver <- fgsea(pathways = miDB_sig4.MLS, stats = obs1.ranks, minSize=0, maxSize=500, eps=0)
res.combo.tumor <- fgsea(pathways = miDB_sig4.MLS, stats = obs2.ranks, minSize=0, maxSize=500, eps=0)

#select columns of interest
res.combo.liver2<-res.combo.liver[,c(1,3,6,8)]
res.combo.tumor2<-res.combo.tumor[,c(1,3,6,8)]

#merge in a single file and save xlsx
merge1 <- merge(x=res.combo.liver2,y=res.combo.tumor2,by="pathway", suffixes = c("-liver", "-tumor"))

setwd(pdir)
write.xlsx(merge1,"GSEA_LandT_SharedexpandedCD8.xlsx",asTable = TRUE,rowNames=T) ####save file

##################GSEA Heatmap of selected#######################################

#rearrange data to plot heatmap
res.combo.liver$origin<-"Liver"
res.combo.tumor$origin<-"Tumor"


# res.OVA$origin<- factor(x = res.OVA$origin, levels = c("OVA.Combo","OVA.Il12","OVA.IFNa","liOVA")) 
merge.results<-rbind(res.combo.liver,res.combo.tumor)

#create variable with star for statistic
merge.results$stats <- ifelse(merge.results$padj < 0.0005, "***",
                              ifelse(merge.results$padj < 0.005, "**",
                                     ifelse(merge.results$padj < 0.05, "*", "")))

#merge.results2 <- merge.results[order(merge.results$padj, decreasing = F)]
#filter_in <- head(unique(merge.results2$pathway), 50)
#merge.results1 <-merge.results[merge.results$pathway %in% filter_in,]

# Manually selected pathways
a1_OLD<-c("Exhaustion_Wherry",
      "GOBP_RESPONSE_TO_TYPE_I_INTERFERON",
      "GOBP_RESPONSE_TO_INTERFERON_GAMMA",
      "GOBP_LEUKOCYTE_MEDIATED_CYTOTOXICITY",
      "GOBP_CELL_CYCLE_G2_M_PHASE_TRANSITION",
      "GOBP_DEFENSE_RESPONSE_TO_VIRUS",
      "GOBP_CELL_KILLING",
      "GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_MIGRATION",
      "GOBP_VASCULAR_ENDOTHELIAL_GROWTH_FACTOR_SIGNALING_PATHWAY",
      "PGE2_RO",
      "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION")

a1<-c("GOBP_RESPONSE_TO_VIRUS",
      "HALLMARK_INTERFERON_GAMMA_RESPONSE",
      "HALLMARK_INTERFERON_ALPHA_RESPONSE",
      "HALLMARK_E2F_TARGETS",
      "GOMF_ANTIGEN_BINDING",
      "GOBP_VASCULATURE_DEVELOPMENT",
      "LPS_RO",
      "GOBP_EPITHELIAL_TO_MESENCHYMAL_TRANSITION",
      "GOBP_MESENCHYMAL_CELL_DIFFERENTIATION",
      "Exhaustion_Wherry")

# Subset selected pathways
merge.results1 <-merge.results[merge.results$pathway %in% a1,]

# Rename pathways
merge.results1[merge.results1 == "GOBP_RESPONSE_TO_VIRUS"] <- "Response to virus"
merge.results1[merge.results1 == "HALLMARK_INTERFERON_ALPHA_RESPONSE"] <- "Response to IFNa"
merge.results1[merge.results1 == "HALLMARK_INTERFERON_GAMMA_RESPONSE"] <- "Response to IFNg"
merge.results1[merge.results1 == "LPS_RO"] <- "LPS signature from Cilenti et al"
merge.results1[merge.results1 == "HALLMARK_E2F_TARGETS"] <- "E2F targets"
merge.results1[merge.results1 == "GOMF_ANTIGEN_BINDING"] <- "Antigen binding"
merge.results1[merge.results1 == "GOBP_EPITHELIAL_TO_MESENCHYMAL_TRANSITION"] <- "Epithelial to mesenchymal transition"
merge.results1[merge.results1 == "GOBP_MESENCHYMAL_CELL_DIFFERENTIATION"] <- "Mesenchymal cell differentiation"
merge.results1[merge.results1 == "GOBP_VASCULATURE_DEVELOPMENT"] <- "Vasculature development"
merge.results1[merge.results1 == "Exhaustion_Wherry"] <- "Exhaustion signature (Wherry)"



# Order in the way we want to see
merge.results1$pathway <- factor(merge.results1$pathway, levels = rev(c("Response to virus", # rev()?
                                                                   "Response to IFNa",
                                                                   "Response to IFNg",
                                                                   "LPS signature from Cilenti et al",
                                                                   "E2F targets",
                                                                   "Antigen binding",
                                                                   "Epithelial to mesenchymal transition",
                                                                   "Mesenchymal cell differentiation",
                                                                   "Vasculature development",
                                                                   "Exhaustion signature (Wherry)")))

#merge.results1$origin <- factor(x = merge.results1$origin, levels = c("Liver","Tumor"))

# Heatmap
y<-ggplot(merge.results1, aes (x=origin, y=pathway, fill=NES)) + 
  geom_tile() +
  scale_fill_distiller(palette="RdBu") +
  geom_text(aes(label = stats), vjust=0.7, angle=0) +
  labs(y='',x='')+
  theme_classic()+
  theme(axis.text.y = element_text(color="Black"),
        axis.text.x = element_text(angle=45, hjust=1, color="Black"),
        axis.line.x=element_line(color="Black"),
        axis.line.y=element_line(color="Black")) +
  ggtitle("GSEA_LandT_putative_Tumor_specific")
ggsave(filename = paste0(pdir, "/GSEA_LandT_putative_Tumorspecific_Fig6E.pdf"), plot=y, width=5, height=2.5) #width=5, height=2.5)

# Unbias heatmap
merge.results <- merge.results[order(merge.results$padj, decreasing = F)]
merge.results_temp <-merge.results[merge.results$padj < 0.001,]
#filter_in <- head(unique(merge.results_temp$pathway), 100)
filter_in <- head(unique(merge.results$pathway), 100)
merge.results2 <-merge.results[merge.results$pathway %in% filter_in,]
merge.results2$origin<- factor(x = merge.results2$origin, levels = c("Liver","Tumor"))

y<-ggplot(merge.results2, aes (x=origin, y=pathway, fill=NES)) + 
  geom_tile() +
  scale_fill_distiller(palette = "RdBu") +
  geom_text(aes(label = stats), vjust = 0.7,angle = 0) +
  labs(y='',x='')+
  theme_classic()+
  theme(axis.text.y = element_text(color = "Black"), text = element_text(size = 15), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, color = "Black"), 
        axis.line.x=element_line(color="Black"),
        axis.line.y=element_line(color="Black"))
ggtitle("GSEA_LandT_putative_Tumorspecific")
ggsave(filename = paste0(pdir, "/GSEA_LandT_putative_Tumorspecific_unbias.pdf"), plot=y, width=15, height=30)


##################################################################################
#############################ALLUVIAL PLOT########################################
##################################################################################

#####Figure 6C#########
df1 <- obs1@meta.data
df1<-df1[df1$NewlablesG == "CD8 T cells",]
#df2<- na.omit(df1[,c(12,37,38)])
df2<- na.omit(df1[,c("cdr3_aa2","RNA_Group","Tissue")])
for (i in unique(df1$RNA_Group)){
  df3 <- df2[df2$RNA_Group==i&df2$Tissue=="Liver",]
  df4 <- df2[df2$RNA_Group==i&df2$Tissue=="Tumor",]
  df3 <- df3 %>%  group_by(cdr3_aa2) %>% mutate(ClonoFreq = n()) %>% unique() %>% ungroup() %>% as.data.frame() 
  df4 <- df4 %>%  group_by(cdr3_aa2) %>% mutate(ClonoFreq = n()) %>% unique() %>% ungroup() %>% as.data.frame() 
  #df5 <- merge(df3,df4,by=1,suffixes = c(".Liver",".Tumor"))
  df6 <- rbind(df3,df4)
  df6 <- df6[df6$ClonoFreq>2,]
  df6 <- df6[order(df6$cdr3_aa2),]
  df6$col1[!duplicated(df6$cdr3_aa2)] <- "#CCCCCC"
  Dtests <- duplicated(df6$cdr3_aa2)|duplicated(df6$cdr3_aa2,fromLast = T)
  df6$col1 [Dtests] <- "#CC6600" #"#3399FF"
  plotA <- ggplot(data = df6,
                  aes(x = Tissue, y=ClonoFreq,alluvium = cdr3_aa2,stratum= cdr3_aa2)) + 
    geom_flow(aes(fill = col1),decreasing = F)+
    geom_stratum(aes(fill = col1,size=0.2),decreasing = F) +
    scale_fill_manual(values = setNames(unique(df6$col), unique(df6$col)))+
    theme_bw() + theme(legend.position = "none")+
    scale_size_identity()+
    ggtitle(i)
  nameR <- paste0(i,".plot") 
  assign(nameR,plotA)}

pdf(paste0(pdir, "/Fluvial.Plot.19.CD8mag2clones_Fig6C.pdf"), width=4, height=2)
grid.arrange(TA33.plot,TA33.Combo.plot, ncol = 2, nrow = 1)
#grid.arrange(sapply(plots.all,get), ncol = 2, nrow = 2)
dev.off()

#####Figure 6G#########
df1 <- obs1@meta.data
df1<-df1[df1$NewlablesG == "CD4 T cells",]
#df2<- na.omit(df1[,c(12,37,38)])
df2<- na.omit(df1[,c("cdr3_aa2","RNA_Group","Tissue")])
for (i in unique(df1$RNA_Group)){
  df3 <- df2[df2$RNA_Group==i&df2$Tissue=="Liver",]
  df4 <- df2[df2$RNA_Group==i&df2$Tissue=="Tumor",]
  df3 <- df3 %>%  group_by(cdr3_aa2) %>% mutate(ClonoFreq = n()) %>% unique() %>% ungroup() %>% as.data.frame() 
  df4 <- df4 %>%  group_by(cdr3_aa2) %>% mutate(ClonoFreq = n()) %>% unique() %>% ungroup() %>% as.data.frame() 
  #df5 <- merge(df3,df4,by=1,suffixes = c(".Liver",".Tumor"))
  df6 <- rbind(df3,df4)
  df6 <- df6[df6$ClonoFreq>1,]
  df6 <- df6[order(df6$cdr3_aa2),]
  df6$col1[!duplicated(df6$cdr3_aa2)] <- "#CCCCCC"
  Dtests <- duplicated(df6$cdr3_aa2)|duplicated(df6$cdr3_aa2,fromLast = T)
  df6$col1 [Dtests] <- "#CC6600" #"#3399FF"
  plotA <- ggplot(data = df6,
                  aes(x = Tissue, y=ClonoFreq,alluvium = cdr3_aa2,stratum= cdr3_aa2)) + 
    geom_flow(aes(fill = col1),decreasing = F)+
    geom_stratum(aes(fill = col1,size=0.2),decreasing = F) +
    scale_fill_manual(values = setNames(unique(df6$col), unique(df6$col)))+
    theme_bw() + theme(legend.position = "none")+
    scale_size_identity()+
    ggtitle(i)
  nameR <- paste0(i,".plot") 
  assign(nameR,plotA)}

pdf(paste0(pdir, "/Fluvial.Plot.19.CD4mag1clones_Fig6G.pdf"), width=4, height=2)
grid.arrange(TA33.plot,TA33.Combo.plot, ncol = 2, nrow = 1)
#grid.arrange(sapply(plots.all,get), ncol = 2, nrow = 2)
dev.off()