---
title: "IMC Analysis"
author: "Lena Cords"
output: html_notebook
---
```{r Load libraries}
library(plyr)
library(dplyr)
library(SingleCellExperiment)
library(scater)
library(scran)
library(ggplot2)
library(Rtsne)
library(RColorBrewer)
library(viridis)
library(tidyverse)
library(tables)
library(devtools)
library(data.table)
library(palr)
library(Seurat)
library(ggplot2)
library(gridExtra)
library(RColorBrewer)
library(ggridges)
library(ComplexHeatmap)
library(ggthemes)
```

```{r Set seed}
set.seed(101100)
```


```{r Folder Setup}
wd <- dirname(getwd())
data_folder <- file.path(wd, 'cpout')
results_folder <-(file.path(wd,'results'))
fig_folder <-(file.path(wd,'final-plots'))
```


```{r Read Data}
#RAWDATA
sce <-readRDS(file=file.path(results_folder, paste("sce_CurioniProject_LC.rds")))

#sce <-readRDS(file=file.path(data_folder, paste("sce_CurioniProject_LC.RDS")))
sce <-readRDS(file=file.path(results_folder, paste("sce-clustering-allsamples-original.RDS")))

sce <-readRDS(file=file.path(results_folder, paste("sce-clustering-allsamples-original_final.RDS")))
unique(sce$PatientID)

#
saveRDS(sce, file=file.path(results_folder, paste("sce-clustering-allsamples-original_final.RDS")))
#saveRDS(sce, file=file.path(results_folder, paste("sce-clustering-sub-samples.RDS")))
```

```{r counts}
cofactor=1
#assay(sce, "exprs") <- asinh(counts(sce)/cofactor)
assay(sce, "exprs") <- asinh((assay(sce,"counts"))/cofactor)

exprs(sce) <-asinh((assay(sce,"counts"))/cofactor)


#zscale data
assay(sce, "zscaled") <- t(scale(t(assay(sce, "exprs"))))

#Therapy Status
sce$TherapyStatus <-recode(sce$TherapyStatus, pre="Response", post="Resitance")
sce$TherapyStatus <-factor(sce$TherapyStatus, levels=c("Response", "Resitance"))
```

```{r, Image quality check, fig.width=7, fig.height=7}
x <- data.frame(t(exprs(sce)),Image = colData(sce)$ImageNumber, Batch = colData(sce)$BatchID, Sample=colData(sce)$SampleID, ROI=colData(sce)$ROI, ImmuneStatus=colData(sce)$ImmuneStatus)
x$cell <-rownames(x)
x_long <-pivot_longer(x, cols = c("Histone.H3":"HLA.ABC"), names_to = "target", values_to = "counts")

p <- ggplot(data=x_long, aes(x=counts, y=Sample))+
  geom_density_ridges(aes(color=Sample))+
  facet_wrap(~target, scales="fixed", ncol=11)+
  theme(strip.background = element_blank(),
        panel.background=element_rect(fill='white', colour = 'black'),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        plot.background=element_blank(),
        legend.key = element_blank(),
        axis.text.x = element_text(angle = 90, hjust = 1))
suppressMessages(plot(p))
ggsave(plot=p,filename = file.path(results_folder, paste("marker density over samples.pdf")), width=12, height=10, dpi=300)
```


```{r run PCA}
#run PCA
sce<-runPCA(sce, exprs_values="exprs")
plotPCA(sce)
dim(reducedDim(sce, "PCA"))

#define number of PC used for PCA based on graph
perc.var <- attr(reducedDim(sce, "PCA"), "percentVar")
plot(perc.var, xlab="PC", ylab="Variance explained (%)")

#define PC here and save into reduced Dimensions under new PCA
reducedDim(sce, "PCA_20") <- reducedDim(sce, "PCA")[,1:20]
plotReducedDim(sce, dimred="PCA", ncomponents=4)

#run tSNE
sce <-runTSNE(sce, dimred="PCA_25", perplexity=50)
p <- plotTSNE(sce,colour_by="SampleID",run_args = list(perplexity=100))
p <-plotTSNE(sce,colour_by="BatchID",run_args = list(perplexity=100))

#run UMAP
sce <-runUMAP(sce, exprs_values="exprs")
p<- plotUMAP(sce, colour_by="SampleID")
p <- plotUMAP(sce, colour_by="BatchID")
plot(p)
```

```{r Specify channels used for clustering}
all.marker <- rownames(sce)
bad.marker <- c("Histone H3", "Iridium", "Iridium2", "CD134")
good.marker <- all.marker[!all.marker %in% bad.marker]
print(good.marker)

#Heatmap colour scheme
hmcol<-rev(brewer.pal(11,"RdBu"))
cluscol <-viridis_pal()
```


```{r clustering all cells k15}
g.jaccard_15.sub <- buildSNNGraph(sce,type="jaccard", k=15, subset.row=good.marker)
clust.louvain_15.sub <- igraph::cluster_louvain(g.jaccard_15.sub)$membership
length(unique(clust.louvain_15.sub))

#add clusters to colData
sce$Cluster_15.sub <-factor(clust.louvain_15.sub)
saveRDS(sce, file=file.path(results_folder, paste("sce-clustering-allsamples-original_final.RDS")))
```

```{r, plot heatmap k15 all cells}
agg_sce_15.sub <-aggregateAcrossCells(sce,ids=sce$Cluster_15.sub, average=TRUE, use_exprs_values="zscaled" )

plotHeatmap(agg_sce_15.sub,
            features = good.marker,
            exprs_values = "zscaled",
            symmetric = FALSE,
            zlim=c(-3,3),
            color=hmcol,
            sortrowss=FALSE,
            show_colnames=TRUE,
            colour_columns_by = "Cluster_15.sub",
            main="Cluster k=15.sub",
            width=10, height=10,
            filename=file.path(results_folder, paste("heatmap_mean_cluster_15-subset.pdf")))
```

