4.1-Characterization_cell_populations.Rmd

---
title: "Characterization of cell populations in wounding experiments"
author: "Massimo Andreatta & Laura Yerly"
output: html_document
date: "2024-02-07"
---

```{r setup, include=FALSE,fig.width=16,fig.height=12}
renv::restore()
library(ggplot2)
library(Seurat)
```

# Set paths

```{r}
data.path <- "cache/WHexp.rds"

seu <- readRDS(data.path)
```

## Set up some helper functions
```{r}
get_means_signature<- function(object,signature){
  means <- sapply(levels(object$seurat_clusters), function(x) {
    mean(object@meta.data[object$seurat_clusters == x, signature])
  })
  means
}

get_medians_signature<- function(object,signature){
  medians <- sapply(levels(object$seurat_clusters), function(x) {
    median(object@meta.data[object$seurat_clusters == x, signature])
  })
  medians
}
```

# Look at the data

```{r fig.width=12}
a <- DimPlot(seu, reduction = "umap", group.by = "Sample") + theme(aspect.ratio = 1)
b <- DimPlot(seu, group.by = "scGate_multi", label = T, repel = T, label.size = 2) +
    theme(aspect.ratio = 1) + ggtitle("scGate annotation")
a | b
```

Simplify annotations
```{r}
seu$scGate_broad <- seu$scGate_multi

seu$scGate_broad[seu$scGate_multi %in% c("Macrophage","Monocyte")] <- "MoMac"
seu$scGate_broad[seu$scGate_multi %in% c("CD4T","CD8T","NK")] <- "Tcell"
```

## Signatures to identify cell populations 

```{r}
library(UCell)
signatures<-list()
signatures$BCC<-c("PTCH1","GLI1","GLI2","HHIP","MYCN","SPON2")
signatures$TCells<-c("CD2","CD3D","CD3E")
signatures$NormKerat<-c("KRT14","KRT5","GATA3","KRTDAP")
signatures$Melanocytes<-c("PMEL","MLANA")
signatures$Cycling<-c("MKI67","ORC1")
signatures$MyeloidCells<-c("LYZ","SPI1")
signatures$BCells<-c("CD79A","MS4A1")
signatures$Fibroblasts<-c("COL1A1","IL6","PDGFRA","CXCL12","CD82","CCL19")
signatures$Pericytes<-c("RGS5","KCNE4","CPE")
signatures$MastCells<-c("CTSG","TPSB2","TPSAB1")
signatures$EndothelialCells<-c("PECAM1","HEY1","CCL14","ACKR1","SELE")

seu<-AddModuleScore_UCell(seu,features = signatures, name = "", ncores=8)
```

View signatures
```{r fig.height=12, fig.width=20}
FeaturePlot(seu, features=names(signatures)) & theme(aspect.ratio = 1,
  axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())
```

```{r}
VlnPlot(seu,features=names(signatures),
        group.by = "seurat_clusters", stack = T,flip = T,assay = "RNA")+ NoLegend()
```

```{r fig.height=12, fig.width=20}
FeaturePlot(seu, features=c("PTCH1","GLI1","GLI2","HHIP","MYCN","SPON2")) & theme(aspect.ratio = 1,
  axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

FeaturePlot(seu, features=c("RGS5","KCNE4","CPE","PTPRC")) & theme(aspect.ratio = 1,
  axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

FeaturePlot(seu, features=c("LYZ","SPI1","CSF1R","CD68","APOE","FCN1","S100A9")) & theme(aspect.ratio = 1,
  axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

FeaturePlot(seu, features=c("KRT1","KRT10","SBSN","KRTDAP","KRT5","KRT14")) & theme(aspect.ratio = 1,
  axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())
```


```{r}
DimPlot(seu, reduction = "umap", group.by = "seurat_clusters", label=T) + theme(aspect.ratio = 1)
```




Clusters by majority vote
```{r}
tab <- table(seu$seurat_clusters, seu$scGate_broad, useNA = "ifany")
tab_norm <- apply(tab, 1, function(x) {x/sum(x)})
tab_norm

majority <- apply(tab_norm, 2, function(x){
  names(which.max(x))
})

seu$annotation <- unname(majority[seu$seurat_clusters])

DimPlot(seu, group.by = "annotation", label = T, repel = T, label.size = 2) +
    theme(aspect.ratio = 1)
```
Focus on T cells and Keratinocytes
```{r}
which.cells <- colnames(seu)[seu$annotation %in% c("Tcell","Epithelial") | is.na(seu$annotation)]
sub <- subset(seu, cells=which.cells)
```

## Dimensionality reduction and clustering

```{r fig.width=12}
sub <- NormalizeData(sub) |> FindVariableFeatures(nfeatures=2000) |>
  ScaleData()

sub <- sub |> RunPCA(npcs=30) |> RunUMAP(reduction = "pca", dims = 1:30)

a <- DimPlot(sub, reduction = "umap", group.by = "Sample") + theme(aspect.ratio = 1)
b <- DimPlot(sub, reduction = "umap", group.by = "annotation") + theme(aspect.ratio = 1)

a | b
```

