## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, 
                       collapse = TRUE, comment = "#>",
                      fig.height = 5, fig.width = 5,
                      fig.retina = 1,
                      dpi = 60)

## ----eval=FALSE---------------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#  install.packages("BiocManager")
# }
# BiocManager::install("CiteFuse")

## -----------------------------------------------------------------------------
library(CiteFuse)
library(scater)
library(SingleCellExperiment)
library(DT)

## -----------------------------------------------------------------------------
data("CITEseq_example", package = "CiteFuse")
names(CITEseq_example)
lapply(CITEseq_example, dim)

## -----------------------------------------------------------------------------
sce_citeseq <- preprocessing(CITEseq_example)
sce_citeseq

## -----------------------------------------------------------------------------
sce_citeseq <- normaliseExprs(sce = sce_citeseq, 
                              altExp_name = "HTO", 
                              transform = "log")

## ----fig.height=6, fig.width=6------------------------------------------------

sce_citeseq <- scater::runTSNE(sce_citeseq, 
                               altexp = "HTO", 
                               name = "TSNE_HTO", 
                               pca = TRUE)


visualiseDim(sce_citeseq,
             dimNames = "TSNE_HTO") + labs(title = "tSNE (HTO)")

sce_citeseq <- scater::runUMAP(sce_citeseq, 
                               altexp = "HTO", 
                               name = "UMAP_HTO")


visualiseDim(sce_citeseq,
             dimNames = "UMAP_HTO") + labs(title = "UMAP (HTO)")

## -----------------------------------------------------------------------------
sce_citeseq <- crossSampleDoublets(sce_citeseq)

## -----------------------------------------------------------------------------
table(sce_citeseq$doubletClassify_between_label)
table(sce_citeseq$doubletClassify_between_class)

## ----fig.height=6, fig.width=6------------------------------------------------
visualiseDim(sce_citeseq, 
             dimNames = "TSNE_HTO", 
             colour_by = "doubletClassify_between_label")

## ----fig.height=8, fig.width=6------------------------------------------------
plotHTO(sce_citeseq, 1:4)

## -----------------------------------------------------------------------------
sce_citeseq <- withinSampleDoublets(sce_citeseq,
                                    minPts = 10)

## -----------------------------------------------------------------------------
table(sce_citeseq$doubletClassify_within_label)
table(sce_citeseq$doubletClassify_within_class)

## ----fig.height=6, fig.width=6------------------------------------------------
visualiseDim(sce_citeseq, 
             dimNames = "TSNE_HTO", 
             colour_by = "doubletClassify_within_label")

## -----------------------------------------------------------------------------
sce_citeseq <- sce_citeseq[, sce_citeseq$doubletClassify_within_class == "Singlet" & sce_citeseq$doubletClassify_between_class == "Singlet"]
sce_citeseq

## -----------------------------------------------------------------------------
sce_citeseq <- scater::logNormCounts(sce_citeseq)
sce_citeseq <- normaliseExprs(sce_citeseq, altExp_name = "ADT", transform = "log")
system.time(sce_citeseq <- CiteFuse(sce_citeseq))

## -----------------------------------------------------------------------------
SNF_W_clust <- spectralClustering(metadata(sce_citeseq)[["SNF_W"]], K = 20)
plot(SNF_W_clust$eigen_values)
which.max(abs(diff(SNF_W_clust$eigen_values)))

## -----------------------------------------------------------------------------
SNF_W_clust <- spectralClustering(metadata(sce_citeseq)[["SNF_W"]], K = 5)
sce_citeseq$SNF_W_clust <- as.factor(SNF_W_clust$labels)

SNF_W1_clust <- spectralClustering(metadata(sce_citeseq)[["ADT_W"]], K = 5)
sce_citeseq$ADT_clust <- as.factor(SNF_W1_clust$labels)

