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

## ----pkg_install, eval=FALSE--------------------------------------------------
# # first check to see if BiocManager is available
# if(!requireNamespace('BiocManager', quietly=TRUE)){
#   install.packages('BiocManager')
# }
# 
# BiocManager::install('carnation')

## ----presetup, eval=FALSE-----------------------------------------------------
# # install optional packages if not present
# pkgs_to_check <- c('airway', 'org.Hs.eg.db', 'DEGreport')
# for(pkg in pkgs_to_check){
#   if(!requireNamespace(pkg, quietly=TRUE)){
#     BiocManager::install(pkg)
#   }
# }

## ----lib, message=FALSE-------------------------------------------------------
library(airway)
library(DESeq2)
library(dplyr)
library(GeneTonic)
library(org.Hs.eg.db)

## ----airway_load--------------------------------------------------------------
data('airway')

## ----extract------------------------------------------------------------------
mat <- assay(airway)
cdata <- colData(airway)

## ----view---------------------------------------------------------------------
dim(mat)

## ----view_cnts----------------------------------------------------------------
head(mat)

## ----view_cdata---------------------------------------------------------------
cdata

## ----annot, message=FALSE-----------------------------------------------------
keytypes <- list('SYMBOL'='SYMBOL', 'ENTREZID'='ENTREZID')

anno_df <- do.call('cbind',
             lapply(keytypes, function(x){
               mapIds(org.Hs.eg.db,
                 column=x,
                 keys=rownames(mat),
                 keytype='ENSEMBL')
               })
             )

# convert to data frame
anno_df <- as.data.frame(anno_df)

## ----view_annot---------------------------------------------------------------
head(anno_df)

## ----dds----------------------------------------------------------------------
dds <- DESeqDataSetFromMatrix(mat,
                              colData=cdata,
                              design=~cell + dex)

## ----view_dds-----------------------------------------------------------------
dds

## ----save_dds-----------------------------------------------------------------
dds_list <- list(main=dds)

## ----vst----------------------------------------------------------------------
rld_list <- lapply(dds_list, function(x) varianceStabilizingTransformation(x, blind=TRUE))

## ----deseq--------------------------------------------------------------------
dds <- DESeq(dds)

## ----results------------------------------------------------------------------
resultsNames(dds)

## ----contrast1, message=FALSE-------------------------------------------------
cell_comparison <- lfcShrink(dds,
                             coef='cell_N61311_vs_N052611',
                             type='normal')

## ----contrast2, message=FALSE-------------------------------------------------
dex_comparison <- lfcShrink(dds,
                            contrast=c('dex', 'trt', 'untrt'),
                            type='normal')

## ----view_contrast------------------------------------------------------------
head(dex_comparison)

## ----res----------------------------------------------------------------------
res_list <- list(
        dex_trt_vs_untrt=list(
            res=dex_comparison,
            dds='main',
            label='dex, treated vs untreated'),
        cell_N61311_vs_N052611=list(
            res=cell_comparison,
            dds='main',
            label='cell, N61311 vs N052611')
        )

## ----bld_res------------------------------------------------------------------
res_list <- lapply(res_list, function(x){
              # save the rownames as a new 'gene' column
              x$res[[ 'gene' ]] <- rownames(x$res)

              # add 'SYMBOL' and 'ENTREZID' columns
              x$res[[ 'SYMBOL' ]] <- anno_df[rownames(x$res), 'SYMBOL']
              x$res[[ 'ENTREZID' ]] <- anno_df[rownames(x$res), 'ENTREZID']

              x
            })

## ----alpha--------------------------------------------------------------------
# padj cutoff
alpha <- 0.01

# log2FoldChange threshold; 1 == 2x difference
lfc_threshold <- 1

# list to save DE genes
de.genes <- lapply(res_list, function(x){
              # changed genes
              idx <- x$res$padj < alpha &
                     !is.na(x$res$padj) &
                     abs(x$res$log2FoldChange) >= lfc_threshold

              # return DE genes as a dataframe
              x$res[idx, c('gene', 'ENTREZID')]
            })

