# This R environment comes with many helpful analytics packages installed
# It is defined by the kaggle/rstats Docker image: https://github.com/kaggle/docker-rstats
# For example, here's a helpful package to load

library(tidyverse) # metapackage of all tidyverse packages

# Input data files are available in the read-only "../input/" directory
# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory

list.files(path = "../input")

# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using "Save & Run All" 
# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session

#Libraries
library(arrow)
library(data.table)
library(Matrix)
library(dplyr)
library(tidyr)
library(Seurat)
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)


#Get the multi-omics data
train <- read_parquet("multiome_train.parquet")
multi <- read.csv("multiome_var_meta.csv")
meta <- read.csv("multiome_obs_meta.csv")


#Separate the multi-omics data in 
#ATAC
nueva_train <- train[grep("^chr", train$location), ]
#RNA
no_chr <- train[!grepl("chr", train$location), ]



train <- nueva_train

#######SPARSE MATRIX PEAKS
customers <- unique(train$location)
products <- unique(train$obs_id)
train$row <- match(train$location
                   , customers)
train$col <- match(train$obs_id
                   , products)
df_sparse2 <- sparseMatrix(
  i = train$row, 
  j = train$col,
  x = train$normalized_count,
  dimnames = list(customers,
                  products)
)
dim(df_sparse2)

################SPARSE RNA MATRIX

train <- no_chr
customers <- unique(train$location)
products <- unique(train$obs_id)
train$row <- match(train$location
                   , customers)
train$col <- match(train$obs_id
                   , products)
df_sparse <- sparseMatrix(
  i = train$row, 
  j = train$col,
  x = train$normalized_count,
  dimnames = list(customers,
                  products)
)

################# ANALYSIS scATAC-seq following signac multiomic
#https://stuartlab.org/signac/articles/pbmc_multiomic


# I don't know if this are the fragments, or the peaks
#Those are extracted from the "train file"
fragpath <- "~/Github/OpenProblems/final_ordenat.tsv.gz"



####### RNA OBJECT
pbmc <- CreateSeuratObject(
  counts = df_sparse,
  assay = "RNA",
  meta.data=meta
)

annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevels(annotation) <- paste0('chr', seqlevels(annotation))



###### PEAKS OBJECT

# create ATAC assay and add it to the object
pbmc[["ATAC"]] <- CreateChromatinAssay(
  counts = df_sparse2,
  sep = c(":", "-"),
  fragments = fragpath,
  annotation = annotation
)


######### Quality control 
DefaultAssay(pbmc) <- "ATAC"

pbmc <- NucleosomeSignal(pbmc)
#See the nucleosomeSignal output, only 1 and Inf.
pbmc <- TSSEnrichment(pbmc)


VlnPlot(
  object = pbmc,
  features = c("nCount_RNA", "nCount_ATAC", "TSS.enrichment"),
  ncol = 4,
  pt.size = 0
)


########################
# filter out low quality cells, changeable numbers, 
#pbmc <- subset(
#  x = pbmc,
#  subset = nCount_ATAC < 100000 &
 #   nCount_RNA < 25000 &
#    nCount_ATAC > 1000 &
 #   nCount_RNA > 1000 &
  #  nucleosome_signal < 2 &   ????? Should we know anything about that output?
  #  TSS.enrichment > 1
#)
#pbmc





################################ PEAKS CALLING????
#Do we have to do a peak calling?



############TREATING THE "COUNTS" AS PEAKS!
DefaultAssay(pbmc) <- "RNA"
pbmc <- SCTransform(pbmc)
pbmc <- RunPCA(pbmc)

DefaultAssay(pbmc) <- "ATAC"
pbmc <- FindTopFeatures(pbmc, min.cutoff = 5)
pbmc <- RunTFIDF(pbmc)
pbmc <- RunSVD(pbmc)





DimPlot(pbmc, group.by = "cell_type")