```{r all cells run tsne and umap}
perpl <- c(seq(from=10, to=100, by=10))
perpl <-c(90)

#tSNE
for (i in perpl){
  sce <-runTSNE(sce,
                exprs_value="exprs",
                PCA=T,
                perplexity=i,
                name=paste0("TSNE_k",i)
                )
  
}
saveRDS(sce, file=file.path(results_folder, paste("sce-clustering-allsamples-original_final.RDS")))

#UMAP
for(i in perpl){
sce <- runUMAP(sce,
                 exprs_values = "exprs",
                 name = paste0("UMAP_p", i),
                 n_neighbors = i)
}
saveRDS(sce, file=file.path(results_folder, paste("sce-clustering-allsamples-original_final.RDS")))
```

```{r cluster assignment all cells}
is.myeloid <- data.frame("Cluster"=c(8,9,10,11,12,15,16,21,23),"Phenotype"="Myeloid")
is.tumour <- data.frame("Cluster"=c(6,7,13,18,19,20,24),"Phenotype"="Tumour")
is.tcells <- data.frame("Cluster"=c(1,2,14,22,25),"Phenotype"="T cell")
is.undefined <- data.frame("Cluster"=c(3),"Phenotype"="Undefined")
is.stroma <- data.frame("Cluster"=c(5,26),"Phenotype"="Stroma")
is.vessel <- data.frame("Cluster"=c(4,17),"Phenotype"="Vessel")

cluster.phenotypes <-rbind(is.tumour,is.myeloid, is.tcells, is.undefined, is.stroma,is.vessel)
cluster.phenotypes <-cluster.phenotypes[order(cluster.phenotypes$Cluster),]
rownames(cluster.phenotypes) <-cluster.phenotypes$Cluster

#merge into sce
levels.cluster15.sub <- sce$Cluster_15.sub
levels(levels.cluster15.sub) <-factor(cluster.phenotypes$Phenotype)
sce$Cluster_levels <- levels.cluster15.sub

sce$Cluster_levels <-factor(sce$Cluster_levels, levels=c("Tumour", "Stroma","Vessel","T cell","Myeloid","Undefined"))
unique(sce$Cluster_levels)
saveRDS(sce, file=file.path(results_folder, paste("sce-clustering-allsamples-original_final.RDS")))
```

```{r, QC marker expression on tsne all cells, fig.width=12, fig.height=12}
dat <-as.data.frame(reducedDims(sce)$`TSNE_k90`)
dat$cell <- rownames(dat)
dat.counts <-as.data.frame(t((assay(sce,"exprs"))))
dat.counts$cell <- rownames(dat.counts)
dat.all <-merge(dat, dat.counts, by.x="cell", by.y="cell")
dat.all.long <-pivot_longer(dat.all, cols = good.marker, names_to = "target", values_to = "counts")


p <-dat.all.long %>%
  ggplot(aes(x=V1, y=V2, color=counts))+
  facet_wrap(~target, scales = "free", ncol = 8)+
  geom_point(alpha=0.5, size=0.2)+
  scale_color_gradientn(colours=rev(brewer.pal(11, 'Spectral')), name='Counts')+
  ggtitle('')+
  theme(strip.background = element_blank(),
        axis.line=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background=element_blank(),
        panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        plot.background=element_blank(),
         strip.text = element_text(size=10)) 
plot(p)
ggsave(filename=file.path(results_folder, paste("all-Marker_tSNE_k90.png")), plot=p, width=12, height=10, dpi=300)
```

```{r, plot clustering outcome all cells and clinical data on tsne}
tsne<- "TSNE_k90"
vp <-c("BatchID","SampleID","Cluster_15.sub","TherapyStatus","PatientID","Cluster_levels")
vp <-c("TherapyStatus")
vp <-c("Cluster_levels")

for( i in vp){
  #No label
  p <- plotReducedDim(sce,
                      paste0(tsne),
                      colour_by=paste0(i))+
      labs(x="tSNE_1", y="tSNE_2")+
      theme(legend.title= element_blank())+
      #scale_colour_manual(values=brewer.pal(8,"Dark2")[3:8])
      scale_colour_manual(values=c("grey","black"))
  plot(p)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_NOlabel.png")), width=7, height=4, dpi=300)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_NOlabel.svg")), width=7, height=4, dpi=300)

 # 
  #With label
  p <- plotReducedDim(sce,
                      paste0(tsne),
                      colour_by=paste0(i),
                      text_by=paste0(i))+
      labs(x="tSNE_1", y="tSNE_2")+
      theme(legend.title= element_blank())+
      #scale_colour_manual(values=brewer.pal(8,"Dark2")[3:8])
      scale_colour_manual(values=c("grey","black"))
  plot(p)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_label.png")), width=7, height=4, dpi=300)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_label.svg")), width=7, height=4, dpi=300)
}
```

```{r, plotting marker expression on TSNE Seurat package, fig.width=18, fig.height=12, dpi=300}
sce.seurat <- as.Seurat(sce, counts = "counts", data = "exprs")
#good.marker 
p<-FeaturePlot(object = sce.seurat, features = good.marker,cols = c("white", "black"), reduction="TSNE_k90", combine=T, label.size = 5, ncol=8) &NoAxes()
plot(p)
ggsave(filename=file.path(fig_folder, paste("AAA_Final_AllMarker_AllCells_Seurat.png")), plot=p, width=24, height=16, dpi = 300)
ggsave(filename=file.path(fig_folder, paste("AAA_Final_AllMarker_AllCells_Seurat.svg")), plot=p, width=24, height=16, dpi = 300)

```

