## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse=TRUE, comment="#>", fig.align="center")

suppressPackageStartupMessages({
    library(BiocStyle)
    library(knitr)
    library(RNAshapeQC)
    library(SummarizedExperiment)
    library(ConsensusClusterPlus)
})

## Test R functions
# devtools::document()
# devtools::load_all() # load development version

## ----LoadPkg, eval=FALSE------------------------------------------------------
# ## Installation and loading
# if (!requireNamespace("BiocManager", quietly=TRUE)) {
#     install.packages("BiocManager")
# }
# BiocManager::install("RNAshapeQC")
# 
# library(RNAshapeQC)

## ----construct_pileup, eval=FALSE---------------------------------------------
# construct_pileup(
#     Gene          = "Gene001",
#     studylist     = "TOY",
#     regionsFile   = "/path/to/hg38.gene.regions.txt",
#     regionsFormat = "gencode.regions",
#     regionsCol    = 2,
#     bamFilesList  = list(TOY=toy_bam_paths),
#     caseIDList    = list(TOY=toy_case_ids),
#     outFile       = "pergene/TOY_Gene001_pileup.RData"
# )

## ----LoadData-----------------------------------------------------------------
data("TOY_mrna_se")
data("TOY_total_se")

## ----Str----------------------------------------------------------------------
TOY_mrna_se
TOY_total_se

## ----get_DIIwt----------------------------------------------------------------
## Using SE input directly
res_dii_se <- get_DIIwt(TOY_mrna_se)

names(res_dii_se)
head(res_dii_se$ds.vec)

## (Optional) Using matrix/vector input
data(TOY_mrna_mat)

res_dii_mat <- get_DIIwt(
    DR         = TOY_mrna_mat$DR, 
    TPM        = TOY_mrna_mat$TPM, 
    genelength = TOY_mrna_mat$gene_length
)

## Compare SE and matrix inputs
all.equal(res_dii_se$DR2, res_dii_mat$DR2)
all.equal(res_dii_se$ds.vec, res_dii_mat$ds.vec)
all.equal(res_dii_se$gene.df, res_dii_mat$gene.df)
all.equal(res_dii_se$s, res_dii_mat$s)

## ----plot_DIIwt, message=FALSE, echo=TRUE-------------------------------------
## Path to a temporary file
outfile <- file.path(tempdir(), "Fig_DIIwt.png")

plot_DIIwt(
    DR        = res_dii_se$DR2,
    DIIresult = res_dii_se,
    outFile   = outfile
)

## Saving DIIwt plot to: /.../Fig_DIIwt.png

## ----Fig_DIIwt, out.width="100%",out.height="100%", fig.align="center", echo=FALSE----
knitr::include_graphics(outfile)

## ----get_SOI------------------------------------------------------------------
## Using SE input directly
## The cutoff value can be chosen by inspecting the PD distribution.
res_soi_se <- get_SOI(TOY_total_se, cutoff=1)

names(res_soi_se)
head(res_soi_se$auc.vec)

## (Optional) Using matrix input
data(TOY_total_mat)

res_soi_mat <- get_SOI(
    MCD    = TOY_total_mat$MCD,
    wCV    = TOY_total_mat$wCV,
    cutoff = 1
)

## Compare SE and matrix inputs
all.equal(res_soi_se$auc.vec, res_soi_mat$auc.vec)
all.equal(res_soi_se$auc.coord, res_soi_mat$auc.coord)
all.equal(res_soi_se$rangeMCD, res_soi_mat$rangeMCD)

## ----plot_SOI, message=FALSE, echo=TRUE---------------------------------------
outfile <- file.path(tempdir(), "Fig_SOI.png")

## The cutoff value can be chosen by inspecting the PD distribution.
plot_SOI(
    SOIresult = res_soi_se,
    cutoff    = 1,
    outFile   = outfile 
)

## Saving SOI plot to: /.../Fig_SOI.png

## ----Fig_SOI, out.width="100%",out.height="100%", fig.align="center", echo=FALSE----
knitr::include_graphics(outfile)

## ----se_mrna, message=FALSE, echo=TRUE----------------------------------------
library(SummarizedExperiment)

## Store QC results in mRNA-seq
se_mrna <- SummarizedExperiment(
    assays   = list(DR=assay(TOY_mrna_se, "DR"), TPM=assay(TOY_mrna_se, "TPM")),
    colData  = DataFrame(res_dii_se$ds.vec),
    metadata = list(s=res_dii_se$s)
)
se_mrna

## Subsetting to keep high-quality samples in mRNA-seq
se_mrna_hq <- se_mrna[, se_mrna$DII == "Intact"]
head(assay(se_mrna_hq, "TPM"))

## ----CC, message=FALSE, echo=TRUE---------------------------------------------
library(ConsensusClusterPlus)

## TPM log transform
tpm_hq  <- log2(assay(se_mrna_hq, "TPM") + 1)

## Set a random seed for reproducibility
set.seed(1)

cc <- ConsensusClusterPlus(
    tpm_hq,
    maxK       = 6,
    reps       = 100,
    pItem      = 0.8,
    pFeature   = 0.8,
    clusterAlg = "hc",
    distance   = "pearson",
    plot       = "png",
    title      = tempdir()
)
cc[[3]]$consensusClass

## ----se_total, message=FALSE, echo=TRUE---------------------------------------
## Store QC results in total RNA-seq
se_total <- SummarizedExperiment(
    assays   = list(MCD=assay(TOY_total_se, "MCD"), wCV=assay(TOY_total_se, "wCV")),
    rowData  = DataFrame(Gene=rownames(assay(TOY_total_se, "MCD"))),
    colData  = DataFrame(res_soi_se$auc.vec),
    metadata = list(auc.coord=res_soi_se$auc.coord, rangeMCD=res_soi_se$rangeMCD)
)
se_total

## Subsetting to keep high-quality samples in total RNA-seq
se_total_hq <- se_total[, se_total$SOI == "Optimal"]
head(assay(se_total_hq, "wCV"))

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