SNF_W2_clust <- spectralClustering(metadata(sce_citeseq)[["RNA_W"]], K = 5)
sce_citeseq$RNA_clust <- as.factor(SNF_W2_clust$labels)

## ----fig.height=8, fig.width=8------------------------------------------------

sce_citeseq <- reducedDimSNF(sce_citeseq,
                             method = "tSNE", 
                             dimNames = "tSNE_joint")

g1 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint", colour_by = "SNF_W_clust") +
  labs(title = "tSNE (SNF clustering)")
g2 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint",  colour_by = "ADT_clust") +
  labs(title = "tSNE (ADT clustering)")
g3 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint",  colour_by = "RNA_clust") +
  labs(title = "tSNE (RNA clustering)")

library(gridExtra)
grid.arrange(g3, g2, g1, ncol = 2)

## ----fig.height=8-------------------------------------------------------------
g1 <- visualiseDim(sce_citeseq, dimNames = "tSNE_joint", 
                   colour_by = "hg19_CD8A",
                   data_from = "assay",
                   assay_name = "logcounts") +
  labs(title = "tSNE: hg19_CD8A (RNA expression)")
g2 <- visualiseDim(sce_citeseq,dimNames = "tSNE_joint", 
                   colour_by = "CD8",
                   data_from = "altExp",
                   altExp_assay_name = "logcounts") +
  labs(title = "tSNE: CD8 (ADT expression)")

grid.arrange(g1, g2, ncol = 2)

## ----fig.height=8, fig.width=8------------------------------------------------
SNF_W_louvain <- igraphClustering(sce_citeseq, method = "louvain")
table(SNF_W_louvain)

sce_citeseq$SNF_W_louvain <- as.factor(SNF_W_louvain)

visualiseDim(sce_citeseq, dimNames = "tSNE_joint", colour_by = "SNF_W_louvain") +
  labs(title = "tSNE (SNF louvain clustering)")

## ----fig.height = 6, fig.width = 6--------------------------------------------
visualiseKNN(sce_citeseq, colour_by = "SNF_W_louvain")

## ----fig.height = 4, fig.width = 8--------------------------------------------
visualiseExprs(sce_citeseq, 
               plot = "boxplot", 
               group_by = "SNF_W_louvain",
               feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))
visualiseExprs(sce_citeseq, 
               plot = "violin", 
               group_by = "SNF_W_louvain",
               feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))
visualiseExprs(sce_citeseq, 
               plot = "jitter", 
               group_by = "SNF_W_louvain",
               feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))
visualiseExprs(sce_citeseq, 
               plot = "density", 
               group_by = "SNF_W_louvain",
               feature_subset = c("hg19_CD2", "hg19_CD4", "hg19_CD8A", "hg19_CD19"))



## ----fig.height = 4, fig.width = 6--------------------------------------------
visualiseExprs(sce_citeseq, 
               altExp_name = "ADT", 
               group_by = "SNF_W_louvain",
               plot = "violin", n = 5)
visualiseExprs(sce_citeseq, altExp_name = "ADT", 
               plot = "jitter", 
               group_by = "SNF_W_louvain", 
               feature_subset = c("CD2", "CD8", "CD4", "CD19"))
visualiseExprs(sce_citeseq, altExp_name = "ADT", 
               plot = "density", 
               group_by = "SNF_W_louvain",
               feature_subset = c("CD2", "CD8", "CD4", "CD19"))


## ----fig.height = 4, fig.width = 8--------------------------------------------
visualiseExprs(sce_citeseq, altExp_name = "ADT", 
               plot = "pairwise", 
               feature_subset = c("CD4", "CD8"))

visualiseExprs(sce_citeseq, altExp_name = "ADT", 
               plot = "pairwise", 
               feature_subset = c("CD45RA", "CD4", "CD8"), 
               threshold = rep(4, 3))