```{r, Plot data per sample config}
#plot_folder=file.path(results_folder, 'plots')
dat.all <-as_tibble(colData(sce))

saveRDS(dat.all, file=file.path(results_folder, paste("data-allsamples-analysis.RDS")))
dat.all<-readRDS(file=file.path(results_folder, paste("data-allsamples-analysis.RDS")))
#Make factor to change order for plots

dat.all$TherapyStatus <-as.factor(dat.all$TherapyStatus)
#dat.all$TherapyStatus <-revalue(dat.all$TherapyStatus, c("pre"="Response", "post"="Resistance"))
#dat.all$TherapyStatus <-factor(dat.all$TherapyStatus, levels=c("Pre", "Post"))

dat.all$ImmuneStatus <-as.factor(dat.all$ImmuneStatus)
dat.all$ImmuneStatus <-factor(dat.all$ImmuneStatus, levels=c("T-cell low", "T-cell high"))

dat.all$PatientID <-as.factor(dat.all$PatientID)
dat.all$PatientID <-factor(dat.all$PatientID, levels=c("2","4","5","9","11","12","14"))


#dat.all$TherapyStatus <-recode(dat.all$TherapyStatus, pre="Pre", post="Post")
dat.all$ImmuneStatus <-recode(dat.all$ImmuneStatus, Tcell_low="T-cell low", Tcell_high="T-cell high")
dat.all$Cluster_levels <-recode(dat.all$Cluster_levels, Tcells ="T cell", Vessels="Vessel")
unique(dat.all$Cluster_levels)
```

```{r, Descriptive bar plots for all cells}
colour.cells <-c("red","grey","blue","orange","darkgreen","yellow")

dat.all$TherapyStatus <-recode(dat.all$TherapyStatus, pre="Response", post="Resistance")
dat.all$TherapyStatus <- factor(dat.all$TherapyStatus, levels = c("Response", "Resistance"))
dat.all$Cluster_levels <-recode(dat.all$Cluster_levels, Tcells="T cell", Vessels="Vessel")

dat.all$Cluster_levels <-factor(dat.all$Cluster_levels, levels=c("Tumour", "Stroma","Vessel","T cell","Myeloid","Undefined"))

#Plot comparison of cell type fraction per Therapy state
library(wesanderson)
colour.cells <- brewer.pal(n=8,name="Dark2")
colour.cells <- brewer.pal(8,"Dark2")[3:8]
p <-
  ggplot(dat.all) +
  aes(x = PatientID, fill = factor(Cluster_levels)) +
  geom_bar(position = "fill")+ coord_flip() +facet_wrap("TherapyStatus") +
  ggtitle("Celltype distribution over therapy status")+
  labs(y="Celltype (%)",x="Patient", fill="Celltype")+
  scale_fill_manual(values=colour.cells)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
p

ggsave(filename=file.path(fig_folder, paste("Final_Celltypes_overPatients_therapyStatus2.pdf")), plot=p)

#Cell Number over patients and therapy
p1 <- ggplot(dat.all)+ aes(x=PatientID, fill=factor(TherapyStatus))+
  geom_bar(position = position_dodge())+ facet_wrap("Cluster_levels")+
  ggtitle("Celltypes over Therapy Status")+
  labs(x="Patient", y="Cell number", fill="Therapy Status")+
  scale_fill_manual(values=c("grey","black"))+
    theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
p1
ggsave(filename=file.path(results_folder, paste("Final_Celltypes over Patients and therapy status_barplot.pdf")), plot=p1)
```

```{r, Subset T cells}
tcell.marker <- c("CD45RO", "FOXP3","PD-1" ,"CD45RA","CD4" ,"HLA-ABC" , "VISTA" , "CD8a","Granzyme B","ICOS", "IDO", "PD-L1" ,"GITR" ,  "CD3", "CD27",  "TCF1 / TCF7" , "LAG-3" , "CD45", "CCL21")
t <- sce[, sce$Cluster_levels == "Tcells"]

t.cells <- sce[, sce$Cluster_levels == "T cell"]
#save T cell sce
saveRDS(t.cells,file=file.path(results_folder, paste("sce-Tcells.RDS")))
t.cells <-readRDS(file=file.path(results_folder, paste("sce-Tcells.RDS")))
```

```{r, Cluster T cells}
g.jaccard_50.tcells <- buildSNNGraph(t.cells,type="jaccard", k=50, subset.row=tcell.marker)
clust.louvain_50.tcells <- igraph::cluster_louvain(g.jaccard_50.tcells)$membership
length(unique(clust.louvain_50.tcells))

#add clusters to colData
t.cells$Cluster_50.tcells <-factor(clust.louvain_50.tcells)

table(clust.louvain_50.tcells)
```

```{r, plot HM T cells}
agg_sce_50.tcells <-aggregateAcrossCells(t.cells,ids=t.cells$Cluster_50.tcells, average=TRUE, use_exprs_values="zscaled" )

plotHeatmap(agg_sce_50.tcells,
            features = tcell.marker,
            exprs_values = "zscaled",
            symmetric = FALSE,
            zlim=c(-3,3),
            color=hmcol,
            sortrowss=FALSE,
            show_colnames=TRUE,
            annotation_legend=T,
            colour_columns_by = "Cluster_50.tcells",
            main="T-cell cluster",
            width=10, height=10)
```