View signatures
```{r fig.height=12, fig.width=20}
FeaturePlot(sub, features=feats) & theme(aspect.ratio = 1,
  axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())
```

Do it by patient?
```{r}
seu.list <- SplitObject(sub, split.by = "Sample")

seu.list <- lapply(seu.list, function(x) {
  x <- NormalizeData(x) |> FindVariableFeatures(nfeatures=2000, verbose = F) |> 
    ScaleData() |> RunPCA(npcs=30, verbose = F) |> RunUMAP(reduction = "pca", dims = 1:30, verbose = F)
  
  x <- FindNeighbors(x, dims=1:30, verbose=F) |> FindClusters(resolution = 2, verbose = F)
  x
})
```

```{r fig.height=12, fig.width=20}
pll <- lapply(seu.list, function(x) {
  FeaturePlot(x, features=names(signatures)) & theme(aspect.ratio = 1,
                                                       axis.title = element_blank(),
                                                       axis.text = element_blank(),
                                                       axis.ticks = element_blank())
})

pll
```

## Get mean expression by cluster
```{r}
BCC_all <- c()
NormKerat_all <- c()
T_all <- c()
Mye_all <- c()

#Define thresholds
for (x in names(seu.list)) {
  sub <- seu.list[[x]]
  BCC_all <- c(BCC_all, get_means_signature(sub,"BCC"))
  NormKerat_all <- c(NormKerat_all, get_means_signature(sub,"NormKerat"))
  T_all <- c(T_all, get_means_signature(sub,"TCells"))
  Mye_all <- c(Mye_all, get_means_signature(sub,"MyeloidCells"))
}

plot(x=BCC_all, y=NormKerat_all)
plot(x=BCC_all, y=T_all)
plot(x=NormKerat_all, y=T_all)
plot(x=NormKerat_all, y=Mye_all)
```

```{r}
#Apply thresholds  
names <- names(seu.list)
seu.list <- lapply(names(seu.list), function(x) {
  
  sub <- seu.list[[x]]
  BCC_means<-get_means_signature(sub,"BCC")
  NormKerat_means<-get_means_signature(sub,"NormKerat")
  TCells_means<-get_means_signature(sub,"TCells")
  Myeloid_means <- get_means_signature(sub,"MyeloidCells")
  
  #Any doublets?
  tcell_pass <- names(TCells_means)[TCells_means>0.1]
  bcc_pass <- names(BCC_means)[BCC_means>0.1]
  ker_pass <- names(NormKerat_means)[NormKerat_means>0.15]
  myeloid_pass <- names(Myeloid_means)[Myeloid_means>0.15]
  
  doublets <- unique(c(intersect(tcell_pass, ker_pass),
                       intersect(tcell_pass, bcc_pass),
                       intersect(myeloid_pass, bcc_pass),
                       intersect(myeloid_pass, ker_pass)))
  
  sub$annotation[sub$seurat_clusters %in% ker_pass] <- "Normal_Kerat"
  sub$annotation[sub$seurat_clusters %in% bcc_pass] <- "Cancer_cells"
  sub$annotation[sub$seurat_clusters %in% doublets] <- "Doublets"
  sub
})  
names(seu.list) <- names
```

```{r fig.height=12, fig.width=20}
pll <- lapply(names(seu.list), function(x) {
  DimPlot(seu.list[[x]], group.by = "annotation", label = T, repel = T, label.size = 2) +
    theme(aspect.ratio = 1) + ggtitle(x)
})
wrap_plots(pll)
```

#Re-merge after annotation
```{r}
sub <- Reduce(f=merge, x=seu.list)

table(sub$Sample, sub$annotation, useNA = "ifany")
```

```{r}
labels.sub <- sub$annotation

seu$annotation[names(labels.sub)] <- unname(labels.sub)

table(seu$Sample, seu$annotation, useNA = "ifany")
```



```{r}
seu$patient_bcc <- seu$Sample
Idents(seu) <- "patient_bcc"
```

```{r fig.width=12}
a <- DimPlot(seu, reduction = "umap", group.by = "Sample") + theme(aspect.ratio = 1)
b <- DimPlot(seu, group.by = "annotation", label = T, repel = T, label.size = 2) +
    theme(aspect.ratio = 1) + ggtitle("Cell type")
a | b
```

# Save

```{r} 
#file <- "M:/DER/LABOS/ONCODER/Laura/1-Experiments/4-Bioinformatic/Code_reproducibility/WHexp_annotated.rds"
file <- "cache/WHexp_annotated.rds"
saveRDS(seu, file)
```

Save cancer cells to file
```{r}
tumor <- subset(seu, subset=annotation == "Cancer_cells")
file_cancer <- 'cache/Tumor_WH.rds'

saveRDS(tumor, file_cancer)
```