04a_analyzeOverlap_seurat_transfer_sct_integrated_pilot.R

##################################################
## Project: Analyze Overlap between Seurat Clusters and Transferred Labels, SCT & Integration
## Date: 26.05.2020
## Author: Nathalie
##################################################

setwd("~/Documents/PostmortemBrain/analysis/markerGeneDefinition")
#setwd("/net/PE1/raid1/LAURA/SC_ANALYSIS/markerGeneDefinition")

library(Seurat)
library(dplyr)
library(reshape2)
library(RColorBrewer)
library(ggplot2)

### PILOT DATA ###################################

dataset <- "pilot"
data <- readRDS(paste0(dataset,"/data_object_sct_integrated.rds"))
subclasses <- read.table("allan_human/subclass_cluster.csv", sep=";", header = TRUE)

## ignore Exclude class
x <- as.numeric(data$seurat_clusters)
y <- as.character(data$predictions_subclass)
i  <- which(y != "Exclude")
x <- x[i]
y <- y[i]
overlap <- table(x, y)

# Plot absolute overlap
overlap_long <- melt(overlap)
overlap_long$x <- as.factor(overlap_long$x)

ggplot(overlap_long, aes(x,y))+
  geom_tile(aes(fill=value))+
  geom_text(aes(label=value))+
  scale_fill_gradient(low="white",
                      high="darkred")+
  theme(panel.grid.major.x=element_blank(), #no gridlines
        panel.grid.minor.x=element_blank(), 
        panel.grid.major.y=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.background=element_rect(fill="white"), # background=white
        axis.text.x = element_text(angle=90, hjust = 1,vjust=1,size = 12,face = "bold"),
        plot.title = element_text(size=15,face="bold"),
        axis.text.y = element_text(size = 12,face = "bold")) + 
  ggtitle("Overlap Louvain Cluster and Transferred Subclasses")+
  theme(legend.title=element_text(face="bold", size=14)) + 
  scale_x_discrete(name="Louvain cluster")+
  scale_y_discrete(name="Transferred Subclass")+
  labs(fill="Overlap")
ggsave(paste0(dataset,"/Overlap_louvain_subclass_sct_integrated.png"))

# Plot relative values from each louvain cluster
rel_ol <- function(m) m / rep(colSums(m), each = nrow(m))
overlap_perc <- rel_ol(t(overlap))

overlap_perc_long <- melt(overlap_perc)
overlap_perc_long$x <- as.factor(overlap_perc_long$x)

ggplot(overlap_perc_long, aes(x,y))+
  geom_tile(aes(fill=value))+
  geom_text(aes(label=round(x = value, digits = 2)))+
  scale_fill_gradient(low="white",
                      high="darkred")+
  theme(panel.grid.major.x=element_blank(), #no gridlines
        panel.grid.minor.x=element_blank(), 
        panel.grid.major.y=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.background=element_rect(fill="white"), # background=white
        axis.text.x = element_text(angle=90, hjust = 1,vjust=1,size = 12,face = "bold"),
        plot.title = element_text(size=15,face="bold"),
        axis.text.y = element_text(size = 12,face = "bold")) + 
  ggtitle("Overlap Louvain Cluster and Transferred Subclasses")+
  theme(legend.title=element_text(face="bold", size=14)) + 
  scale_x_discrete(name="Louvain cluster")+
  scale_y_discrete(name="Transferred Subclass")+
  labs(fill="Overlap")
ggsave(paste0(dataset,"/Overlap_louvain_subclass_sct_integrated.png"))

# Assign labels to louvain clusters according to highest overlap
ass_cluster <- apply(overlap_perc, 2, function(t) 
  rownames(overlap_perc)[which.max(t)])
ass_vec <- ass_cluster[as.character(as.numeric(data$seurat_clusters))]

data <- AddMetaData(
  object = data,
  metadata = ass_vec,
  col.name = "seurat_transfer_subclass"
)

png(paste0(dataset,"/UMAP_labelTransfer_seurat_transfer_subclass_sct_integrated.png"), height = 600, width = 600)
DimPlot(data, reduction = "umap", group.by = "seurat_transfer_subclass", 
        cols = colorRampPalette(brewer.pal(9, "Set1"))(nlevels(as.factor(data$seurat_transfer_subclass))))
dev.off()



# Subset cells from neuronal classes and refine assignment
subclusters_IT <- subclasses[subclasses$subclass_label=="IT",]  # refine only in clusters of IT
data_exc <- subset(data, subset = seurat_transfer_subclass == "IT")
# use only IT clusters
x <- as.numeric(data_exc$seurat_clusters)
y <- as.character(data_exc$predictions_cluster)
i  <- which(y %in% subclusters_IT[,1])
x <- x[i]
y <- y[i]
overlap <- table(x, y)

overlap_perc <- rel_ol(t(overlap))
overlap_perc_long <- melt(overlap_perc)
overlap_perc_long$x <- as.factor(overlap_perc_long$x)

