# 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="character"),
  make_option(c("--input2"), type="character", default=NULL,
              help="dataset2 file name", metavar="character")
  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="character")
);

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)){
  print_help(opt_parser)
  stop("At least one argument must be supplied (input file).n", call.=FALSE)
}

# Read cellranger input
cellranger.input1 <- opt$input1
cellranger.input2 <- 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 = t(cellranger.input1$X), meta.data = cellranger.input1$obs, project = "wholerootinoc")
cellranger.seurat2 <- CreateSeuratObject(counts = t(cellranger.input2$X), meta.data = cellranger.input2$obs, project = "wholerootinoc")

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

seurat.mock.wholeroot.qc <- CreateSeuratObject(counts = t(mock.wholeroot.qc$X), meta.data = mock.wholeroot.qc$obs, project = "wholerootmock")

# Wholeroot QC combined
wholeroot.combined.qc <- merge(seurat.input1, y = seurat.mock.wholeroot.qc, add.cell.ids = c("wholerootinoc", "wholerootmock"), project = "wholeroot")
str(wholeroot.combined.qc)

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