```{r, T cell Type and Category}
types_Tcells.k50 <- data.frame("Type"=c("naive T cell", "naive T cell", "CD8+ T cell", "regulatory T cell", "exhausted CD4+ T cell","naive T cell", "activated CD4+ T cell", "naive T cell", "exhausted CD8+ T cell", "exhausted CD8+ T cell","exhausted CD8+ T cell"), "Cluster"=c(1:11))
types_Tcells.k50$Type <-factor(types_Tcells.k50$Type, levels=c("naive T cell", "CD8+ T cell","activated CD4+ T cell", "regulatory T cell", "exhausted CD4+ T cell","exhausted CD8+ T cell"))

#Tcells.k50_category <- data.frame("Category"=c("CD4+ T cell","CD8+ T cell", "CD8+ T cell","CD4+", "CD4+","T cell", "CD4+", "T cell","CD8+ T cell","CD8+ T cell","CD8+ T cell"), "Cluster"=c(1:11))

Tcells.k50_category <- data.frame("Category"=c("CD4+ T cell","CD8+ T cell", "CD8+ T cell","CD4+ T cell", "CD4+ T cell","T cell", "CD4+", "T cell","CD8+ T cell","CD8+ T cell","CD8+ T cell"), "Cluster"=c(1:11))

Tcells.k50_category$Category <-factor(Tcells.k50_category$Category)
```

```{r, HM t cells}
agg_sce <-aggregateAcrossCells(t.cells,ids=t.cells$Cluster_50.tcells, average=TRUE, use_exprs_values="zscaled" )
agg_sce <- agg_sce[rownames(agg_sce) %in% tcell.marker,]

hm.genes <- rownames(agg_sce)
hm.val <- assay(agg_sce,"zscaled")[hm.genes,]
censor = T
#option 1
if (censor==T){
# Z-score data
# Censor z-score at 2
hm.val[hm.val > 2] =2
hm.val[hm.val < -2] =-2
}

ncells.stroma <-as.data.frame(table(t.cells$Cluster_50.tcells))
colnames(ncells.stroma) <-c("ids","Frequency")

col.viridis = list(Cluster = viridis(length(unique(ncells.stroma$ids))), Type= c("#D6DAE7","#A202E6","#0F43FC","#FF9201","#1BFF00","#FF180D"))

#col.viridis = list(Cluster = viridis(length(unique(ncells.stroma$ids))), Type= tcell.colour)
tcell.colour <-brewer.pal(n = 6, name = "Dark2")
tcell.colour <-c("#D6DAE7","#A202E6","#0F43FC","#FF9201","#1BFF00","#FF180D")
tcell.colour <- palette("Okabe-Ito")[2:7]
col.viridis = list(Cluster = viridis(length(unique(ncells.stroma$ids))), Type= tcell.colour)
names(col.viridis$Cluster) <- unique(ncells.stroma$ids)
names(col.viridis$Type) <- unique(types_Tcells.k50$Type)


col.viridis
#with cluster line
ht1<-Heatmap(hm.val, name=" ", col=hmcol,clustering_method_rows = "ward.D2",clustering_method_columns = "ward.D2",top_annotation = HeatmapAnnotation(Frequency=anno_barplot(ncells.stroma$Frequency, gp=gpar(fill = "grey")),Cluster=ncells.stroma$ids,col=col.viridis, Cluster=types_Tcells.k50$Cluster,Type=types_Tcells.k50$Type))
draw(ht1)

#without cluster line
ht1<-Heatmap(hm.val, name=" ", col=hmcol,clustering_method_rows = "ward.D2",clustering_method_columns = "ward.D2",top_annotation = HeatmapAnnotation(Frequency=anno_barplot(ncells.stroma$Frequency, gp=gpar(fill = "grey")),col=col.viridis ,Type=types_Tcells.k50$Type))
draw(ht1)

#Save complex heatmap as pdf
pdf(file=file.path(fig_folder, paste("Final_HM_Tcells_k50_Types_woCluster.pdf")), width=7, height=7)
draw(ht1)
dev.off()
```

```{r, Tcell Types to sce}
#Add T cell types to tcell sce
type.tcell <- t.cells$Cluster_50.tcells
levels(type.tcell) <-factor(types_Tcells.k50$Type)
t.cells$Tcell_Type <-type.tcell

#Add T cell category to tcell sce
type.tcell <- t.cells$Cluster_50.tcells
levels(type.tcell) <-factor(Tcells.k50_category$Category)
t.cells$Tcell_Category <-type.tcell
```

```{r run tsne and umap tcells}
perpl <- c(seq(from=10, to=100, by=10))
perpl <-c(60)

#tSNE
for (i in perpl){
  t.cells <-runTSNE(t.cells,
                exprs_value="exprs",
                PCA=T,
                perplexity=i,
                name=paste0("tcell-TSNE_k",i)
                )
  
}
saveRDS(t.cells,file=file.path(results_folder, paste("sce-Tcells.RDS")))

#UMAP
for(i in perpl){
t.cells <- runUMAP(t.cells,
                 exprs_values = "exprs",
                 name = paste0("cell_UMAP_k", i),
                 n_neighbors = i)
}
saveRDS(t.cells,file=file.path(results_folder, paste("sce-Tcells.RDS")))
```

```{r t cell config}
display.brewer.pal(n = 6, name = 'Dark2')
tcell.colour <-brewer.pal(n = 6, name = "Dark2")
tcell.colour <-c("#D6DAE7","#A202E6","#0F43FC","#FF9201","#1BFF00","#FF180D")
tcell.colour <- palette("R4")[1:6]

t.cells$TherapyStatus <- revalue(t.cells$TherapyStatus, c("pre"="Response", "post"="Resistance"))
t.cells$TherapyStatus <-factor(t.cells$TherapyStatus, levels=c("Response", "Resistance"))
t.cells$Tcell_Type <-factor(t.cells$Tcell_Type, levels=c("naive T cell", "CD8+ T cell","activated CD4+ T cell", "regulatory T cell", "exhausted CD4+ T cell","exhausted CD8+ T cell"))
t.cells$PatientID <-as.factor(t.cells$PatientID)
```

