library(Seurat)
library(SeuratWrappers)
library(openxlsx)

wdir <- "Birocchi_SciTraMed2022/results"
full_min <- readRDS(paste(wdir, "Full", "Full_minimal.rds", sep = "/"))

# Compute mito % 
full_min[["percent.mt"]] <- PercentageFeatureSet(full_min, pattern = "^mt-")

png(filename = paste(wdir, "Full", "01-fastMNN_2", paste("Full", "VlnQCplots.png", sep = "_"), sep = "/"),
    units = "in", width = 12, height = 9, res = 100)
VlnPlot(full_min, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.2)
dev.off()

# Filter cells
full_min <- subset(x = full_min, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 10)
# Normalize
full_min <- NormalizeData(full_min)
full_min <- FindVariableFeatures(full_min, selection.method = "vst",
                                 nfeatures = ceiling(dim(full_min)[1] * 0.2),
                                 verbose = T)
# Run FastMNN
full_min <- RunFastMNN(object.list = SplitObject(full_min, split.by = "RNA_Group"))
full_min <- RunUMAP(full_min, reduction = "mnn", dims = 1:50)
full_min <- FindNeighbors(full_min, reduction = "mnn", dims = 1:50)
for (res in c(1.0)) {
    full_min <- FindClusters(full_min, resolution = res)
}

png(filename = paste(wdir, "Full", "01-fastMNN_2", "Full_RNA_Group_res1.2_umap_plot.png", sep = "/"),
    width = 14, height = 8, units = "in", res = 300)
DimPlot(full_min,
        group.by = c("RNA_Group", "RNA_snn_res.1.2"),
        ncol = 2, label = T)
dev.off()

png(filename = paste(wdir, "Full", "01-fastMNN_2", "Full_RNA_Group_res1.8_umap_plot.png", sep = "/"),
    width = 14, height = 8, units = "in", res = 300)
DimPlot(full_min,
        group.by = c("RNA_Group", "RNA_snn_res.1.8"),
        ncol = 2, label = T)
dev.off()

