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')


eryidx <- read.table(file='../notebooks/output/20231128_bthalcombo_v6c_eryplus_traj_cell_idx.txt', header = FALSE)
eryidx <- as.vector(eryidx$V1)


pseudotime <- ann[eryidx]$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[eryidx]$obs$donor, rownames(conv)) ]


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

sce <- sce[,eryidx]


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=5,
                 verbose=T
                 #parallel=TRUE, 
                 #BPPARAM = BPPARAM
                )

Sys.time()

saveRDS(file="../output/20231129_sceGAM_erytroid_plus_traj_conditions.rds", sceGAM)