## ----enrich-------------------------------------------------------------------
# fe results from dex comparison
data(eres_dex, package='carnation')

# fe results from dex comparison
data(eres_cell, package='carnation')

# compile into a list
go_list <- list(dex_trt_vs_untrt=eres_dex, cell_N61311_vs_N052611=eres_cell)

# list to save functional enrichment results
enrich_list <- list()

# list to save a converted object for GeneTonic plots
genetonic <- list()

for(comp in names(res_list)){
    # NOTE: this is the command used to generate the functional
    # enrichment results
    #go.res <- clusterProfiler::enrichGO(
    #        gene=de.genes[[comp]][['ENTREZID']],
    #        keyType='ENTREZID',
    #        OrgDb=org.Hs.eg.db,
    #        ont='BP',
    #        pvalueCutoff=1, qvalueCutoff=1,
    #        readable = TRUE)

    # we use the precomputed results here instead
    go.res <- go_list[[ comp ]]

    enrich_list[[ comp ]] <- list(
                               res=comp,
                               changed=list( BP=as.data.frame(go.res) )
                             )

    genetonic[[ comp ]] <- list(
                             res=comp,
                             changed=list(
                               BP=carnation::enrich_to_genetonic(go.res, res_list[[comp]]$res)
                             )
                           )

}

## ----deg_metadata-------------------------------------------------------------
# extract normalized data & metadata
ma <- assay(rld_list[['main']])
colData.i <- colData(rld_list[['main']])

# only keep data from DE genes
idx <- rownames(ma) %in% de.genes[['dex_trt_vs_untrt']][['gene']]
ma.i <- ma[idx,]

# remove any genes with 0 variance
ma.i <- ma.i[rowVars(ma.i) != 0, ]

## ----degpatterns--------------------------------------------------------------
# NOTE: This is the command used to perform pattern analysis
# degpatterns_dex <- DEGreport::degPatterns(
#                      ma.i,
#                      colData.i,
#
#                      time='cell',
#                      col='dex',
#
#                      # NOTE: reduce and merge cutoff----------------------------------------
#                      #   Reduce will merge clusters that are similar; similarity determined
#                      #   by cutoff
#                      reduce=TRUE,
#
#                      plot=FALSE
#                    )

# We use the pre-computed results here instead
data(degpatterns_dex, package='carnation')

## ----deg_extract--------------------------------------------------------------
# extract normalized slot and add symbol column
p_norm <- degpatterns_dex$normalized
p_norm[[ 'SYMBOL' ]] <- anno_df[p_norm[['genes']], 'SYMBOL']

# save pattern analysis results
degpatterns <- list(dex_by_cell=p_norm)

## ----compose------------------------------------------------------------------
combined <- list(res_list=res_list,
                 dds_list=dds_list,
                 rld_list=rld_list,
                 enrich_list=enrich_list,
                 genetonic=genetonic,
                 degpatterns_list=degpatterns)
saveRDS(combined, 'carnation_vignette.rds', compress=FALSE)

## ----first_load---------------------------------------------------------------
library(carnation)

## ----install, eval=FALSE------------------------------------------------------
# install_carnation()  # Installs plotly and kaleido for PDF export

## ----run_app, eval=FALSE------------------------------------------------------
# run_carnation()

## ----run_app2, eval=FALSE-----------------------------------------------------
# run_carnation(options=list(port=12345, launch.browser=FALSE))

## ----server, eval=FALSE-------------------------------------------------------
# # Create user database
# credentials <- data.frame(
#   user = c('shinymanager'),
#   password = c('12345'),
#   admin = c(TRUE),
#   stringsAsFactors = FALSE
# )
# 
# # Initialize the database
# shinymanager::create_db(
#   credentials_data = credentials,
#   sqlite_path = 'credentials.sqlite',
#   passphrase = 'admin_passphrase'
# )

## ----run_server, eval=FALSE---------------------------------------------------
# run_carnation(credentials='credentials.sqlite', passphrase='admin_passphrase')

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