# Run seurat to combine h5ad datasets

# load libraries

library(scater)
library(Seurat)
library(cowplot)
library(reticulate)
library(SeuratData)
library(SeuratDisk)
library(anndata)
library(optparse)
library(data.table)
 
#Rscript soynod.R --input1 ../Nodule_1_gex_only/outs/filtered_feature_bc_matrix/ --input2 ../Nodule_2_gex_only/outs/filtered_feature_bc_matrix/ --qc1 ./Nodule_1_3pct_contam.qc.h5ad --qc2 ./Nodule_2_3pct_contam.qc.h5ad

# Read cellranger input
cellranger.input1 <- Read10X(data.dir = "../Nodule_1_gex_only/outs/filtered_feature_bc_matrix/")
cellranger.input2 <- Read10X(data.dir = "../Nodule_2_gex_only/outs/filtered_feature_bc_matrix/")
cellranger.input3 <- Read10X(data.dir = "../Nodule_3_gex_only/outs/filtered_feature_bc_matrix/")

# Read QC data
qc.input1 <- read_h5ad("./Nodule_1_3pct_contam.qc.h5ad")
qc.input2 <- read_h5ad("./Nodule_2_3pct_contam.qc.h5ad")
qc.input3 <- read_h5ad("./Nodule_2_3pct_contam.qc.h5ad")

# load objects into seurat
cellranger.seurat1 <- CreateSeuratObject(counts = cellranger.input1, project = "Nodule1_cellranger")
cellranger.seurat2 <- CreateSeuratObject(counts = cellranger.input2, project = "Nodule2_cellranger")
cellranger.seurat3 <- CreateSeuratObject(counts = cellranger.input3, project = "Nodule3_cellranger")

# load qc objects into seurat
qc.seurat1 <- CreateSeuratObject(counts = t(qc.input1$X), meta.data = qc.input1$obs, project = "Nodule1_qc")
qc.seurat2 <- CreateSeuratObject(counts = t(qc.input2$X), meta.data = qc.input2$obs, project = "Nodule2_qc")
qc.seurat3 <- CreateSeuratObject(counts = t(qc.input3$X), meta.data = qc.input3$obs, project = "Nodule3_qc")

#str(cellranger.seurat1)
#str(cellranger.seurat2)
#str(qc.seurat1)
#str(qc.seurat2)

# Integrate objects
cellranger.combined <- merge(cellranger.seurat1, y = c(cellranger.seurat2, cellranger.seurat3), add.cell.ids = c("Nodule1_cellranger", "Nodule2_cellranger", "Nodule3_cellranger"), project = "Nodule_cellranger")
qc.combined <- merge(qc.seurat1, y = c(qc.seurat2, qc.seurat3), add.cell.ids = c("Nodule1_qc", "Nodule2_qc", "Nodule3_qc"), project = "Nodule_qc")
#str(wholeroot.combined.qc)

#str(SplitObject(wholeroot.combined.qc, split.by = "orig.ident"))

str(cellranger.combined)
str(qc.combined)

# split the datasets into a list of two seurat objects
cellranger.list <- SplitObject(cellranger.combined, split.by = "orig.ident")
qc.list <- SplitObject(qc.combined, split.by = "orig.ident")

# make new lists for SCtransform

cellranger.list.sc <- SplitObject(cellranger.combined, split.by = "orig.ident")
qc.list.sc <- SplitObject(qc.combined, split.by = "orig.ident")

# normalize and identify variable features for each dataset independently
cellranger.list <- lapply(X = cellranger.list, FUN = function(x) {
                            x <- NormalizeData(x)
                            x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
})

qc.list <- lapply(X = qc.list, FUN = function(x) {
                    x <- NormalizeData(x)
                    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
})

cellranger.list.sc <- lapply(X = cellranger.list, FUN = function(x) {
                            x <- SCTransform(x)
                            x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
})

qc.list.sc <- lapply(X = qc.list, FUN = function(x) {
                    x <- SCTransform(x)
                    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000)
})

# Standard Integration changes default features to 3000

cellranger.features <- SelectIntegrationFeatures(object.list = cellranger.list, nfeatures = 3000)
qc.features <- SelectIntegrationFeatures(object.list = qc.list, nfeatures = 3000)

cellranger.anchors <-  FindIntegrationAnchors(object.list = cellranger.list, anchor.features = cellranger.features)
qc.anchors <-  FindIntegrationAnchors(object.list = qc.list, anchor.features = qc.features)

cellranger.combined <- IntegrateData(anchorset = cellranger.anchors)
qc.combined <- IntegrateData(anchorset = qc.anchors)