ggplot(overlap_perc_long, aes(x,y))+
  geom_tile(aes(fill=value))+
  geom_text(aes(label=round(x = value, digits = 2)))+
  scale_fill_gradient(low="white",
                      high="darkred")+
  theme(panel.grid.major.x=element_blank(), #no gridlines
        panel.grid.minor.x=element_blank(), 
        panel.grid.major.y=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.background=element_rect(fill="white"), # background=white
        axis.text.x = element_text(angle=90, hjust = 1,vjust=1,size = 12,face = "bold"),
        plot.title = element_text(size=15,face="bold"),
        axis.text.y = element_text(size = 12,face = "bold")) + 
  ggtitle("Overlap Louvain Cluster and Transferred Clusters (IT Subclass)")+
  theme(legend.title=element_text(face="bold", size=14)) + 
  scale_x_discrete(name="Louvain cluster")+
  scale_y_discrete(name="Transferred Subclass")+
  labs(fill="Overlap")
ggsave(paste0(dataset,"/Overlap_louvain_cluster_IT_sct_integrated.png"), width = 12, height = 18)

# Assign labels to louvain clusters according to highest overlap
ass_cluster <- apply(overlap_perc, 2, function(t) 
  rownames(overlap_perc)[which.max(t)])
ass_vec <- ass_cluster[as.character(as.numeric(data_exc$seurat_clusters))]

data_exc <- AddMetaData(
  object = data_exc,
  metadata = ass_vec,
  col.name = "seurat_transfer_cluster_IT"
)



# # Subset cells from neuronal classes and refine assignment
# subclusters_VIP <- subclasses[subclasses$subclass_label=="VIP",]  # refine only in clusters of IT
# data_inh <- subset(data, subset = seurat_transfer_subclass == "VIP")
# # use only IT clusters
# x <- as.numeric(data_inh$seurat_clusters)
# y <- as.character(data_inh$predictions_cluster)
# i  <- which(y %in% subclusters_VIP[,1])
# x <- x[i]
# y <- y[i]
# overlap <- table(x, y)
# 
# overlap_perc <- rel_ol(t(overlap))
# overlap_perc_long <- melt(overlap_perc)
# overlap_perc_long$x <- as.factor(overlap_perc_long$x)
# 
# ggplot(overlap_perc_long, aes(x,y))+
#   geom_tile(aes(fill=value))+
#   geom_text(aes(label=round(x = value, digits = 2)))+
#   scale_fill_gradient(low="white",
#                       high="darkred")+
#   theme(panel.grid.major.x=element_blank(), #no gridlines
#         panel.grid.minor.x=element_blank(), 
#         panel.grid.major.y=element_blank(), 
#         panel.grid.minor.y=element_blank(),
#         panel.background=element_rect(fill="white"), # background=white
#         axis.text.x = element_text(angle=90, hjust = 1,vjust=1,size = 12,face = "bold"),
#         plot.title = element_text(size=15,face="bold"),
#         axis.text.y = element_text(size = 12,face = "bold")) + 
#   ggtitle("Overlap Louvain Cluster and Transferred Clusters (VIP Subclass)")+
#   theme(legend.title=element_text(face="bold", size=14)) + 
#   scale_x_discrete(name="Louvain cluster")+
#   scale_y_discrete(name="Transferred Subclass")+
#   labs(fill="Overlap")
# ggsave(paste0(dataset,"/Overlap_louvain_cluster_VIP_sct_integrated.png"), width = 12, height = 18)
# 
# # Assign labels to louvain clusters according to highest overlap
# ass_cluster <- apply(overlap_perc, 2, function(t) 
#   rownames(overlap_perc)[which.max(t)])
# ass_vec <- ass_cluster[as.character(as.numeric(data_inh$seurat_clusters))]
# 
# data_inh <- AddMetaData(
#   object = data_inh,
#   metadata = ass_vec,
#   col.name = "seurat_transfer_cluster_VIP"
# )

# Merge cluster assignments of IT and VIP to data object
data$sub_cluster <- as.character(data$seurat_transfer_subclass)

data$sub_cluster[Cells(data_exc)] <- data_exc$seurat_transfer_cluster_IT
# data$sub_cluster[Cells(data_inh)] <- data_inh$seurat_transfer_cluster_VIP
png(paste0(dataset,"/UMAP_labelTransfer_sub_cluster_sct_integrated.png"), height = 600, width = 600)
DimPlot(data, reduction = "umap", group.by = "sub_cluster", 
        cols = colorRampPalette(brewer.pal(9, "Set1"))(nlevels(as.factor(data$sub_cluster))))
dev.off()

saveRDS(data, paste0(dataset,"/data_object_sct_integrated.rds"))


# Plot UMAP with donor colour
png(paste0(dataset,"/UMAP_donor_sct_integrated.png"), height = 600, width = 600)
DimPlot(data, reduction = "umap", group.by = "orig.ident", 
        cols = colorRampPalette(brewer.pal(9, "Set1"))(nlevels(as.factor(data$orig.ident))))
dev.off()

# Write cell type annotations to file
ct_list <- data.frame("cellID" = dimnames(data)[[2]],
                      "celltype" = data$sub_cluster)
write.table(ct_list, file = paste0(dataset,"/celltypes.csv"), sep = "\t", quote = FALSE,
            row.names = FALSE)