# Run seurat to combine h5ad datasets

# load libraries

library(scater)
library(Seurat)
library(cowplot)
library(reticulate)
library(SeuratData)
library(SeuratDisk)
library(anndata)
library(optparse)
 
option_list = list(
  make_option(c("--input1"), type="character", default=NULL,
              help="dataset1 file name", metavar=NULL),
  make_option(c("--input2"), type="character", default=NULL,
              help="dataset2 file name", metavar=NULL),
  make_option(c("--qc1"), type="character", default=NULL,
              help="qc input1 from anndata.h5ad", metavar=NULL),
  make_option(c("--qc2"), type="character", default=NULL,
              help="qc input2 from anndata.h5ad", metavar=NULL)
);

opt_parser = OptionParser(option_list=option_list);
opt = parse_args(opt_parser);

if (is.null(opt$input1) || is.null(opt$input2) || is.null(opt$qc1) || is.null(opt$qc2)){
#if (is.null(opt$input1) || is.null(opt$qc1)){
  print_help(opt_parser)
  stop("At least one argument must be supplied (input file).n", call.=FALSE)
}

# Read cellranger input
cellranger.input1 <- Read10X(data.dir = opt$input1)
cellranger.input2 <- Read10X(data.dir = opt$input2)

# Read QC data
qc.input1 <- read_h5ad(opt$qc1)
qc.input2 <- read_h5ad(opt$qc2)

# load objects into seurat
cellranger.seurat1 <- CreateSeuratObject(counts = cellranger.input1, project = "Nodule1_cellranger")
cellranger.seurat2 <- CreateSeuratObject(counts = cellranger.input2, project = "Nodule2_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")

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

# Integrate objects
cellranger.combined <- merge(cellranger.seurat1, y = cellranger.seurat2, add.cell.ids = c("Nodule1_cellranger", "Nodule2_cellranger"), project = "Nodule_cellranger")
qc.combined <- merge(qc.seurat1, y = qc.seurat2, add.cell.ids = c("Nodule1_qc", "Nodule2_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)
})

## 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)

# 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)

str(cellranger.features)
str(qc.features)

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)

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