```{r plot tsne and umap for T cells with clustering and clinical data}
tsne<- "tcell-TSNE_k70"
vp <-c("BatchID","SampleID","Cluster_50.tcells","TherapyStatus","PatientID","Tcell_Category", "Tcell_Type")
vp <- "PatientID"
vp <- "Tcell_Type"
vp <- "TherapyStatus"

for( i in vp){
  #No label
  p <- plotReducedDim(t.cells,
                      paste0(tsne),
                      colour_by=paste0(i))+
      labs(x="tSNE_1", y="tSNE_2")+
      theme(legend.title= element_blank())+
    #scale_colour_brewer(palette="Set1")
    #scale_colour_manual(values=tcell.colour)
  scale_colour_manual(values=c("grey","black"))
  plot(p)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_NOlabel.svg")), width=7, height=4, dpi=300)
    ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_NOlabel.png")), width=7, height=4, dpi=300)

  #With label
  p <- plotReducedDim(t.cells,
                      paste0(tsne),
                      colour_by=paste0(i),
                      text_by=paste0(i))+
      labs(x="tSNE_1", y="tSNE_2")+
      theme(legend.title= element_blank())+
    #scale_colour_brewer(palette="Set1")
   #scale_colour_manual(values=tcell.colour)
    scale_colour_manual(values=c("grey","black"))

  plot(p)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_label.svg")), width=7, height=4, dpi=300)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_label.png")), width=7, height=4, dpi=300)
}
```

```{r Plot T cell data}
dat.tcells <-as_tibble(colData(t.cells))

dat.tcells$TherapyStatus <-as.factor(dat.tcells$TherapyStatus)

dat.tcells$ImmuneStatus <-as.factor(dat.tcells$ImmuneStatus)
dat.tcells$ImmuneStatus <-factor(dat.tcells$ImmuneStatus, levels=c("Tcell_low", "Tcell_high"))

dat.tcells$PatientID <-as.factor(dat.tcells$PatientID)
dat.tcells$PatientID <-factor(dat.tcells$PatientID, levels=c("2","4","5","9","11","12","14"))

dat.tcells$TherapyStatus <-recode(dat.tcells$TherapyStatus, pre="Pre", post="Post")
dat.tcells$ImmuneStatus <-recode(dat.tcells$ImmuneStatus, Tcell_low="T-cell low", Tcell_high="T-cell high")
dat.tcells$Cluster_levels <-recode(dat.tcells$Cluster_levels, Tcells ="T cell", Vessels="Vessel")
unique(dat.tcells$Cluster_levels)

#Subsetting T cell Data into high and low to be able to compare pre and Resistance between them
dat.tcells.high <-filter(dat.tcells, ImmuneStatus=="T-cell high")
dat.tcells.low <-filter(dat.tcells, ImmuneStatus=="T-cell low")

dat.tcells$Tcell_Type <-factor(dat.tcells$Tcell_Type, levels=c("naive T cell", "CD8+ T cell","activated CD4+ T cell", "regulatory T cell", "exhausted CD4+ T cell","exhausted CD8+ T cell"))

#Cluster
p <-
  ggplot(dat.tcells) +
  aes(x = PatientID, fill = factor(Cluster_50.tcells)) +
  geom_bar(position = "fill")+ coord_flip()+
  facet_wrap("TherapyStatus")+
  ggtitle("T cell cluster over patients and therapy status")+
  labs(y="Cluster (%)",x="Patient", fill="Cluster")+
  #scale_fill_viridis("discrete"=T)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
p

ggsave(filename = file.path(fig_folder, paste("Final_T cell cluster over Patients and Therapy Status_all.PDF")), plot=p)

#Type
p <-
  ggplot(dat.tcells) +
  aes(x = PatientID, fill = factor(Tcell_Type)) +
  geom_bar(position = "fill")+ coord_flip()+
  facet_wrap("TherapyStatus")+
  ggtitle("T-cell type over patients and therapy status")+
  labs(y="Cell Type (%)",x="Patient", fill="Type")+
  scale_fill_manual(values=tcell.colour)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
p

ggsave(filename = file.path(fig_folder, paste("Final_T cell Type over Patients and Therapy Status_all_2.PDF")), plot=p)

#Category
p <-
  ggplot(dat.tcells) +
  aes(x = PatientID, fill = factor(Tcell_Category)) +
  geom_bar(position = "fill")+ coord_flip()+
  facet_wrap("TherapyStatus")+
  ggtitle("T-cell category over patients and therapy status")+
  labs(y="Cell Category (%)",x="Patient", fill="Category")+
  scale_fill_viridis("discrete"=T)
p

ggsave(filename = file.path(results_folder, paste("T cell Category over Patients and Therapy Status_all.PDF")), plot=p)



#Cell Number over patients and therapy
p1 <- ggplot(dat.tcells)+ aes(x=PatientID, fill=factor(TherapyStatus))+
  geom_bar(position = position_dodge())+ facet_wrap("Tcell_Type", scales="free")+
  ggtitle("Celltypes over Therapy Status")+
  labs(x="Patient", y="Cell number", fill="Therapy Status")+
  scale_fill_manual(values=c("grey","black"))+
    theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
p1
ggsave(filename=file.path(results_folder, paste("Final_TcellTypes_overPatients-TherapyStatus_barplot_freescales.pdf")), plot=p1)
```

