library(Seurat)
library(devtools)
library(loomR)
library(SCopeLoomR)
library(AUCell)
library(SCENIC)
library(dplyr)
library(KernSmooth)
library(RColorBrewer)
library(plotly)
library(BiocParallel)
library(grid)
library(ComplexHeatmap)
library(data.table)
library(openxlsx)

wdir <- "./SCENIC_CART"

scenicLoomPath <- file.path(paste(wdir, "Full_CART.from_tsv.loom", sep = "/"))
motifEnrichmentFile <- file.path(paste(wdir, "Full_CART.adjacences.tsv", sep = "/"))

loom <- open_loom(scenicLoomPath)
# Read information from loom file:
regulons_incidMat <- get_regulons(loom, column.attr.name = 'Regulons')
regulons <- regulonsToGeneLists(regulons_incidMat)
regulonAUC <- get_regulons_AUC(loom, column.attr.name = 'RegulonsAUC')
regulonAucThresholds <- get_regulon_thresholds(loom)
embeddings <- get_embeddings(loom)
close_loom(loom)

# Get metadata info
cellClusters <- read.xlsx(xlsxFile = paste(wdir, "Full_CART_metadata.xlsx", sep = "/"), rowNames = T)

### Fig6 B - SCENIC CART Heatmap ###
selectedResolution <- "PROP.Condition"
# Split the cells by cluster:
cellsPerCluster <- split(rownames(cellClusters), cellClusters[,selectedResolution]) 
regulonAUC <- regulonAUC[onlyNonDuplicatedExtended(rownames(regulonAUC)),]
# Calculate average expression:
regulonActivity_byCellType <- sapply(cellsPerCluster, function(cells) rowMeans(getAUC(regulonAUC)[,cells]))
# Save regulon activity
write.xlsx(x = list("RegulonActivity" = regulonActivity_byCellType), rowNames = T,
           file = paste(wdir, "Full_CART_HM_regulonActivity.xlsx", sep = "/"))

# Scale expression:
regulonActivity_byCellType_Scaled <- t(scale(t(regulonActivity_byCellType), center = T, scale=T))

# plot:
gg_index <- which(rownames(regulonActivity_byCellType_Scaled) %in%
                      c("Jund(+)","Foxo1(+)","Junb(+)","Klf2(+)","Fos(+)","Irf7(+)","Stat1(+)","Stat5b(+)","Klf9(+)","Irf4(+)","Jun(+)",
                        "Stat5a(+)","Nfkb1(+)","Fosb(+)","Ikzf1(+)","Tfdp2(+)","Myc(+)","Irf5(+)","Bcl3(+)","Stat2(+)","Irf9(+)",
                        "Nelfe(+)","Zscan29(+)","Zfp454(+)","Zic1(+)","Thap11(+)","Batf3(+)","Myb(+)","Tbx21(+)","Foxp1(+)","Eomes(+)"))
gg_lab <- rownames(regulonActivity_byCellType_Scaled)[gg_index]
hm_aanot <- rowAnnotation(link = anno_mark(at = gg_index, labels = gg_lab, labels_gp = gpar(fontsize = 10)))

options(repr.plot.width = 8, repr.plot.height = 10)
png(file = paste(wdir, "Full_CART_HM.png", sep = "/"))
draw(ComplexHeatmap::Heatmap(matrix = regulonActivity_byCellType_Scaled,
                             name = "Regulon activity",
                             column_dend_side = "bottom",
                             column_names_side = "top",
                             column_names_rot = 0,
                             show_row_names = FALSE,
                             right_annotation = hm_aanot)) # row font size
dev.off()