## ---- echo=FALSE, results="hide", message=FALSE-------------------------------
require(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)
library(BiocStyle)

## -----------------------------------------------------------------------------
library(scater)
library(scran)
library(scRNAseq)
sce <- KotliarovPBMCData()
sce <- sce[,1:1000] # subset for speed.
sce

## -----------------------------------------------------------------------------
stats <- perCellQCMetrics(sce, subsets=list(Mito=grep("^MT-", rownames(sce))))
filter <- quickPerCellQC(stats, sub.fields="subsets_Mito_percent")
sce <- sce[,!filter$discard]
sce <- logNormCounts(sce)
dec <- modelGeneVar(sce)
set.seed(10000)
sce <- runPCA(sce, ncomponents=25, subset_row=getTopHVGs(dec, n=5000))

## -----------------------------------------------------------------------------
# TODO: sync this with the book.
sceA <- altExp(sce)
statsA <- perCellQCMetrics(sceA)
keep <- statsA$detected > 50 
sceA <- sceA[,keep]

library(DropletUtils)
amb <- inferAmbience(assay(sceA))
sceA <- computeMedianFactors(sceA, reference=amb)
sceA <- logNormCounts(sceA)

set.seed(10000)
sceA <- runPCA(sceA, ncomponents=15)

## -----------------------------------------------------------------------------
sce2 <- sce[,keep]
altExp(sce2) <- sceA

## -----------------------------------------------------------------------------
library(mumosa)
output <- rescaleByNeighbors(list(reducedDim(sce2), reducedDim(altExp(sce2))))
dim(output)

## -----------------------------------------------------------------------------
set.seed(100001)
library(bluster)
sce2$combined.clustering <- clusterRows(output, NNGraphParam())

reducedDim(sce2, "combined") <- output
sce2 <- runTSNE(sce2, dimred="combined")
plotTSNE(sce2, colour_by="combined.clustering")

## -----------------------------------------------------------------------------
g.rna <- buildSNNGraph(sce2, use.dimred="PCA")
rna.clusters <- igraph::cluster_walktrap(g.rna)$membership
table(rna.clusters)

g.adt <- buildSNNGraph(altExp(sce2), use.dimred='PCA')
adt.clusters <- igraph::cluster_walktrap(g.adt)$membership
table(adt.clusters)

## -----------------------------------------------------------------------------
intersected <- paste0(adt.clusters, ".", rna.clusters)
table(intersected)

## -----------------------------------------------------------------------------
intersected2 <- intersectClusters(list(rna.clusters, adt.clusters),
    list(reducedDim(sce2), reducedDim(altExp(sce2))))
table(intersected2)

## -----------------------------------------------------------------------------
set.seed(100002)
umap.out <- calculateMultiUMAP(list(reducedDim(sce2), reducedDim(altExp(sce2))), 
    n_components=20)
dim(umap.out)

## -----------------------------------------------------------------------------
library(bluster)
sce2$umap.clustering <- clusterRows(umap.out, NNGraphParam())

## -----------------------------------------------------------------------------
set.seed(100002)
sce2 <- runMultiUMAP(sce2, dimred="PCA", extras=list(reducedDim(altExp(sce2))))
plotReducedDim(sce2, "MultiUMAP", colour_by="umap.clustering")

## -----------------------------------------------------------------------------
g.com <- intersectGraphs(g.rna, g.adt)

## -----------------------------------------------------------------------------
com.clusters <- igraph::cluster_walktrap(g.com)$membership
table(com.clusters)

## -----------------------------------------------------------------------------
cor.all <- computeCorrelations(sce2, altExp(sce2)[1:5,])
cor.all[order(cor.all$p.value),]

## -----------------------------------------------------------------------------
b <- rep(1:3, length.out=ncol(sce2))
cor.all.batch <- computeCorrelations(sce2, altExp(sce2)[1:5,], block=b)

## -----------------------------------------------------------------------------
set.seed(100001) # for IRLBA.
top.cor <- findTopCorrelations(sce2[1:100,], y=altExp(sce2), number=10)
top.cor
top.cor$positive

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