```{r Plot t cell data seurat}
tcell.seurat <- as.Seurat(t.cells, counts = "counts", data = "exprs")
#good.marker 
p<-FeaturePlot(object = tcell.seurat, features = c("CD3","CD4","CD8a","GITR","IDO","LAG-3"),cols = c("white", "black"), reduction="tcell.TSNE_k70", combine=T, label.size = 10, ncol=3) &NoAxes()

#p<-FeaturePlot(object = tcell.seurat, features = tcell.marker,cols = c("white", "black"),reduction="tcell.TSNE_k60", combine=T, label.size = 10, ncol=5) &NoLegend() &NoAxes()

plot(p)
ggsave(plot=p, file=file.path(fig_folder, paste("Final_T_cell_AllMarker_TSNEk70.pdf")), width=12, height=7, dpi=300)
ggsave(plot=p, file=file.path(fig_folder, paste("Final_T_cell_AllMarker_TSNEk70.svg")), width=12, height=7, dpi=300)
```

```{r, QC marker expression on tsne tcells, fig.width=8, fig.height=8}
dat <-as.data.frame(reducedDims(t.cells)$`tcell-TSNE_k70`)
dat$cell <- rownames(dat)
dat.counts <-as.data.frame(t((assay(t.cells,"exprs"))))
dat.counts$cell <- rownames(dat.counts)
dat.all <-merge(dat, dat.counts, by.x="cell", by.y="cell")
dat.all.long <-pivot_longer(dat.all, cols = tcell.marker, names_to = "target", values_to = "counts")


p <-dat.all.long %>%
  ggplot(aes(x=V1, y=V2, color=counts))+
  facet_wrap(~target, scales = "free", ncol = 5)+
  geom_point(alpha=0.5, size=0.2)+
  scale_color_gradientn(colours=rev(brewer.pal(11, 'Spectral')), name='Counts')+
  ggtitle('')+
  theme(strip.background = element_blank(),
        axis.line=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background=element_blank(),
        panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        plot.background=element_blank(),
         strip.text = element_text(size=10)) 
plot(p)
ggsave(filename=file.path(results_folder, paste("tcell-Marker_tcells_tSNE_k90.png")), plot=p, width=12, height=10, dpi=300)
```
```{r Subset myeloids}
myeloid.sce <- sce[, sce$Cluster_levels == "Myeloid"]

#save Myeloid sce
saveRDS(myeloid.sce,file=file.path(results_folder, paste("sce-myeloid.RDS")))
myeloid.sce <-readRDS(file=file.path(results_folder, paste("sce-myeloid.RDS")))
```

```{r Myeloid Marker}
myeloid.marker <- c("PD-L1" , "CD45", "CD68",  "CD163","CD15", "CD14", "CD38" ,"Siglec-1 / CD169", "CD303" ,"HLA-DR" ,"CD11b", "PD-1", "CD206 (MMR)","TMEM173","CD204")
```


```{r Subclustering on myeloids k25}
g.jaccard_25.myeloid <- buildSNNGraph(myeloid.sce,type="jaccard", k=25, subset.row=myeloid.marker)
clust.louvain_25.myeloid <- igraph::cluster_louvain(g.jaccard_25.myeloid)$membership
length(unique(clust.louvain_25.myeloid))

#add clusters to colData
myeloid.sce$Cluster_25.myeloid <-factor(clust.louvain_25.myeloid)

#plot heatmap
agg_sce_25.myeloid <-aggregateAcrossCells(myeloid.sce,ids=myeloid.sce$Cluster_25.myeloid, average=TRUE, use_exprs_values="zscaled" )
plotHeatmap(agg_sce_25.myeloid, features = myeloid.marker, exprs_values = "zscaled",symmetric = FALSE,zlim=c(-3,3),color=hmcol,sortrowss=FALSE, show_colnames=TRUE,colour_columns_by = "Cluster_25.myeloid", main="Cluster k=25.myeloid",width=10, height=10,filename=file.path(results_folder, paste("heatmap_mean_cluster_25-myeloid.pdf")))

saveRDS(myeloid.sce,file=file.path(results_folder, paste("sce-myeloid.RDS")))
myeloid.sce <-readRDS(file=file.path(results_folder, paste("sce-myeloid.RDS")))
```


```{r add myeloid types}
is.myeloid <- data.frame("Cluster"=c(1,13), "Phenotype"="Myeloid")
is.monocytes <- data.frame("Cluster"=c(2,8), "Phenotype"="Monocyte")
is.undefined <-data.frame("Cluster"=c(3,10),"Phenotype"="Undefined")
is.macrophage <- data.frame("Cluster"=c(4,5,6,12,14,15),"Phenotype"="Macrophage")
is.DC <- data.frame("Cluster"=c(7,9),"Phenotype"="Dendritic Cell")
is.gran <-data.frame("Cluster"=c(11),"Phenotype"="Granulocyte")
is.pDC <-data.frame("Cluster"=c(16),"Phenotype"="pDC")

myeloid.cluster.k25 <-rbind(is.myeloid, is.monocytes, is.undefined, is.macrophage, is.DC, is.gran, is.pDC)
myeloid.cluster.k25 <-myeloid.cluster.k25[order(myeloid.cluster.k25$Cluster),]
rownames(myeloid.cluster.k25) <-myeloid.cluster.k25$Cluster

levels.myeloid.cluster.25 <- myeloid.sce$Cluster_25.myeloid
levels(levels.myeloid.cluster.25) <-factor(myeloid.cluster.k25$Phenotype)
myeloid.sce$Myeloid_Phenotype_k25 <-levels.myeloid.cluster.25
```

