# Seurat based analysis of GM12878 paired data scRNA component
# Preparation for integration with ArchR
# References
# https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html
# https://satijalab.org/seurat/v3.1/merge_vignette.html
# https://learn.gencore.bio.nyu.edu/single-cell-rnaseq/loading-your-own-data-in-seurat-reanalyze-a-different-dataset/
rm(list = ls(all.names = TRUE))  # Clear environment

library(dplyr)
library(Seurat)

option_list = list(
  make_option(c("-i", "--input"), type="character", default=NULL,
              help="Input csv file", metavar="character"),
  make_option(c("-p", "--project"), type="character", default="proj",
              help="Project name", metavar="character"),
  make_option(c("-o", "--output"), type="character", default=NULL, 
              help="output file name", metavar = "character")
);
opt_parser = OptionParser(option_list=option_list);
opt = parse_args(opt_parser);

raw_counts <- read.table(file=opt$input, sep=",")
sobj <- CreateSeuratObject(raw.data = raw_counts, min.cells = 3, min.genes = 200, project = opt$project)

sobj[["percent.mt"]] <- PercentageFeatureSet(sobj, pattern = "^MT-")
VlnPlot(sobj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# Remove some unwanted cells
sobj <- subset(sobj, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 25)

# Normalize
sobj <- NormalizeData(sobj, normalization.method = "LogNormalize", scale.factor = 10000)

# Find variable features
sobj <- FindVariableFeatures(sobj, selection.method = "vst", nfeatures = 2000)

# Scale
sobj <- ScaleData(sobj, features = rownames(sobj))

# PCA
sobj <- RunPCA(sobj, features = VariableFeatures(object = sobj))
DimPlot(sobj, reduction="pca")
ElbowPlot(sobj)

# Cluster the cells
sobj <- FindNeighbors(sobj, dims=1:10)  # Appears to be the same for 20 or 10 dimensions
sobj <- FindClusters(sobj, resolution=1)

# Plot
sobj <- RunUMAP(sobj, dims=1:20)
DimPlot(sobj, reduction="umap")

# Save
saveRDS(sobj, opt$output)