## -----------------------------------------------------------------------------
# DE will be performed for RNA if altExp_name = "none" 
sce_citeseq <- DEgenes(sce_citeseq,
                       altExp_name = "none", 
                       group = sce_citeseq$SNF_W_louvain,
                       return_all = TRUE,
                       exprs_pct = 0.5)


sce_citeseq <- selectDEgenes(sce_citeseq,
                             altExp_name = "none")
datatable(format(do.call(rbind, metadata(sce_citeseq)[["DE_res_RNA_filter"]]), 
                 digits = 2))

## -----------------------------------------------------------------------------
sce_citeseq <- DEgenes(sce_citeseq,
                       altExp_name = "ADT", 
                       group = sce_citeseq$SNF_W_louvain,
                       return_all = TRUE,
                       exprs_pct = 0.5)


sce_citeseq <- selectDEgenes(sce_citeseq,
                             altExp_name = "ADT")
datatable(format(do.call(rbind, metadata(sce_citeseq)[["DE_res_ADT_filter"]]), 
                 digits = 2))

## ----fig.height = 10, fig.width = 10------------------------------------------
rna_DEgenes <- metadata(sce_citeseq)[["DE_res_RNA_filter"]]
adt_DEgenes <- metadata(sce_citeseq)[["DE_res_ADT_filter"]]

rna_DEgenes <- lapply(rna_DEgenes, function(x){
  x$name <- gsub("hg19_", "", x$name)
  x})
DEbubblePlot(list(RNA = rna_DEgenes, ADT = adt_DEgenes))

## ----fig.height = 8, fig.width = 8--------------------------------------------
rna_list <- c("hg19_CD4",
              "hg19_CD8A",
              "hg19_HLA-DRB1",
              "hg19_ITGAX",
              "hg19_NCAM1",
              "hg19_CD27",
              "hg19_CD19")

adt_list <- c("CD4", "CD8", "MHCII (HLA-DR)", "CD11c", "CD56", "CD27", "CD19")

rna_DEgenes_all <- metadata(sce_citeseq)[["DE_res_RNA"]]
adt_DEgenes_all <- metadata(sce_citeseq)[["DE_res_ADT"]]

feature_list <- list(RNA = rna_list, ADT = adt_list)
de_list <- list(RNA = rna_DEgenes_all, ADT = adt_DEgenes_all)

DEcomparisonPlot(de_list = de_list,
                 feature_list = feature_list)

## ----fig.height = 8, fig.width = 8--------------------------------------------
set.seed(2020)
sce_citeseq <- importanceADT(sce_citeseq, 
                             group = sce_citeseq$SNF_W_louvain,
                             subsample = TRUE)

visImportance(sce_citeseq, plot = "boxplot")
visImportance(sce_citeseq, plot = "heatmap")

sort(metadata(sce_citeseq)[["importanceADT"]], decreasing = TRUE)[1:20]

## -----------------------------------------------------------------------------
subset_adt <- names(which(metadata(sce_citeseq)[["importanceADT"]] > 5))
subset_adt

system.time(sce_citeseq <- CiteFuse(sce_citeseq,
                                    ADT_subset = subset_adt,
                                    metadata_names = c("W_SNF_adtSubset1",
                                                       "W_ADT_adtSubset1",
                                                       "W_RNA")))

SNF_W_clust_adtSubset1 <- spectralClustering(metadata(sce_citeseq)[["W_SNF_adtSubset1"]], K = 5)
sce_citeseq$SNF_W_clust_adtSubset1 <- as.factor(SNF_W_clust_adtSubset1$labels)


library(mclust)
adjustedRandIndex(sce_citeseq$SNF_W_clust_adtSubset1, sce_citeseq$SNF_W_clust)

## ----fig.height = 8, fig.width = 8--------------------------------------------
RNA_feature_subset <- unique(as.character(unlist(lapply(rna_DEgenes_all, "[[", "name"))))
ADT_feature_subset <- unique(as.character(unlist(lapply(adt_DEgenes_all, "[[", "name"))))

