##########################################################################
####################### Load relevant libraries ##########################
##########################################################################

library(dplyr)
library(tidyr)
library(deMULTIplex)
library(ggplot2)
library(stringr)
library(Seurat)
library(SeuratObject)
library(SingleCellExperiment)
library(scater)
library(scran)
library(tibble)
library(sjmisc)

## read in and filter away Undetermined

counts_mat <- read.csv("C:/Users/Cathal/Documents/Nova_30k/Combined_Nova-LMOfilteredOUT-30kset-_RSEC_MolsPerCell.csv", skip = 7)
# 1 - SMK Calls file
smk_calls <- read.csv(file = "C:/Users/Cathal/Documents/Nova_30k/Nova-LMOfilteredOUT-30kset-_Sample_Tag_Calls.csv", skip=7)
head(smk_calls)

# 2 - Immune Classification cell type file
Immune_class <- read.csv(file = "C:/Users/Cathal/Documents/Nova_30k/Combined_Nova-LMOfilteredOUT-30kset-_cell_type_experimental.csv")
head(Immune_class)

#### merge SMK calls file and Immune class file
df_meta_SevenBridges <- full_join(smk_calls, Immune_class, by = "Cell_Index")
head(df_meta_SevenBridges)

## merge meta to count mat
count_mat_full <- full_join(counts_mat, df_meta_SevenBridges, by = "Cell_Index")

## remove rows that contain "Undetermined" or "Multiplet"
count_mat_filtered <- count_mat_full[!(count_mat_full$Sample_Name=="Undetermined" | count_mat_full$Sample_Name=="Multiplet"),]

## reduced count mat 
count_mat_filtered_only_counts <- count_mat_filtered[,1:428]

## reduced meta data
SB_meta_filtered <- select(count_mat_filtered, Cell_Index, Sample_Tag:Cell_Type_Experimental)
# count_mat_filtered[,429:431]

## add LMO bartable to BD metadata
## Read in LMO barTable
barTable_all <- read.csv("C:/Users/Cathal/Downloads/bar_table_full.csv")
head(barTable_all)

# reduce
barTable_all <- select(barTable_all,-contains("nUMI"))
barTable_all <- select(barTable_all,-contains("label"))
head(barTable_all)

## Join sev bridges data with LMO barTable
#dim(df_meta_SevenBridges)
dim(barTable_all)

##
df_meta_all <- left_join(SB_meta_filtered, barTable_all, by="Cell_Index")


##########

## count matrix
## set rownames
rownames(count_mat_filtered_only_counts) <- count_mat_filtered_only_counts$Cell_Index
count_mat_filtered_only_counts$Cell_Index <- NULL
# check size of mat
dim(count_mat_filtered_only_counts)

# View top of mat
count_mat_filtered_only_counts[1:5,1:5]

## tt - transpose
count_mat_filtered_only_counts_tt <- rotate_df(count_mat_filtered_only_counts)

## Ensure that cells and genes are in the correct position for Seurat
head(colnames(count_mat_filtered_only_counts_tt))

dim(count_mat_filtered_only_counts_tt)

## read in as is

# counts_mat <- read.csv("C:/Users/Cathal/Documents/Nova_30k/Combined_Nova-LMOfilteredOUT-30kset-_RSEC_MolsPerCell.csv", skip = 7)
# 
# ## set rownames
# rownames(counts_mat) <- counts_mat$Cell_Index
# counts_mat$Cell_Index <- NULL
# # check size of mat
# dim(counts_mat)
# # View top of mat
# counts_mat[1:5,1:5]
# ## tt - transpose
# counts_mat_tt <- rotate_df(counts_mat)
# 
# ## Ensure that cells and genes are in the correct position for Seurat
# head(colnames(counts_mat_tt))
# dim(counts_mat_tt)

##########################################################################
######################## Make SC obj #####################################
##########################################################################


# Seurat Filtered
seurat_filtered <- CreateSeuratObject(counts = count_mat_filtered_only_counts_tt)
## View
seurat_filtered

head(rownames(seurat_filtered))

