prefix <- format(as.Date(Sys.time()), "%Y%m%d")
prefix

basedir <- '.'


library(tradeSeq)
library(RColorBrewer)
library(SingleCellExperiment)


ann <- anndata::read_h5ad('../h5ad/20231122_bthalcombo_v6c_cr2_no_B_palantir_pseudotime_kernel.h5ad')


monoidx <- read.table(file='../notebooks/output/20231207_bthalcombo_v6c_monoplus_traj_cell_idx_CLEAN.txt', header = FALSE)
monoidx <- as.vector(monoidx$V1)


pseudotime <- ann[monoidx]$obs$palantir_pseudotime


conv <- t( data.frame("bthal005" = "bthal", "bthal006" = "bthal", "healthy_CTRL" = "healthy", "bthal010" = "bthal",
          "bthal009" = "bthal", "bthal007" = "bthal", "bthal008" = "bthal", "P185" = "healthy",
          "P181" = "healthy", "P257" = "healthy", "paedBM1" = "healthy", "paedBM2" = "healthy") )

condition <- conv[ match(ann[monoidx]$obs$donor, rownames(conv)) ]


sce <- readRDS( paste0(basedir, 'h5ad/', '20230815_bthalcombo_v6b_counts_sf_SCE_obj.rds') )

sce <- sce[,monoidx]


cw <- rep(1,ncol(sce)) # cell weights

rm(ann)

set.seed(42)


Sys.time()

sceGAM <- fitGAM(counts = as.matrix( counts(sce) ),
                 conditions = factor(condition),
                 pseudotime=pseudotime,
                 cellWeights=cw,
                 nknots=6,
                 verbose=T
                 #parallel=TRUE, 
                 #BPPARAM = BPPARAM
                )

Sys.time()

saveRDS(file="../output/20231207_sceGAM_mono_plus_traj_conditions.rds", sceGAM)