geneADTnetwork(sce_citeseq,
               RNA_feature_subset = RNA_feature_subset,
               ADT_feature_subset = ADT_feature_subset,
               cor_method = "pearson",
               network_layout = igraph::layout_with_fr)

## -----------------------------------------------------------------------------
data("lr_pair_subset", package = "CiteFuse")
head(lr_pair_subset)

sce_citeseq <- normaliseExprs(sce = sce_citeseq, 
                              altExp_name = "ADT", 
                              transform = "zi_minMax")

sce_citeseq <- normaliseExprs(sce = sce_citeseq, 
                              altExp_name = "none", 
                              exprs_value = "logcounts",
                              transform = "minMax")

sce_citeseq <- ligandReceptorTest(sce = sce_citeseq,
                                  ligandReceptor_list = lr_pair_subset,
                                  cluster = sce_citeseq$SNF_W_louvain,
                                  RNA_exprs_value = "minMax",
                                  use_alt_exp = TRUE,
                                  altExp_name = "ADT",
                                  altExp_exprs_value = "zi_minMax",
                                  num_permute = 1000) 

## ----fig.height=10, fig.width=8-----------------------------------------------
visLigandReceptor(sce_citeseq, 
                  type = "pval_heatmap",
                  receptor_type = "ADT")

## ----fig.height=12, fig.width=6-----------------------------------------------
visLigandReceptor(sce_citeseq, 
                  type = "pval_dotplot",
                  receptor_type = "ADT")

## ----fig.height=8, fig.width=8------------------------------------------------
visLigandReceptor(sce_citeseq, 
                  type = "group_network",
                  receptor_type = "ADT")

## ----fig.height=8, fig.width=8------------------------------------------------
visLigandReceptor(sce_citeseq, 
                  type = "group_heatmap",
                  receptor_type = "ADT")

## ----fig.height=8, fig.width=8------------------------------------------------
visLigandReceptor(sce_citeseq,
                  type = "lr_network",
                  receptor_type = "ADT")

## -----------------------------------------------------------------------------
data("sce_ctcl_subset", package = "CiteFuse")

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

visualiseExprsList(sce_list = list(control = sce_citeseq,
                                   ctcl = sce_ctcl_subset),
                   plot = "boxplot",
                   altExp_name = "none",
                   exprs_value = "logcounts",
                   feature_subset = c("hg19_S100A10", "hg19_CD8A"),
                   group_by = c("SNF_W_louvain", "SNF_W_louvain"))

visualiseExprsList(sce_list = list(control = sce_citeseq,
                                   ctcl = sce_ctcl_subset),
                   plot = "boxplot",
                   altExp_name = "ADT", 
                   feature_subset = c("CD19", "CD8"),
                   group_by = c("SNF_W_louvain", "SNF_W_louvain"))

## -----------------------------------------------------------------------------
de_res <- DEgenesCross(sce_list = list(control = sce_citeseq,
                                       ctcl = sce_ctcl_subset),
                       colData_name = c("SNF_W_louvain", "SNF_W_louvain"),
                       group_to_test = c("2", "6"))
de_res_filter <- selectDEgenes(de_res = de_res)
de_res_filter

## ----eval = FALSE-------------------------------------------------------------
# tmpdir <- tempdir()
# download.file("http://cf.10xgenomics.com/samples/cell-exp/3.1.0/connect_5k_pbmc_NGSC3_ch1/connect_5k_pbmc_NGSC3_ch1_filtered_feature_bc_matrix.tar.gz", file.path(tmpdir, "/5k_pbmc_NGSC3_ch1_filtered_feature_bc_matrix.tar.gz"))
# untar(file.path(tmpdir, "5k_pbmc_NGSC3_ch1_filtered_feature_bc_matrix.tar.gz"),
#       exdir = tmpdir)
# sce_citeseq_10X <- readFrom10X(file.path(tmpdir, "filtered_feature_bc_matrix/"))
# sce_citeseq_10X

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