# SCtransform integration
cellranger.list.sc <- lapply(X = cellranger.list.sc, FUN = SCTransform)
cellranger.features.sc <- SelectIntegrationFeatures(object.list = cellranger.list.sc, nfeatures = 3000)
cellranger.list.sc <- PrepSCTIntegration(object.list = cellranger.list.sc, anchor.features = cellranger.features.sc)
cellranger.sc.anchors <- FindIntegrationAnchors(object.list = cellranger.list.sc, normalization.method = "SCT",
                                                anchor.features = cellranger.features.sc)
cellranger.combined.sct <- IntegrateData(anchorset = cellranger.sc.anchors, normalization.method = "SCT")

qc.list.sc <- lapply(X = qc.list.sc, FUN = SCTransform)
qc.features.sc <- SelectIntegrationFeatures(object.list = qc.list.sc, nfeatures = 3000)
qc.list.sc <- PrepSCTIntegration(object.list = qc.list.sc, anchor.features = qc.features.sc)
qc.sc.anchors <- FindIntegrationAnchors(object.list = qc.list.sc, normalization.method = "SCT",
                                                anchor.features = qc.features.sc)
qc.combined.sct <- IntegrateData(anchorset = qc.sc.anchors, normalization.method = "SCT")

#cellranger.combined.sct <- RunPCA(cellranger.combined.sct, verbose = FALSE)
#cellranger.combined.sct <- RunUMAP(cellranger.combined.sct, reduction = "pca", dims = 1:30)

#str(cellranger.combined)
#str(qc.combined)

fwrite(x = list(qc.features.sc), file = "./qc.features.sc.txt")
fwrite(x = list(qc.features), file = "./qc.features.txt")
fwrite(x = list(cellranger.features.sc), file = "./cellranger.features.sc.txt")
fwrite(x = list(cellranger.features), file = "./cellranger.features.txt")

DefaultAssay(cellranger.combined) <- "integrated"
DefaultAssay(qc.combined) <- "integrated"
DefaultAssay(cellranger.combined.sct) <- "integrated"
DefaultAssay(qc.combined.sct) <- "integrated"

cellranger.combined <- ScaleData(cellranger.combined, verbose = FALSE)
cellranger.combined <- RunPCA(cellranger.combined, npcs = 30, verbose = FALSE)
cellranger.combined <- RunUMAP(cellranger.combined, reduction = "pca", dims = 1:30)
cellranger.combined <- FindNeighbors(cellranger.combined, reduction = "pca", dims = 1:30)
cellranger.combined <- FindClusters(cellranger.combined, resolution = 0.5)

qc.combined <- ScaleData(qc.combined, verbose = FALSE)
qc.combined <- RunPCA(qc.combined, npcs = 30, verbose = FALSE)
qc.combined <- RunUMAP(qc.combined, reduction = "pca", dims = 1:30)
qc.combined <- FindNeighbors(qc.combined, reduction = "pca", dims = 1:30)
qc.combined <- FindClusters(qc.combined, resolution = 0.5)

cellranger.combined.sct <- ScaleData(cellranger.combined.sct, verbose = FALSE)
cellranger.combined.sct <- RunPCA(cellranger.combined.sct, npcs = 30, verbose = FALSE)
cellranger.combined.sct <- RunUMAP(cellranger.combined.sct, reduction = "pca", dims = 1:30)
cellranger.combined.sct <- FindNeighbors(cellranger.combined.sct, reduction = "pca", dims = 1:30)
cellranger.combined.sct <- FindClusters(cellranger.combined.sct, resolution = 0.5)

qc.combined.sct <- ScaleData(qc.combined.sct, verbose = FALSE)
qc.combined.sct <- RunPCA(qc.combined.sct, npcs = 30, verbose = FALSE)
qc.combined.sct <- RunUMAP(qc.combined.sct, reduction = "pca", dims = 1:30)
qc.combined.sct <- FindNeighbors(qc.combined.sct, reduction = "pca", dims = 1:30)
qc.combined.sct <- FindClusters(qc.combined.sct, resolution = 0.5)

# write all the files

write.csv(qc.combined@reductions$umap@cell.embeddings, file = "./qc.strict.umap.csv", quote=F)
write.csv(qc.combined.sct@reductions$umap@cell.embeddings, file = "./qc.strict.sct.umap.csv", quote=F)
write.csv(cellranger.combined@reductions$umap@cell.embeddings, file = "./cellranger.umap.csv", quote=F)
write.csv(cellranger.combined.sct@reductions$umap@cell.embeddings, file = "./cellranger.sct.umap.csv", quote=F)

write.csv(FetchData(qc.combined, "ident"), "./qc.strict.clusters.csv", quote=F)
write.csv(FetchData(qc.combined.sct, "ident"), "./qc.strict.sct.clusters.csv", quote=F)
write.csv(FetchData(cellranger.combined, "ident"), "./cellranger.clusters.csv", quote=F)
write.csv(FetchData(cellranger.combined.sct, "ident"), "./cellranger.sct.clusters.csv", quote=F)