# Seurat
# seurat <- CreateSeuratObject(counts = counts_mat_tt)
# ## View
# seurat
# head(rownames(seurat))

#
seurat_filtered <- AddMetaData(object = seurat_filtered, metadata = df_meta_all)
#
seurat_filtered@meta.data$Cell_Index <- as.character(df_meta_all$Cell_Index)
seurat_filtered@meta.data$Sample_Tag <- as.character(df_meta_all$Sample_Tag)
seurat_filtered@meta.data$Sample_Name <- as.character(df_meta_all$Sample_Name)
seurat_filtered@meta.data$Cell_Type_Experimental <- as.character(df_meta_all$Cell_Type_Experimental)
seurat_filtered@meta.data$Bar1 <- as.character(df_meta_all$Bar1)
seurat_filtered@meta.data$Bar2 <- as.character(df_meta_all$Bar2)
#
head(seurat_filtered)

##########################################################################
################################ QC ######################################
##########################################################################


# nFeature_RNA is the number of genes detected in each cell. nCount_RNA is the total number of molecules detected within a cell

# Visualize QC metrics as a violin plot
VlnPlot(seurat_filtered, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)

# scatter
plot <- FeatureScatter(seurat_filtered, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot

seurat_filtered_lower <- subset(seurat_filtered, subset = nFeature_RNA < 50)
VlnPlot(seurat_filtered_lower, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)

seurat_filtered_Upper <- subset(seurat_filtered, subset = nFeature_RNA > 50)
VlnPlot(seurat_filtered_Upper, features = c("nFeature_RNA", "nCount_RNA"), ncol = 2)

dim(seurat_filtered_lower)

dim(seurat_filtered_Upper)

# # 1 - SMK Calls file
# smk_calls <- read.csv(file = "C:/Users/Cathal/Documents/Nova_30k/Nova-LMOfilteredOUT-30kset-_Sample_Tag_Calls.csv", skip=7)
# head(smk_calls)
# 
# ####################################
# # 2 - Immune Classification cell type file
# Immune_class <- read.csv(file = "C:/Users/Cathal/Documents/Nova_30k/Combined_Nova-LMOfilteredOUT-30kset-_cell_type_experimental.csv")
# head(Immune_class)
# 
# #### merge SMK calls file and Immune class file
# df_meta_SevenBridges <- full_join(smk_calls, Immune_class, by = "Cell_Index")
# head(df_meta_SevenBridges)
# 
# ## Read in LMO barTable
# barTable_all <- read.csv("C:/Users/Cathal/Downloads/bar_table_full.csv")
# head(barTable_all)
# # reduce
# barTable_all <- select(barTable_all,-contains("nUMI"))
# barTable_all <- select(barTable_all,-contains("label"))
# 
# head(barTable_all)
# 
# ## Join sev bridges data with LMO barTable
# 
# dim(df_meta_SevenBridges)
# dim(barTable_all)
# ##
# df_meta_all <- left_join(df_meta_SevenBridges, barTable_all, by="Cell_Index")
# 
# head(df_meta_all)
# 
# # meta-data reduced
# 
# ## 
# # df_meta <- read.csv("C:/Users/Cathal/Desktop/Nova_seq/Nova_Seq_reduced/df_metadata_reduced.csv")[ ,c('Cell_Index', 'Sample_Tag', 'Sample_Name', 'Cell_Type_Experimental')]
# # 
# # head(df_meta)
# 
# # Add metadata to seurat object
# seurat <- AddMetaData(object = seurat, metadata = df_meta_all)
# ## as.character()
# seurat@meta.data$Cell_Type_Experimental <- as.character(df_meta_all$Cell_Type_Experimental)
# seurat@meta.data$Sample_Name <- as.character(df_meta_all$Sample_Name)
# seurat@meta.data$Sample_Tag <- as.character(df_meta_all$Sample_Tag)
# seurat@meta.data$Cell_Index <- as.character(df_meta_all$Cell_Index)
# seurat@meta.data$Bar1 <- as.character(df_meta_all$Bar1)
# seurat@meta.data$Bar2 <- as.character(df_meta_all$Bar2)
# # 
# head(seurat)

##########################################################################
####################### Seurat Normalise #################################
##########################################################################

seurat_filtered <- NormalizeData(object = seurat_filtered, normalization.method = "LogNormalize")
seurat_filtered_lower <- NormalizeData(object = seurat_filtered_lower, normalization.method = "LogNormalize")
seurat_filtered_Upper <- NormalizeData(object = seurat_filtered_Upper, normalization.method = "LogNormalize")

## Lower

########### Set features as Abseqs for Clustering
Abseqs <- rownames(x = seurat_filtered[1:30])
Abseqs

#all_genes <- rownames(seurat)

### Scale
seurat_filtered <- ScaleData(seurat_filtered)
### Run PCA with no of PC'S not specified, then check the Elbow plot. Re-run with the inflection number.
##seurat <- RunPCA(seurat, features = Abseqs)
#ElbowPlot(seurat)

### PCA
seurat_filtered <- RunPCA(seurat_filtered, features = Abseqs, npcs = 8)
### Find Neighbours
seurat_filtered <- FindNeighbors(seurat_filtered, dims = 1:8, k.param = 8)
### FindClusters
seurat_filtered <- FindClusters(seurat_filtered, resolution = 0.5)

### Run TSNE
seurat_filtered <- RunTSNE(seurat_filtered, features = Abseqs, check_duplicates = FALSE)
### Run UMAP
seurat_filtered <- RunUMAP(seurat_filtered, dims = 1:8)

## examine clustering results

print(seurat_filtered[["pca"]], dims = 1:5, nfeatures = 5)

VizDimLoadings(seurat_filtered, dims = 1, reduction = "pca")

VizDimLoadings(seurat_filtered, dims = 2, reduction = "pca")

seurat_filtered

seurat_filtered[["pca"]]

seurat_filtered[["tsne"]]

head(Loadings(seurat_filtered, reduction = "pca")[])

DimPlot(seurat_filtered, reduction = 'tsne') + ggtitle("tSNE filtered")

DimPlot(seurat_filtered, reduction = 'umap') + ggtitle("UMAP filtered")

## Lower

########### Set features as Abseqs for Clustering
Abseqs <- rownames(x = seurat_filtered_lower[1:30])
Abseqs

#all_genes <- rownames(seurat)

### Scale
seurat_filtered_lower <- ScaleData(seurat_filtered_lower)
### Run PCA with no of PC'S not specified, then check the Elbow plot. Re-run with the inflection number.
##seurat <- RunPCA(seurat, features = Abseqs)
#ElbowPlot(seurat)

### PCA
seurat_filtered_lower <- RunPCA(seurat_filtered_lower, features = Abseqs, npcs = 8)
### Find Neighbours
seurat_filtered_lower <- FindNeighbors(seurat_filtered_lower, dims = 1:8, k.param = 8)
### FindClusters
seurat_filtered_lower <- FindClusters(seurat_filtered_lower, resolution = 0.5)

### Run TSNE
seurat_filtered_lower <- RunTSNE(seurat_filtered_lower, features = Abseqs, check_duplicates = FALSE)
### Run UMAP
seurat_filtered_lower <- RunUMAP(seurat_filtered_lower, dims = 1:8)

## examine clustering results

print(seurat_filtered_lower[["pca"]], dims = 1:5, nfeatures = 5)

VizDimLoadings(seurat_filtered_lower, dims = 1, reduction = "pca")

VizDimLoadings(seurat_filtered_lower, dims = 2, reduction = "pca")

seurat_filtered_lower

seurat_filtered_lower[["pca"]]

seurat_filtered_lower[["tsne"]]

head(Loadings(seurat_filtered_lower, reduction = "pca")[])

DimPlot(seurat_filtered_lower, reduction = 'tsne') + ggtitle("tSNE Lower group")

DimPlot(seurat_filtered_lower, reduction = 'umap') + ggtitle("UMAP Lower group")

##########################################################################
####################### Seurat Clustering ################################
##########################################################################

## Upper

########### Set features as Abseqs for Clustering
Abseqs <- rownames(x = seurat_filtered_Upper[1:30])
Abseqs

#all_genes <- rownames(seurat)

### Scale
seurat_filtered_Upper <- ScaleData(seurat_filtered_Upper)
### Run PCA with no of PC'S not specified, then check the Elbow plot. Re-run with the inflection number.
##seurat <- RunPCA(seurat, features = Abseqs)
#ElbowPlot(seurat)
### PCA
#seurat <- RunPCA(seurat, features = Abseqs)
seurat_filtered_Upper <- RunPCA(seurat_filtered_Upper, npcs = 8, features = Abseqs)
### Find Neighbours
#seurat <- FindNeighbors(seurat, dims = 1:8)
seurat_filtered_Upper <- FindNeighbors(seurat_filtered_Upper, dims = 1:8, k.param = 8)
### FindClusters
seurat_filtered_Upper <- FindClusters(seurat_filtered_Upper, resolution = 0.5)

### Run TSNE
seurat_filtered_Upper <- RunTSNE(seurat_filtered_Upper, features = Abseqs, check_duplicates = FALSE)
### Run UMAP
seurat_filtered_Upper <- RunUMAP(seurat_filtered_Upper, dims = 1:8)

## examine clustering results
print(seurat_filtered_Upper[["pca"]], dims = 1:5, nfeatures = 5)

VizDimLoadings(seurat_filtered_Upper, dims = 1, reduction = "pca")

VizDimLoadings(seurat_filtered_Upper, dims = 2, reduction = "pca")

seurat_filtered_Upper

seurat_filtered_Upper[["pca"]]

seurat_filtered_Upper[["tsne"]]

head(Loadings(seurat_filtered_Upper, reduction = "pca")[])

DimPlot(seurat_filtered_Upper, reduction = 'tsne') + ggtitle("tSNE Upper group")

DimPlot(seurat_filtered_Upper, reduction = 'umap') + ggtitle("UMAP Upper group")

##########################################################################
############################ Dotplots ####################################
##########################################################################

########## Dotplot

# RNA list
rna_list <- c('CD3E',
              'CD4',
              'CD8A',
              'MS4A1',
              'CD14',
              'FCGR3A',
              'ITGAM',
              'CD5',
              'FCER2',
              'IL2RA',
              'CD69',
              'CD22',
              'CD79B')


## seurat filtered
## Clusters Idents
DotPlot(seurat_filtered, features = rna_list) + RotatedAxis() + ggtitle("Dotplot - filtered all. RNA list")

## calls file Idents
DotPlot(seurat_filtered, features = Abseqs, group.by = 'Cell_Type_Experimental') + RotatedAxis() + ggtitle("Dotplot - filtered all. Abseqs")

## Lower
## Clusters Idents
DotPlot(seurat_filtered_lower, features = rna_list) + RotatedAxis() + ggtitle("Dotplot - Lower. RNA list")

## calls file Idents
DotPlot(seurat_filtered_lower, features = Abseqs, group.by = 'Cell_Type_Experimental') + RotatedAxis() + ggtitle("Dotplot - Lower. Abseqs")

## Upper
## Clusters Idents
DotPlot(seurat_filtered_Upper, features = rna_list) + RotatedAxis() + ggtitle("Dotplot - Upper. RNA list")

## calls file Idents
DotPlot(seurat_filtered_Upper, features = Abseqs, group.by = 'Cell_Type_Experimental') + RotatedAxis() + ggtitle("Dotplot - Upper. Abseqs")

##########################################################################
####################### Subset SMK & LMO's and plot ######################
##########################################################################

### P1
# cond 1
LMO_1_cells_P1 <- WhichCells(seurat_filtered, expression = Bar1 >= 1 & Sample_Name == 'P1')
# cond 2
LMO_2_cells_P1 <- WhichCells(seurat_filtered, expression = Bar2 >= 1 & Sample_Name == 'P1')


DimPlot(seurat_filtered, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P1, LMO_2_cells_P1),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 1 -- SMK2") + theme_bw() + theme(legend.position = 'top')

# 
# 
# ### P2 
# # cond 1
LMO_1_cells_P2 <- WhichCells(seurat_filtered, expression = Bar1 >= 1 & Sample_Name == 'P2')
# cond 2
LMO_2_cells_P2 <- WhichCells(seurat_filtered, expression = Bar2 >= 1 & Sample_Name == 'P2')

DimPlot(seurat_filtered, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P2, LMO_2_cells_P2),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 2 -- SMK3") + theme_bw() + theme(legend.position = 'top')

# 
# 
# 
# 
# ### P3 
# # cond 1
LMO_1_cells_P3 <- WhichCells(seurat_filtered, expression = Bar1 >= 1 & Sample_Name == 'P3')
# cond 2
LMO_2_cells_P3 <- WhichCells(seurat_filtered, expression = Bar2 >= 1 & Sample_Name == 'P3')

DimPlot(seurat_filtered, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P3, LMO_2_cells_P3),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 3 -- SMK4") + theme_bw() + theme(legend.position = 'top')

# 
# 
# 
# 
# ### P4 
# # cond 1
LMO_1_cells_P4 <- WhichCells(seurat_filtered, expression = Bar1 >= 1 & Sample_Name == 'P4')
# cond 2
LMO_2_cells_P4 <- WhichCells(seurat_filtered, expression = Bar2 >= 1 & Sample_Name == 'P4')

DimPlot(seurat_filtered, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P4, LMO_2_cells_P4),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 4 -- SMK5") + theme_bw() + theme(legend.position = 'top')

# 
# 
# ### P5 
# # cond 1
LMO_1_cells_P5 <- WhichCells(seurat_filtered, expression = Bar1 >= 1 & Sample_Name == 'P5')
# cond 2
LMO_2_cells_P5 <- WhichCells(seurat_filtered, expression = Bar2 >= 1 & Sample_Name == 'P5')

DimPlot(seurat_filtered, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P5, LMO_2_cells_P5),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 5 -- SMK6") + theme_bw() + theme(legend.position = 'top')

# by Calls file
DimPlot(seurat_filtered, reduction = 'tsne', group.by = 'Cell_Type_Experimental')

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

##########################################################################
####################### Subset SMK & LMO's and plot ######################
##########################################################################

### P1
# cond 1
LMO_1_cells_P1 <- WhichCells(seurat_filtered_lower, expression = Bar1 >= 1 & Sample_Name == 'P1')
# cond 2
LMO_2_cells_P1 <- WhichCells(seurat_filtered_lower, expression = Bar2 >= 1 & Sample_Name == 'P1')


DimPlot(seurat_filtered_lower, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P1, LMO_2_cells_P1),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 1 -- SMK2") + theme_bw() + theme(legend.position = 'top')

# 
# 
# ### P2 
# # cond 1
LMO_1_cells_P2 <- WhichCells(seurat_filtered_lower, expression = Bar1 >= 1 & Sample_Name == 'P2')
# cond 2
LMO_2_cells_P2 <- WhichCells(seurat_filtered_lower, expression = Bar2 >= 1 & Sample_Name == 'P2')

DimPlot(seurat_filtered_lower, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P2, LMO_2_cells_P2),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 2 -- SMK3") + theme_bw() + theme(legend.position = 'top')

# 
# 
# 
# 
# ### P3 
# # cond 1
LMO_1_cells_P3 <- WhichCells(seurat_filtered_lower, expression = Bar1 >= 1 & Sample_Name == 'P3')
# cond 2
LMO_2_cells_P3 <- WhichCells(seurat_filtered_lower, expression = Bar2 >= 1 & Sample_Name == 'P3')

DimPlot(seurat_filtered_lower, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P3, LMO_2_cells_P3),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 3 -- SMK4") + theme_bw() + theme(legend.position = 'top')

# 
# 
# 
# 
# ### P4 
# # cond 1
LMO_1_cells_P4 <- WhichCells(seurat_filtered_lower, expression = Bar1 >= 1 & Sample_Name == 'P4')
# cond 2
LMO_2_cells_P4 <- WhichCells(seurat_filtered_lower, expression = Bar2 >= 1 & Sample_Name == 'P4')

DimPlot(seurat_filtered_lower, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P4, LMO_2_cells_P4),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 4 -- SMK5") + theme_bw() + theme(legend.position = 'top')

# 
# 
# ### P5 
# # cond 1
LMO_1_cells_P5 <- WhichCells(seurat_filtered_lower, expression = Bar1 >= 1 & Sample_Name == 'P5')
# cond 2
LMO_2_cells_P5 <- WhichCells(seurat_filtered_lower, expression = Bar2 >= 1 & Sample_Name == 'P5')

DimPlot(seurat_filtered_lower, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P5, LMO_2_cells_P5),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 5 -- SMK6") + theme_bw() + theme(legend.position = 'top')

# by Calls file
DimPlot(seurat_filtered_lower, reduction = 'tsne', group.by = 'Cell_Type_Experimental')

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

##########################################################################
####################### Subset SMK & LMO's and plot ######################
##########################################################################

### P1
# cond 1
LMO_1_cells_P1 <- WhichCells(seurat_filtered_Upper, expression = Bar1 >= 1 & Sample_Name == 'P1')
# cond 2
LMO_2_cells_P1 <- WhichCells(seurat_filtered_Upper, expression = Bar2 >= 1 & Sample_Name == 'P1')


DimPlot(seurat_filtered_Upper, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P1, LMO_2_cells_P1),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 1 -- SMK2") + theme_bw() + theme(legend.position = 'top')

# 
# 
# ### P2 
# # cond 1
LMO_1_cells_P2 <- WhichCells(seurat_filtered_Upper, expression = Bar1 >= 1 & Sample_Name == 'P2')
# cond 2
LMO_2_cells_P2 <- WhichCells(seurat_filtered_Upper, expression = Bar2 >= 1 & Sample_Name == 'P2')

DimPlot(seurat_filtered_Upper, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P2, LMO_2_cells_P2),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 2 -- SMK3") + theme_bw() + theme(legend.position = 'top')

# 
# 
# 
# 
# ### P3 
# # cond 1
LMO_1_cells_P3 <- WhichCells(seurat_filtered_Upper, expression = Bar1 >= 1 & Sample_Name == 'P3')
# cond 2
LMO_2_cells_P3 <- WhichCells(seurat_filtered_Upper, expression = Bar2 >= 1 & Sample_Name == 'P3')

DimPlot(seurat_filtered_Upper, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P3, LMO_2_cells_P3),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 3 -- SMK4") + theme_bw() + theme(legend.position = 'top')

# 
# 
# 
# 
# ### P4 
# # cond 1
LMO_1_cells_P4 <- WhichCells(seurat_filtered_Upper, expression = Bar1 >= 1 & Sample_Name == 'P4')
# cond 2
LMO_2_cells_P4 <- WhichCells(seurat_filtered_Upper, expression = Bar2 >= 1 & Sample_Name == 'P4')

DimPlot(seurat_filtered_Upper, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P4, LMO_2_cells_P4),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 4 -- SMK5") + theme_bw() + theme(legend.position = 'top')

# 
# 
# ### P5 
# # cond 1
LMO_1_cells_P5 <- WhichCells(seurat_filtered_Upper, expression = Bar1 >= 1 & Sample_Name == 'P5')
# cond 2
LMO_2_cells_P5 <- WhichCells(seurat_filtered_Upper, expression = Bar2 >= 1 & Sample_Name == 'P5')

DimPlot(seurat_filtered_Upper, reduction = "tsne",
        cells.highlight = list(LMO_1_cells_P5, LMO_2_cells_P5),
        cols.highlight = c("darkblue", "darkorange"), cols = "grey", sizes.highlight = 0.9) +
    ggtitle("Patient 5 -- SMK6") + theme_bw() + theme(legend.position = 'top')

# by Calls file
DimPlot(seurat_filtered_Upper, reduction = 'tsne', group.by = 'Cell_Type_Experimental')

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

# seu filtered

# ## Pre-process Seurat object (standard) --------------------------------------------------------------------------------------
# #seu_kidney <- CreateSeuratObject(kidney.data)
# #seu_kidney <- NormalizeData(seu_kidney)
# seurat_filtered2 <- FindVariableFeatures(seurat_filtered, selection.method = "vst", nfeatures = 200)
# #seu_kidney <- ScaleData(seu_kidney)
# #seu_kidney <- RunPCA(seu_kidney)
# #seu_kidney <- RunUMAP(seu_kidney, dims = 1:10)
# 
# 
# ## Pre-process Seurat object (sctransform) -----------------------------------------------------------------------------------
# #seu_kidney <- CreateSeuratObject(kidney.data)
# seurat_filtered2 <- SCTransform(seurat_filtered2)
# #seu_kidney <- RunPCA(seu_kidney)
# #seu_kidney <- RunUMAP(seu_kidney, dims = 1:10)
# 
# 
# ## pK Identification (no ground-truth) ---------------------------------------------------------------------------------------
# sweep.res.list_filtered <- paramSweep_v3(seurat_filtered2, PCs = 1:8, sct = FALSE)
# sweep.stats_filtered <- summarizeSweep(sweep.res.list_filtered, GT = FALSE)
# bcmvn_filtered <- find.pK(sweep.stats_filtered)
# 
# 
# 
# ## pK Identification (ground-truth) ------------------------------------------------------------------------------------------
# sweep.res.list_filtered <- paramSweep_v3(seurat_filtered2, PCs = 1:8, sct = FALSE)
# gt.calls <- seurat_filtered2@meta.data[rownames(sweep.res.list_filtered[[1]]), "GT"].   ## GT is a vector containing "Singlet" and "Doublet" calls recorded using sample multiplexing classification and/or in silico geneotyping results 
# sweep.stats_filtered <- summarizeSweep(sweep.res.list_filtered, GT = TRUE, GT.calls = gt.calls)
# bcmvn_filtered <- find.pK(sweep.stats_filtered)
# 
# 
# ## Homotypic Doublet Proportion Estimate -------------------------------------------------------------------------------------
# homotypic.prop <- modelHomotypic(seurat_filtered2@meta.data$seurat_clusters)           ## ex: annotations <- seu_kidney@meta.data$ClusteringResults
# nExp_poi <- round(0.075*nrow(seurat_filtered2@meta.data))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
# nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
# 
# ## Run DoubletFinder with varying classification stringencies ----------------------------------------------------------------
# seurat_filtered2 <- doubletFinder_v3(seurat_filtered2, PCs = 1:8, pN = 0.25, pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
# seurat_filtered2 <- doubletFinder_v3(seurat_filtered2, PCs = 1:8, pN = 0.25, pK = 0.09, nExp = nExp_poi.adj, reuse.pANN = "pANN_0.25_0.09_500", sct = FALSE)
# 
# 
# ## plot
# 
# seurat_filtered2 <- doubletFinder_v3(seurat_filtered2, PCs = 1:8, pN = 0.25, pK = 0.09, nExp = nExp_poi.adj, reuse.pANN = "DF.classifications_0.25_0.09_500", sct = FALSE)
# 
# seurat_filtered2@meta.data[,"CellTypes_DF"] <- seurat_filtered2@meta.data$CT.Park
# 
# seurat_filtered2@meta.data$CellTypes_DF[which(seurat_filtered2@meta.data$DF.classifications_0.25_0.09_500 == "Doublet")] <- "Doublet"
# 
# DimPlot(seurat_filtered2, group.by="CellTypes_DF", reduction="tsne", pt.size=0.5, order=c("Coll.Duct.TC","Doublet"), cols=c("#66C2A5","#FFD92F","#8DA0CB","#A6D854","#E78AC3","#B3B3B3","#E5C494","black","#FC8D62"))