##### download sample files

url = "https://genome.med.nyu.edu/results/external/iCellR/example7_Spatial_Transcriptomic/V1_Mouse_Brain_Sagittal_Anterior_Section_2_filtered_feature_bc_matrix.tar.gz"

# download the file
download.file(url = url,
     destfile = "V1_Mouse_Brain_Sagittal_Anterior_Section_2_filtered_feature_bc_matrix.tar.gz",
     method = "auto")


url ="https://genome.med.nyu.edu/results/external/iCellR/example7_Spatial_Transcriptomic/V1_Mouse_Brain_Sagittal_Anterior_Section_2_spatial.tar.gz"

# download the file
download.file(url = url,
     destfile = "V1_Mouse_Brain_Sagittal_Anterior_Section_2_spatial.tar.gz",
     method = "auto")

url ="https://genome.med.nyu.edu/results/external/iCellR/example7_Spatial_Transcriptomic/V1_Mouse_Brain_Sagittal_Posterior_Section_2_filtered_feature_bc_matrix.tar.gz"

# download the file
download.file(url = url,
     destfile = "V1_Mouse_Brain_Sagittal_Posterior_Section_2_filtered_feature_bc_matrix.tar.gz",
     method = "auto")

url ="https://genome.med.nyu.edu/results/external/iCellR/example7_Spatial_Transcriptomic/V1_Mouse_Brain_Sagittal_Posterior_Section_2_spatial.tar.gz"

# download the file
download.file(url = url,
     destfile = "V1_Mouse_Brain_Sagittal_Posterior_Section_2_spatial.tar.gz",
     method = "auto")

##### untar
untar("V1_Mouse_Brain_Sagittal_Anterior_Section_2_filtered_feature_bc_matrix.tar.gz")
untar("V1_Mouse_Brain_Sagittal_Anterior_Section_2_spatial.tar.gz")

file.rename("spatial","spatial_Anterior2")
file.rename("filtered_feature_bc_matrix","filtered_feature_bc_matrix_Anterior2")

untar("V1_Mouse_Brain_Sagittal_Posterior_Section_2_filtered_feature_bc_matrix.tar.gz")
untar("V1_Mouse_Brain_Sagittal_Posterior_Section_2_spatial.tar.gz")

file.rename("spatial","spatial_Posterior2")
file.rename("filtered_feature_bc_matrix","filtered_feature_bc_matrix_Posterior2")

######

library(iCellR)

Anterior2 <- load10x("filtered_feature_bc_matrix_Anterior2",gene.name = 2)
Posterior2 <- load10x("filtered_feature_bc_matrix_Posterior2",gene.name = 2)

# if you want to analyze both samples
Samples <- c("Anterior2","Posterior2")
my.data <- data.aggregation(samples = Samples, condition.names = Samples)

# if you want to analyze 1 sample
# my.data <- load10x("filtered_feature_bc_matrix_Posterior2",gene.name = 2)

my.obj <- make.obj(my.data)


Anterior2 <- capture.image.10x("spatial_Anterior2")
Posterior2 <- capture.image.10x("spatial_Posterior2")

# if you want to analyze both samples
Samples <- c("Anterior2","Posterior2")
my.obj <- add.10x.image(my.obj,
          image.data.list = Samples, condition.names = Samples)

# if one sample
# My.image <- image.capture.10x("Post2_spatial")
# my.obj <- add.10x.image(my.obj, image.data.list = "My.image")

my.obj

############################# rest is like regular scRNA-Seq
my.obj <- qc.stats(my.obj,
which.data = "raw.data")

head(my.obj@stats)
summary(my.obj@stats)

png('plot1_QC_stats.png',width = 15, height = 6, units = 'in', res = 300)
stats.plot(my.obj,
    out.name = "UMI-plot",
    interactive = FALSE,
    cell.color = "slategray3",
    cell.size = 1,
    cell.transparency = 0.5,
    box.color = "red",
    box.line.col = "green")
dev.off()

png('plot2_scatter.mito.umi.png',width = 6, height = 6, units = 'in', res = 300)
stats.plot(my.obj, plot.type = "point.mito.umi",
    interactive = F, out.name = "plot2_scatter.mito.umi")
dev.off()

png('plot3_scatter.gene.umi.png',width = 6, height = 6, units = 'in', res = 300)
stats.plot(my.obj, plot.type = "point.gene.umi",
    interactive = F, out.name = "plot3_scatter.gene.umi")
