## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE
)

## ----results='hide', message=FALSE, eval=FALSE--------------------------------
#  if (!require("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  
#  BiocManager::install("epiregulon")
#  

## ----setup, message=FALSE, eval=FALSE-----------------------------------------
#  devtools::install_github(repo ='xiaosaiyao/epiregulon')

## -----------------------------------------------------------------------------
library(epiregulon)

## -----------------------------------------------------------------------------
# load the MAE object
library(scMultiome)
library(beachmat.hdf5)

mae <- scMultiome::reprogramSeq()

# extract peak matrix
PeakMatrix <- mae[["PeakMatrix"]]

# extract expression matrix
GeneExpressionMatrix <- mae[["GeneExpressionMatrix"]]
rownames(GeneExpressionMatrix) <- rowData(GeneExpressionMatrix)$name

# define the order of hash_assigment
GeneExpressionMatrix$hash_assignment <- 
  factor(as.character(GeneExpressionMatrix$hash_assignment),
         levels = c("HTO10_GATA6_UTR", "HTO2_GATA6_v2", "HTO8_NKX2.1_UTR", "HTO3_NKX2.1_v2",
                    "HTO1_FOXA2_v2", "HTO4_mFOXA1_v2", "HTO6_hFOXA1_UTR", "HTO5_NeonG_v2"))
# extract dimensional reduction matrix
reducedDimMatrix <- reducedDim(mae[['TileMatrix500']], "LSI_ATAC")

# transfer UMAP_combined from TileMatrix to GeneExpressionMatrix
reducedDim(GeneExpressionMatrix, "UMAP_Combined") <- 
  reducedDim(mae[['TileMatrix500']], "UMAP_Combined")

## -----------------------------------------------------------------------------

scater::plotReducedDim(GeneExpressionMatrix, 
                       dimred = "UMAP_Combined", 
                       text_by = "Clusters", 
                       colour_by = "Clusters")


## ----getTFMotifInfo-----------------------------------------------------------
grl <- getTFMotifInfo(genome = "hg38")
grl

## ----eval=FALSE---------------------------------------------------------------
#  grl.motif <- getTFMotifInfo(genome = "hg38",
#                              mode = "motif",
#                              peaks = rowRanges(PeakMatrix))

## ----calculateP2G-------------------------------------------------------------
set.seed(1010)
p2g <- calculateP2G(peakMatrix = PeakMatrix, 
                    expMatrix = GeneExpressionMatrix, 
                    reducedDim = reducedDimMatrix,
                    exp_assay = "normalizedCounts",
                    peak_assay = "counts")

p2g


## ----addTFMotifInfo-----------------------------------------------------------

overlap <- addTFMotifInfo(grl = grl, p2g = p2g, 
                          peakMatrix = PeakMatrix)
head(overlap)

## ----warning=FALSE, getRegulon------------------------------------------------

regulon <- getRegulon(p2g = p2g, overlap = overlap, 
                      aggregate = FALSE)
regulon


## ----network pruning, results = "hide", message = FALSE-----------------------

pruned.regulon <- pruneRegulon(expMatrix = GeneExpressionMatrix,
                               exp_assay = "normalizedCounts",
                               peakMatrix = PeakMatrix,
                               peak_assay = "counts",
                               test = "chi.sq",
                               regulon = regulon[regulon$tf %in% 
                                         c("NKX2-1","GATA6","FOXA1","FOXA2", "AR"),],
                               clusters = GeneExpressionMatrix$Clusters,
                               prune_value = "pval",
                               regulon_cutoff = 0.05
                               )

pruned.regulon

## ----addWeights, results = "hide", warning = FALSE, message = FALSE-----------

regulon.w <- addWeights(regulon = pruned.regulon,
                        expMatrix  = GeneExpressionMatrix,
                        exp_assay  = "normalizedCounts",
                        peakMatrix = PeakMatrix,
                        peak_assay = "counts",
                        clusters = GeneExpressionMatrix$Clusters,
                        method = "wilcox")

regulon.w


## ----addMotifScore------------------------------------------------------------

regulon.w.motif <- addMotifScore(regulon = regulon.w,
                                 peaks = rowRanges(PeakMatrix),
                                 species = "human",
                                 genome = "hg38")

# if desired, set weight to 0 if no motif is found
regulon.w.motif$weight[regulon.w.motif$motif == 0] <- 0

regulon.w.motif

## ----add logFC----------------------------------------------------------------
regulon.w <- addLogFC(regulon = regulon.w,
                      clusters = GeneExpressionMatrix$hash_assignment,
                      expMatrix  = GeneExpressionMatrix,
                      assay.type  = "normalizedCounts",
                      pval.type = "any",
                      logFC_condition = c("HTO10_GATA6_UTR", 
                                          "HTO2_GATA6_v2", 
                                          "HTO8_NKX2.1_UTR"),
                      logFC_ref = "HTO5_NeonG_v2")

regulon.w

## ----calculateActivity--------------------------------------------------------
score.combine <- calculateActivity(expMatrix = GeneExpressionMatrix,
                                   regulon = regulon.w, 
                                   mode = "weight", 
                                   method = "weightedMean", 
                                   exp_assay = "normalizedCounts",
                                   normalize = FALSE)

score.combine[1:5,1:5]

## -----------------------------------------------------------------------------
sessionInfo()