```{r subset macrophages}
macrophage.marker <- c("PD-L1" , "CD45", "CD68",  "CD163", "CD14", "CD38" ,"Siglec-1 / CD169","HLA-DR", "PD-1", "CD206 (MMR)","TMEM173","CD204" )

macrophage.sce <- myeloid.sce[, myeloid.sce$Myeloid_Phenotype_k25 == "Macrophage"]
#save T cell sce
saveRDS(macrophage.sce,file=file.path(results_folder, paste("sce-Macrophages.RDS")))
macrophage.sce <- readRDS(file=file.path(results_folder, paste("sce-Macrophages.RDS")))
```

```{r Subclustering on Macrophages k85}
g.jaccard_85.macrophage <- buildSNNGraph(macrophage.sce,type="jaccard", k=85, subset.row=macrophage.marker)
clust.louvain_85.macrophage <- igraph::cluster_louvain(g.jaccard_85.macrophage)$membership
length(unique(clust.louvain_85.macrophage))

macrophage.sce$Cluster_85.macrophage <-factor(clust.louvain_85.macrophage)
length(unique(macrophage.sce$Cluster_85.macrophage))
```

```{r HM  Macrophages k85}
agg_sce_85.macrophage <-aggregateAcrossCells(macrophage.sce,ids=macrophage.sce$Cluster_85.macrophage, average=TRUE, use_exprs_values="zscaled" )

plotHeatmap(agg_sce_85.macrophage,
            features = macrophage.marker,
            exprs_values = "zscaled",
            symmetric = FALSE,
            zlim=c(-3,3),
            color=hmcol,
            sortrowss=FALSE,
            show_colnames=TRUE,
            colour_columns_by = "Cluster_85.macrophage",
            main="Cluster k=85.macrophage",
            width=10, height=10,
            filename=file.path(results_folder, paste("heatmap_mean_cluster_85-macrophage.pdf")))
```

```{r, Macrophage Types}
is.M1 <-data.frame("Cluster"=c(1,3,6,7,8,10),"Type"="M1")
is.M2 <-data.frame("Cluster"=c(2,4,5,9),"Type"="M2")
macrophage.types <-rbind(is.M1, is.M2)
macrophage.types <-macrophage.types[order(macrophage.types$Cluster),]
rownames(macrophage.types) <-macrophage.types$Cluster

type.macro <- macrophage.sce$Cluster_85.macrophage
levels(type.macro) <-factor(macrophage.types$Type)
macrophage.sce$Macrophage_Type_k85 <-type.macro
```

```{r, HM Macrophage Types}
agg_sce <-aggregateAcrossCells(macrophage.sce,ids=macrophage.sce$Cluster_85.macrophage, average=TRUE, use_exprs_values="zscaled" )
agg_sce <- agg_sce[rownames(agg_sce) %in% macrophage.marker,]

hm.genes <- rownames(agg_sce)
hm.val <- assay(agg_sce,"zscaled")[hm.genes,]
censor = T
#option 1
if (censor==T){
# Z-score data
# Censor z-score at 2
hm.val[hm.val > 2] =2
hm.val[hm.val < -2] =-2
}

ncells.stroma <-as.data.frame(table(macrophage.sce$Cluster_85.macrophage))
colnames(ncells.stroma) <-c("ids","Frequency")

col.viridis = list(Cluster = viridis(length(unique(ncells.stroma$ids))), Type= c("#A9A9A9","#000000"))
col.viridis = list(Cluster = viridis(length(unique(ncells.stroma$ids))), Type= pal_jco()(2))

names(col.viridis$Cluster) <- unique(ncells.stroma$ids)
names(col.viridis$Type) <- unique(macrophage.types$Type)

col.viridis

#with cluster row
ht1<-Heatmap(hm.val, name=" ", col=hmcol,clustering_method_rows = "ward.D2",clustering_method_columns = "ward.D2",top_annotation = HeatmapAnnotation(Frequency=anno_barplot(ncells.stroma$Frequency,gp=gpar(fill = "grey")),Cluster=ncells.stroma$ids,col=col.viridis, Type=macrophage.types$Type)) 

#without cluster row
ht1<-Heatmap(hm.val, name=" ", col=hmcol,clustering_method_rows = "ward.D2",clustering_method_columns = "ward.D2",top_annotation = HeatmapAnnotation(Frequency=anno_barplot(ncells.stroma$Frequency,gp=gpar(fill = "grey")),col=col.viridis, Type=macrophage.types$Type)) 
draw(ht1)

#Save complex heatmap as pdf
pdf(file=file.path(fig_folder, paste("Final_HM_MacrophageTypes_k85_wCluster.pdf")), width=7, height=7)
draw(ht1)
dev.off()
```

```{r, macrophage markers seurat tsne}
macrophage.sce <-runUMAP(macrophage.sce, ncomponents=2, n_neighbors=20, exprs_values="exprs", name="UMAP_k20")
macrophage.seurat <- as.Seurat(macrophage.sce, counts = "counts", data = "exprs")
#good.marker 

p<-FeaturePlot(object = macrophage.seurat, features = macrophage.marker,cols = c("white", "black"),reduction="UMAP_k20", combine=T, label.size = 10, ncol=4) &NoAxes()

plot(p)
ggsave(plot=p, file=file.path(fig_folder, paste("Final_Macrophages_AllMarker_UMAPk20.png")), width=12, height=7, dpi=300)
ggsave(plot=p, file=file.path(fig_folder, paste("Final_Macrophages_AllMarker_UMAPk20.svg")), width=12, height=7, dpi=300)

```