# Compute Markers
df_markers <- list()
for (md in c("RNA_snn_res.1")) { # "RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"
    Idents(full_min) <- md
    png(paste(wdir, "Full", "01-fastMNN_2", paste0("Full_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = full_min, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = full_min, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "Full", "01-fastMNN_2", paste("Full_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
full_min@misc[["markers"]] <- df_markers
# Save final object
saveRDS(object = full_min, file = paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))

##########################
### Filter Macrophages ###
##########################
Idents(full_min) <- "RNA_snn_res.1"
full_min_macro <- subset(x = full_min, idents = c(0, 1, 4, 5, 12, 14, 18, 19, 22))

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "Full_Macrophages_subset_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap", group.by = "RNA_snn_res.1", label = T)
dev.off()

full_min_macro <- FindVariableFeatures(full_min_macro, selection.method = "vst",
                                       nfeatures = ceiling(dim(full_min_macro)[1] * 0.2),
                                       verbose = T)
full_min_macro <- RunFastMNN(object.list = SplitObject(full_min_macro, split.by = "RNA_Group"))
full_min_macro <- RunUMAP(full_min_macro, reduction = "mnn", dims = 1:30)
full_min_macro <- FindNeighbors(full_min_macro, reduction = "mnn", dims = 1:30)
for (res in c(0.6, 1.2, 1.8)) {
    full_min_macro <- FindClusters(full_min_macro, resolution = res)
}

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
    width = 14, height = 12, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
        group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"),
        ncol = 2, label = T)
dev.off()

df_markers <- list()
for (md in c("RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
    Idents(full_min_macro) <- md
    png(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", paste0("MacrophagesFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = full_min_macro, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = full_min_macro, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2",
                              paste("MacrophagesFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
full_min_macro@misc[["markers_full"]] <- df_markers

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
    width = 14, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap", group.by = c("RNA_Group", "orig.ident"), ncol = 2, label = T)
dev.off()

# cell-cycle scoring
cc.genes <- readRDS(file = "regev_lab_cell_cycle_mouse.rds")
s.genes <- cc.genes[["s.genes"]]
g2m.genes <- cc.genes[["g2m.genes"]]
full_min_macro@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
full_min_macro <- CellCycleScoring(object = full_min_macro, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
full_min_macro@meta.data$CC.Difference <- full_min_macro@meta.data$S.Score - full_min_macro@meta.data$G2M.Score

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_CCPhase_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
        group.by = c("Phase"),
        ncol = 2)
dev.off()

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
        group.by = c("Phase"), split.by = "RNA_Group",
        ncol = 3)
dev.off()

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_macro, reduction = "umap",
        group.by = c("Phase"), split.by = "Phase",
        ncol = 3)
dev.off()

# Save final object
saveRDS(object = full_min_macro,
        file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))

## Filter Monocytes
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
Idents(macro_obj) <- "RNA_snn_res.0.6"
monoc <- subset(macro_obj, idents = c(7, 10))

monoc <- FindVariableFeatures(monoc, selection.method = "vst", nfeatures = ceiling(dim(monoc)[1] * 0.2), verbose = T)
monoc <- RunFastMNN(object.list = SplitObject(monoc, split.by = "RNA_Group"))
monoc <- RunUMAP(monoc, reduction = "mnn", dims = 1:20)
monoc <- FindNeighbors(monoc, reduction = "mnn", dims = 1:20)
for (res in c(0.6, 1.2, 1.8)) {
    monoc <- FindClusters(monoc, resolution = res)
}

png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
    width = 14, height = 12, units = "in", res = 200)
DimPlot(monoc, reduction = "umap",
        group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"),
        ncol = 2, label = T)
dev.off()

df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
    Idents(monoc) <- md
    png(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", paste0("MonocytesFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = monoc, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = monoc, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2",
                              paste("MonocytesFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
monoc@misc[["markers_full"]] <- df_markers

png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
    width = 14, height = 8, units = "in", res = 300)
DimPlot(monoc, reduction = "umap",
        group.by = c("RNA_Group", "orig.ident"),
        ncol = 2, label = T)
dev.off()

# Save object
saveRDS(object = monoc, file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_fastMNN_final.rds", sep = "/"))

# Filter out Clu 3
Idents(monoc) <- "RNA_snn_res.0.6"
monoc_filt <- subset(monoc, idents = c(0, 1, 2))

monoc_filt <- FindVariableFeatures(monoc_filt,
                                   selection.method = "vst", nfeatures = ceiling(dim(monoc_filt)[1] * 0.2), verbose = T)
monoc_filt <- RunFastMNN(object.list = SplitObject(monoc_filt, split.by = "RNA_Group"))
monoc_filt <- RunUMAP(monoc_filt, reduction = "mnn", dims = 1:20)
monoc_filt <- FindNeighbors(monoc_filt, reduction = "mnn", dims = 1:20)
for (res in c(0.6, 1.2, 1.8)) {
    monoc_filt <- FindClusters(monoc_filt, resolution = res)
}

png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
    width = 14, height = 12, units = "in", res = 200)
DimPlot(monoc_filt, reduction = "umap",
        group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8"),
        ncol = 2, label = T)
dev.off()

df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
    Idents(monoc_filt) <- md
    png(paste(wdir, "MonocytesFiltered", "01-fastMNN_2", paste0("MonocytesFiltered_No3_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = monoc_filt, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = monoc_filt, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2",
                              paste("MonocytesFiltered_No3_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
monoc_filt@misc[["markers_full"]] <- df_markers

png(filename = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
    width = 14, height = 8, units = "in", res = 300)
DimPlot(monoc_filt, reduction = "umap",
        group.by = c("RNA_Group", "orig.ident"),
        ncol = 2, label = T)
dev.off()

saveRDS(object = monoc_filt,
        file = paste(wdir, "MonocytesFiltered", "01-fastMNN_2", "MonocytesFiltered_No3_fastMNN_final.rds", sep = "/"))

######################
### Filter T cells ###
######################
Idents(full_min) <- "RNA_snn_res.1"
full_min_tcell <- subset(x = full_min, idents = c(2, 3, 6, 7, 11, 17))

png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Full_Tcell_subset_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_tcell, reduction = "umap", group.by = "RNA_snn_res.1", label = T)
dev.off()

full_min_tcell <- FindVariableFeatures(full_min_tcell, selection.method = "vst",
                                       nfeatures = ceiling(dim(full_min_tcell)[1] * 0.2),
                                       verbose = T)
full_min_tcell <- RunFastMNN(object.list = SplitObject(full_min_tcell, split.by = "RNA_Group"))
full_min_tcell <- RunUMAP(full_min_tcell, reduction = "mnn", dims = 1:30)
full_min_tcell <- FindNeighbors(full_min_tcell, reduction = "mnn", dims = 1:30)
for (res in c(0.6, 1.2, 1.8)) {
    full_min_tcell <- FindClusters(full_min_tcell, resolution = res)
}

png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_Cd4_Cd8_umap_plot.png", sep = "/"),
    width = 12, height = 8, units = "in", res = 300)
FeaturePlot(full_min_tcell, reduction = "umap", features = c("Cd4", "Cd8a"))
dev.off()

for (md in c("RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8", "RNA_snn_res.2.2", "RNA_snn_res.2.4")) {
    Idents(full_min_tcell) <- md
    png(paste(wdir, "TcellFiltered", "01-fastMNN_2", paste0("TcellFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = full_min_tcell, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    png(paste(wdir, "TcellFiltered", "01-fastMNN_2", paste0("TcellFiltered_fastMNN_", md, "_vln_plot.png"), sep = "/"),
        width = 12, height = 8, units = "in", res = 300)
    print(VlnPlot(object = full_min_tcell, features = c("Cd4", "Cd8a"), pt.size = 0.2))
    dev.off()
}

df_markers <- list()
for (md in c("RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
    Idents(full_min_tcell) <- md
    png(paste(wdir, "TcellFiltered", "01-fastMNN_2", paste0("TcellFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = full_min_macro, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = full_min_tcell, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "TcellFiltered", "01-fastMNN_2",
                              paste("TcellFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
full_min_tcell@misc[["markers_full"]] <- df_markers

saveRDS(object = full_min_tcell,
        file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "TcellFiltered_fastMNN_final.rds", sep = "/"))


### Subset of Cd4 and Cd8
full_min_tcell <- FindClusters(full_min_tcell, resolution = 8)

## Cd4
Idents(full_min_tcell) <- "RNA_snn_res.8"
cd4_sub <- subset(x = full_min_tcell, idents = c(0, 1, 3, 4, 5, 7, 9, 10, 13, 14, 16, 17, 18, 21, 22, 24, 27, 30, 31, 33, 34, 36, 38, 43, 44, 45, 46, 47, 48, 50, 51, 52, 53, 55, 56, 57))
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub_fastMNN_Cd4_Cd8_umap_plot.png", sep = "/"),
    width = 12, height = 8, units = "in", res = 300)
FeaturePlot(cd4_sub, reduction = "umap", features = c("Cd4", "Cd8a"))
dev.off()
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub_fastMNN_umap_plot.png", sep = "/"),
    width = 12, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap", pt.size = 1.5, label = T)
dev.off()

cd4_sub <- NormalizeData(cd4_sub)
cd4_sub <- FindVariableFeatures(cd4_sub, selection.method = "vst",
                                nfeatures = ceiling(dim(cd4_sub)[1] * 0.2),
                                verbose = T)
cd4_sub <- RunFastMNN(object.list = SplitObject(cd4_sub, split.by = "RNA_Group"))
cd4_sub <- RunUMAP(cd4_sub, reduction = "mnn", dims = 1:15)
cd4_sub <- FindNeighbors(cd4_sub, reduction = "mnn", dims = 1:15)
for (res in c(0.6, 1.2, 1.8)) {
    cd4_sub <- FindClusters(cd4_sub, resolution = res)
}

df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
    Idents(cd4_sub) <- md
    png(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", paste0("Cd4sub_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = cd4_sub, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = cd4_sub, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub",
                              paste("Cd4sub_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
cd4_sub@misc[["markers"]] <- df_markers

# cell-cycle scoring
cd4_sub@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
cd4_sub <- CellCycleScoring(object = cd4_sub, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
cd4_sub@meta.data$CC.Difference <- cd4_sub@meta.data$S.Score - cd4_sub@meta.data$G2M.Score

png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_CCPhase_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap",
        group.by = c("Phase"),
        ncol = 2)
dev.off()

png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap",
        group.by = c("Phase"), split.by = "RNA_Group",
        ncol = 3)
dev.off()

png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(cd4_sub, reduction = "umap",
        group.by = c("Phase"), split.by = "Phase",
        ncol = 3)
dev.off()

# Save object
saveRDS(object = cd4_sub,
        file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd4sub", "Cd4sub_fastMNN_final.rds", sep = "/"))


## Cd8
cd8_sub <- subset(x = full_min_tcell, idents = c(2, 6, 8, 11, 12, 15, 19, 20, 23, 25, 26, 28, 29, 32, 35, 37, 39, 40, 41, 42, 49, 54, 58, 59, 60, 61))
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub_fastMNN_Cd4_Cd8_umap_plot.png", sep = "/"),
    width = 12, height = 8, units = "in", res = 300)
FeaturePlot(cd8_sub, reduction = "umap", features = c("Cd4", "Cd8a"), order = T)
dev.off()
png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub_fastMNN_umap_plot.png", sep = "/"),
    width = 12, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap", pt.size = 1.5, label = T, order= T)
dev.off()

cd8_sub <- NormalizeData(cd8_sub)
cd8_sub <- FindVariableFeatures(cd8_sub, selection.method = "vst",
                                nfeatures = ceiling(dim(cd8_sub)[1] * 0.2),
                                verbose = T)
cd8_sub <- RunFastMNN(object.list = SplitObject(cd8_sub, split.by = "RNA_Group"))
cd8_sub <- RunUMAP(cd8_sub, reduction = "mnn", dims = 1:15)
cd8_sub <- FindNeighbors(cd8_sub, reduction = "mnn", dims = 1:15)
for (res in c(0.6, 1.2, 1.8)) {
    cd8_sub <- FindClusters(cd8_sub, resolution = res)
}

df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
    Idents(cd8_sub) <- md
    png(paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", paste0("Cd8sub_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = cd8_sub, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = cd8_sub, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub",
                              paste("Cd8sub_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
cd8_sub@misc[["markers"]] <- df_markers

# cell-cycle scoring
cd8_sub@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
cd8_sub <- CellCycleScoring(object = cd8_sub, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
cd8_sub@meta.data$CC.Difference <- cd8_sub@meta.data$S.Score - cd8_sub@meta.data$G2M.Score

png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_CCPhase_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap",
        group.by = c("Phase"),
        ncol = 2)
dev.off()

png(filename = paste(wdir,"TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap",
        group.by = c("Phase"), split.by = "RNA_Group",
        ncol = 3)
dev.off()

png(filename = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(cd8_sub, reduction = "umap",
        group.by = c("Phase"), split.by = "Phase",
        ncol = 3)
dev.off()

# Save object
saveRDS(object = cd8_sub,
        file = paste(wdir, "TcellFiltered", "01-fastMNN_2", "Cd8sub", "Cd8sub_fastMNN_final.rds", sep = "/"))

##############################
### Filter Denditric cells ###
##############################
full_obj <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
Idents(full_obj) <- "RNA_snn_res.1.2"
DimPlot(object = full_obj, reduction = "umap", label = T)
full_min_dc <- subset(x = full_obj, idents = c(10, 11, 16, 22))

png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "Full_Dendritic_subset_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap", group.by = "RNA_snn_res.1.2", label = T)
dev.off()

full_min_dc <- NormalizeData(full_min_dc)
full_min_dc <- FindVariableFeatures(full_min_dc, selection.method = "vst",
                                    nfeatures = ceiling(dim(full_min_dc)[1] * 0.2), verbose = T)
full_min_dc <- RunFastMNN(object.list = SplitObject(full_min_dc, split.by = "RNA_Group"))
full_min_dc <- RunUMAP(full_min_dc, reduction = "mnn", dims = 1:15)
full_min_dc <- FindNeighbors(full_min_dc, reduction = "mnn", dims = 1:15)
for (res in c(0.6, 1.2, 1.8)) {
    full_min_dc <- FindClusters(full_min_dc, resolution = res)
}

df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2", "RNA_snn_res.1.8")) {
    Idents(full_min_dc) <- md
    png(paste(wdir,  "DendriticFilteredNEW", "01-fastMNN_2", paste0("DendriticFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = full_min_dc, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = full_min_dc, min.pct = 0.1, logfc.threshold = 0.25, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1e-06)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2",
                              paste("DendriticFilteredNEW_fastMNN", md, "markers.xlsx", sep = "_"), sep = "/"))
}
full_min_dc@misc[["markers"]] <- df_markers

# cell-cycle scoring
full_min_dc@misc[["cc.genes"]] <- c(s.genes, g2m.genes)
full_min_dc <- CellCycleScoring(object = full_min_dc, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
full_min_dc@meta.data$CC.Difference <- full_min_dc@meta.data$S.Score - full_min_dc@meta.data$G2M.Score

png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_CCPhase_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap",
        group.by = c("Phase"),
        ncol = 2)
dev.off()

png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_CCPhase_split_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap",
        group.by = c("Phase"), split.by = "RNA_Group",
        ncol = 3)
dev.off()

png(filename = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_CCPhase_splitCC_umap_plot.png", sep = "/"),
    width = 16, height = 8, units = "in", res = 300)
DimPlot(full_min_dc, reduction = "umap",
        group.by = c("Phase"), split.by = "Phase",
        ncol = 3)
dev.off()

# Save object
saveRDS(object = full_min_dc,
        file = paste(wdir, "DendriticFilteredNEW", "01-fastMNN_2", "DendriticFilteredNEW_fastMNN_final.rds", sep = "/"))



######################
### Filter B cells ###
######################
fmnn_full <- readRDS(paste(wdir, "Full", "01-fastMNN_2", "Full_fastMNN_final.rds", sep = "/"))
Idents(fmnn_full) <- "RNA_snn_res.0.6"
fmnn_bcells <- subset(fmnn_full, idents = c(12))

fmnn_bcells <- FindVariableFeatures(fmnn_bcells, selection.method = "vst", nfeatures = ceiling(dim(fmnn_bcells)[1] * 0.2), verbose = T)
fmnn_bcells <- RunFastMNN(object.list = SplitObject(fmnn_bcells, split.by = "RNA_Group"))
fmnn_bcells <- RunUMAP(fmnn_bcells, reduction = "mnn", dims = 1:20)
fmnn_bcells <- FindNeighbors(fmnn_bcells, reduction = "mnn", dims = 1:20)
for (res in c(0.6, 1.2)) {
    fmnn_bcells <- FindClusters(fmnn_bcells, resolution = res)
}

png(filename = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_Group_AllRes_umap_plot.png", sep = "/"),
    width = 14, height = 12, units = "in", res = 200)
DimPlot(fmnn_bcells, reduction = "umap",
        group.by = c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2"),
        ncol = 2, label = T)
dev.off()

df_markers <- list()
for (md in c("RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1.2")) {
    Idents(fmnn_bcells) <- md
    png(paste(wdir, "BcellFiltered", "01-fastMNN_2", paste0("BcellFiltered_fastMNN_", md, "_umap_plot.png"), sep = "/"),
        width = 10, height = 8, units = "in", res = 300)
    print(DimPlot(object = fmnn_bcells, reduction = "umap", pt.size = 1.5, label = T))
    dev.off()
    df_markers[[md]] <- FindAllMarkers(object = fmnn_bcells, min.pct = 0.1, logfc.threshold = 0, only.pos = FALSE,
                                       min.cells.group = 10, test.use = "wilcox", return.thresh = 1)
    wb_markers <- createWorkbook()
    for (clu in unique(df_markers[[md]]$cluster)) {
        df_markers_clu <- subset(df_markers[[md]], cluster == clu)
        addWorksheet(wb_markers, clu)
        writeData(wb_markers, sheet = clu, x = df_markers_clu)
    }
    saveWorkbook(wb = wb_markers, overwrite = T,
                 file = paste(wdir, "BcellFiltered", "01-fastMNN_2",
                              paste("BcellFiltered_fastMNN", md, "markers_full.xlsx", sep = "_"), sep = "/"))
}
fmnn_bcells@misc[["markers_full"]] <- df_markers

png(filename = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_Group_Orig_umap_plot.png", sep = "/"),
    width = 14, height = 8, units = "in", res = 300)
DimPlot(fmnn_bcells, reduction = "umap",
        group.by = c("RNA_Group", "orig.ident"),
        ncol = 2, label = T)
dev.off()

# Save object
saveRDS(object = fmnn_bcells, file = paste(wdir, "BcellFiltered", "01-fastMNN_2", "BcellFiltered_fastMNN_final.rds", sep = "/"))