## ----maits_load, message=FALSE, warning=FALSE------------------------------
suppressPackageStartupMessages({
  library(MAST)
  library(singleCellTK)
  library(xtable)
})

data(maits, package="MAST")
maits_sce <- createSCE(assayFile = t(maits$expressionmat),
                       annotFile = maits$cdat,
                       featureFile = maits$fdat,
                       assayName = "logtpm",
                       inputDataFrames = TRUE,
                       createLogCounts = FALSE)
rm(maits)

## ----maits_summarize, results='asis'---------------------------------------
knitr::kable(summarizeTable(maits_sce, useAssay = "logtpm"))

## ----maits_colnames--------------------------------------------------------
colnames(colData(maits_sce))
table(colData(maits_sce)$ourfilter)

## ----maits_filter----------------------------------------------------------
maits_subset <- maits_sce[, colData(maits_sce)$ourfilter]
table(colData(maits_subset)$ourfilter)

## ----maits_filter_table, results='asis'------------------------------------
knitr::kable(summarizeTable(maits_subset, useAssay = "logtpm"))

## ----maits_availablereduceddims--------------------------------------------
reducedDims(maits_subset)

## ----maits_getpcatsne------------------------------------------------------
maits_subset <- getPCA(maits_subset, useAssay = "logtpm",
                       reducedDimName = "PCA_logtpm")
maits_subset <- getTSNE(maits_subset, useAssay = "logtpm",
                        reducedDimName = "TSNE_logtpm")
reducedDims(maits_subset)

## ----maits_pca-------------------------------------------------------------
plotPCA(maits_subset, reducedDimName = "PCA_logtpm", colorBy = "condition")

## ----maits_tsne------------------------------------------------------------
plotTSNE(maits_subset, reducedDimName = "TSNE_logtpm", colorBy = "condition")

## ----maits_convert_symbols, message=FALSE----------------------------------
suppressPackageStartupMessages({
  library(org.Hs.eg.db)
})
maits_entrez <- maits_subset
maits_subset <- convertGeneIDs(maits_subset, inSymbol = "ENTREZID",
                               outSymbol = "SYMBOL", database = "org.Hs.eg.db")
#to remove confusion for MAST about the gene name:
rowData(maits_subset)$primerid <- NULL

## ----maits_thresh, fig.height=8, message=FALSE-----------------------------
thresholds <- thresholdGenes(maits_subset, useAssay = "logtpm")
par(mfrow = c(5, 4))
plot(thresholds)
par(mfrow = c(1, 1))

## ----maits_MAST, message=FALSE---------------------------------------------
mast_results <- MAST(maits_subset, condition = "condition", useThresh = TRUE,
                     useAssay = "logtpm")

## ----maits_violin, fig.height=8, message=FALSE-----------------------------
MASTviolin(maits_subset, useAssay = "logtpm", fcHurdleSig = mast_results,
           threshP = TRUE, condition = "condition")

## ----maits_lm, fig.height=8, message=FALSE---------------------------------
MASTregression(maits_subset, useAssay = "logtpm", fcHurdleSig = mast_results,
               threshP = TRUE, condition = "condition")

## ----maits_heatmap---------------------------------------------------------
plotDiffEx(maits_subset, useAssay = "logtpm", condition = "condition",
           geneList = mast_results$Gene[1:100], annotationColors = "auto",
           displayRowLabels = FALSE, displayColumnLabels = FALSE)

## ----maits_gsva, message=FALSE---------------------------------------------
gsvaRes <- gsvaSCE(maits_entrez, useAssay = "logtpm",
                   "MSigDB c2 (Human, Entrez ID only)",
                   c("KEGG_PROTEASOME",
                     "REACTOME_VIF_MEDIATED_DEGRADATION_OF_APOBEC3G",
                     "REACTOME_P53_INDEPENDENT_DNA_DAMAGE_RESPONSE",
                     "BIOCARTA_PROTEASOME_PATHWAY",
                     "REACTOME_METABOLISM_OF_AMINO_ACIDS",
                     "REACTOME_REGULATION_OF_ORNITHINE_DECARBOXYLASE",
                     "REACTOME_CYTOSOLIC_TRNA_AMINOACYLATION",
                     "REACTOME_STABILIZATION_OF_P53",
                     "REACTOME_SCF_BETA_TRCP_MEDIATED_DEGRADATION_OF_EMI1"),
                    parallel.sz=1)

set.seed(1234)
gsvaPlot(maits_subset, gsvaRes, "Violin", "condition")
gsvaPlot(maits_subset, gsvaRes, "Heatmap", "condition")

## ----load_bladderbatch, message=FALSE--------------------------------------
library(bladderbatch)
data(bladderdata)
dat <- bladderEset

pheno <- pData(dat)
edata <- exprs(dat)
bladder_sctke <- createSCE(assayFile = edata,
                           annotFile = pheno,
                           assayName = "microarray",
                           inputDataFrames = TRUE,
                           createLogCounts = FALSE)

## ----plot_var_microarray, message=FALSE------------------------------------
plotBatchVariance(bladder_sctke, useAssay="microarray",
                  batch="batch", condition = "cancer")

## ----run_combat, message=FALSE---------------------------------------------
assay(bladder_sctke, "combat") <- ComBatSCE(inSCE = bladder_sctke,
                                            batch = "batch",
                                            useAssay = "microarray",
                                            covariates = "cancer")

## ----plot_var_postcombat, message=FALSE------------------------------------
plotBatchVariance(bladder_sctke, useAssay="combat",
                  batch="batch", condition = "cancer")

## ----sessionInfo, echo=FALSE-----------------------------------------------
sessionInfo()