#TSNE for subset of macrophages
```{r Check clusters on tsne and UMAP}
perpl <- c(seq(from=10, to=100, by=10))
perpl <-c(10,20)

#tSNE
for (i in perpl){
  macrophage.sce <-runTSNE(macrophage.sce,
                exprs_value="exprs",
                PCA=T,
                perplexity=i,
                name=paste0("macrophage-TSNE_k",i)
                )
  
}
saveRDS(macrophage.sce,file=file.path(results_folder, paste("sce-Macrophages.RDS")))

#UMAP
for(i in perpl){
macrophage.sce <- runUMAP(macrophage.sce,
                 exprs_values = "exprs",
                 name = paste0("macrophage_UMAP_k", i),
                 n_neighbors = i)
}
saveRDS(macrophage.sce,file=file.path(results_folder, paste("sce-Macrophages.RDS")))
```


```{r compare tsne & UMAP + clustering outcomes for macrophages}
tsne<- "UMAP_k20"
#tsne<- "macrophage-TSNE_10"
macrophage.sce$PatientID <-factor(macrophage.sce$PatientID)
vp <-c("BatchID","SampleID","Cluster_85.macrophage","TherapyStatus","PatientID","Macrophage_Type_k85")
vp <-c("PatientID")
vp <-c("TherapyStatus")
vp <-c("Macrophage_Type_k85")
for( i in vp){
  #No label
  p <- plotReducedDim(macrophage.sce,
                      paste0(tsne),
                      colour_by=paste0(i))+
      labs(x="tSNE_1", y="tSNE_2")+
      theme(legend.title= element_blank())+
    #scale_colour_brewer(palette="Set1")
    #scale_colour_manual(values=c("black","grey"))
    scale_color_jco()
  plot(p)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_NOlabel.png")), width=7, height=4, dpi=300)
    ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_NOlabel.svg")), width=7, height=4, dpi=300)

  #With label
  p <- plotReducedDim(macrophage.sce,
                      paste0(tsne),
                      colour_by=paste0(i),
                      text_by=paste0(i))+
      labs(x="tSNE_1", y="tSNE_2")+
      theme(legend.title= element_blank())+
  #scale_colour_brewer(palette="Set1")
    #scale_colour_manual(values=c("black","grey"))¨
    scale_color_jco()
  plot(p)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_label.png")), width=7, height=4, dpi=300)
  ggsave(plot=p, file=file.path(fig_folder, paste0(tsne,"_",i,"_label.svg")), width=7, height=4, dpi=300)
}
```

```{r, Macropahge plots}
dat.macrophage <-as_tibble(colData(macrophage.sce))

dat.macrophage$ImmuneStatus <-as.factor(dat.macrophage$ImmuneStatus)
dat.macrophage$ImmuneStatus <-factor(dat.macrophage$ImmuneStatus, levels=c("Tcell_low", "Tcell_high"))

dat.macrophage$PatientID <-as.factor(dat.macrophage$PatientID)
dat.macrophage$PatientID <-factor(dat.macrophage$PatientID, levels=c("2","4","5","9","11","12","14"))

dat.macrophage$TherapyStatus <-recode(dat.macrophage$TherapyStatus, pre="Pre", post="Post")
dat.macrophage$ImmuneStatus <-recode(dat.macrophage$ImmuneStatus, Tcell_low="T-cell low", Tcell_high="T-cell high")
dat.macrophage$Cluster_levels <-recode(dat.macrophage$Cluster_levels, Tcells ="T cell", Vessels="Vessel")
unique(dat.macrophage$Cluster_levels)


#Subsetting T cell Data into high and low to be able to compare pre and Resistance between them
dat.macrophage.high <-filter(dat.macrophage, ImmuneStatus=="T-cell high")
dat.macrophage.low <-filter(dat.macrophage, ImmuneStatus=="T-cell low")

#Macrophages Cluster % Distribution over Therapystatus 
p <-
  ggplot(dat.macrophage) +
  aes(x = PatientID, fill = factor(Cluster_85.macrophage)) +
  geom_bar(position = "fill")+ coord_flip()+
  facet_wrap("TherapyStatus")+
  ggtitle("Macrophage cluster over patients and therapy status")+
  labs(y="Cluster Type (%)",x="Patient", fill="Cluster")+
  #scale_fill_viridis("discrete"=T)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
p

ggsave(filename = file.path(fig_folder, paste("Final_Macrophage cluster over Patients and Therapy Status_all.pdf")), plot=p)

#Macrophage Type distibutiuon over patient
p1 <- ggplot(dat.macrophage, aes(x=PatientID, fill=Macrophage_Type_k85))+
  geom_bar(position = "fill")+
  coord_flip()+
  facet_wrap("TherapyStatus")+
  labs(title="Macropahge Type Distribution over Therapy Status", fill="Type", x="Patient", y="Proportion")+
  scale_fill_jco()+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
  #scale_fill_manual(values=c("#A9A9A9","#000000"))
p1
ggsave(filename=file.path(fig_folder, paste("Macrophage_Type-over_TherapyStatus.pdf")), plot=p1)


#Cell Number over patients and therapy
p1 <- ggplot(dat.macrophage)+ aes(x=PatientID, fill=factor(TherapyStatus))+
  geom_bar(position = position_dodge())+ facet_wrap("Macrophage_Type_k85")+
  ggtitle("Celltypes over Therapy Status")+
  labs(x="Patient", y="Cell number", fill="Therapy Status")+
  scale_fill_manual(values=c("grey","black"))+
    theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(strip.background = element_rect(fill="white"),
              strip.text = element_text(colour="black"))
p1
ggsave(filename=file.path(fig_folder, paste("Final_Macrophage1_2_overPatients-TherapyStatus_barplot.pdf")), plot=p1)
```