dev.off()

####################################### Filter
#If you don't want to filter
# my.obj@main.data <- my.obj@raw.data

my.obj <- cell.filter(my.obj,
min.mito = 0,
max.mito = 0.50,
min.genes = 200,
max.genes = 10000,
min.umis = 0,
max.umis = Inf)

message(my.obj@my.filters)

################ Normalize
my.obj <- norm.data(my.obj, norm.method = "ranked.glsf", top.rank = 500)

################################################################# make model
my.obj <- gene.stats(my.obj, which.data = "main.data")
#head(my.obj@gene.data[order(my.obj@gene.data$numberOfCells, decreasing = T),])
my.gene.data <- my.obj@gene.data
head(my.gene.data)
write.table((my.gene.data), file="gene.data.tsv", sep="\t", row.names =F)

png('plot4_gene.model.png',width = 6, height = 6, units = 'in', res = 300)
make.gene.model(my.obj, my.out.put = "plot",
    dispersion.limit = 1.5,
    base.mean.rank = 500,
    no.mito.model = T,
    mark.mito = T,
    interactive = F,
    no.cell.cycle = T,
    out.name = "gene.model")
dev.off()

my.obj <- make.gene.model(my.obj, my.out.put = "data",
    dispersion.limit = 1.5,
    base.mean.rank = 500,
    no.mito.model = T,
    mark.mito = T,
    interactive = F,
    no.cell.cycle = T,
    out.name = "gene.model")

head(my.obj@gene.model)
length(my.obj@gene.model)

####### dim redux
my.obj <- run.pca(my.obj, method = "gene.model", gene.list = my.obj@gene.model,data.type = "main")

# if batch alignment is required replace this with PCA
# my.obj <- iba(my.obj,dims = 1:30, k = 10,ba.method = "CPCA", method = "gene.model", gene.list = my.obj@gene.model)

my.obj <- run.pc.tsne(my.obj, dims = 1:10)

my.obj <- run.umap(my.obj, dims = 1:10)

my.obj <- run.knetl(my.obj, dims = 1:20, zoom =200, dim.redux = "umap")

###### cluster
# if you want to cluster based on KNeTL
# my.obj <- run.phenograph(my.obj,k = 100, data.type = "knetl")

# if you want to cluster based on PCA
my.obj <- iclust(my.obj, sensitivity = 100 ,dims = 1:10, data.type = "pca")

##### Order cluster numbers
# how.to.order = "random"
how.to.order = "distance"
my.obj <- clust.ord(my.obj,top.rank = 500, how.to.order = how.to.order)

####
G0 <- readLines(system.file('extdata', 'G0.txt', package = 'iCellR'))
G1S <- readLines(system.file('extdata', 'G1S.txt', package = 'iCellR'))
G2M <- readLines(system.file('extdata', 'G2M.txt', package = 'iCellR'))
M <- readLines(system.file('extdata', 'M.txt', package = 'iCellR'))
MG1 <- readLines(system.file('extdata', 'MG1.txt', package = 'iCellR'))
S <- readLines(system.file('extdata', 'S.txt', package = 'iCellR'))

my.obj <- cell.cycle(my.obj,scoring.List = c("G0","G1S","G2M","M","MG1","S"), scoring.method = "tirosh")

# save your object
save(my.obj, file = "my.obj.Robj")

##### dimentionality plots

A=spatial.plot(my.obj,col.by = "clusters",conds.to.plot = "Anterior2",interactive= F)
B=spatial.plot(my.obj,col.by = "clusters",conds.to.plot = "Posterior2",interactive= F)
C= cluster.plot(my.obj,plot.type = "tsne",interactive = F,cell.size = 0.5,cell.transparency = 1, anno.clust=T)
D= cluster.plot(my.obj,plot.type = "tsne",col.by = "conditions",interactive = F,cell.size = 0.5,cell.transparency = 1, anno.clust=T)
E=spatial.plot(my.obj,col.by = "gene", gene = c("Cd4"), conds.to.plot = "Anterior2",interactive= F, scaleValue = TRUE)
F=spatial.plot(my.obj,col.by = "gene", gene = c("Cd4"), conds.to.plot = "Posterior2",interactive= F, scaleValue = TRUE)

library(gridExtra)
png('AllClusts.png', width = 8, height = 8, units = 'in', res = 300)
grid.arrange(A,B,C,D,E,F)
dev.off()