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

## ----install------------------------------------------------------------------
# library(devtools)
# install_github("fusk-kpu/SRscore", build_vignettes = TRUE)
library(SRscore)

## ----metadata-----------------------------------------------------------------
library(tibble)
data("MetadataABA")
tibble(MetadataABA)

## ----transcriptome------------------------------------------------------------
data("TranscriptomeABA")
tibble(TranscriptomeABA)

## ----srga---------------------------------------------------------------------
data("SRGA")
tibble(SRGA)

## ----comparison---------------------------------------------------------------
grp <- "Series"
var1 <- "control_sample"
var2 <- "treated_sample"

ebg <- expand_by_group(
  .data = MetadataABA,
  grp = grp,
  var1 = var1,
  var2 = var2
)

unique_series <- unique(MetadataABA$Series)
unique_series

lapply(unique_series, function(x) subset(ebg, Series == x))

## ----SRratio------------------------------------------------------------------
SRratio <- calcSRratio(
  .data = TranscriptomeABA,
  var1 = var1,
  var2 = var2,
  pair = ebg,
  is.log2 = TRUE
)

tibble(SRratio)

## ----conventional SRratio-----------------------------------------------------
conventional_SRratio <- calcSRratio(
  TranscriptomeABA,
  var1 = var1,
  var2 = var2,
  pair = MetadataABA,
  is.log2 = TRUE
)

tibble(conventional_SRratio)

## ----SRscore------------------------------------------------------------------
SRscore <- calcSRscore(SRratio, threshold = c(-2, 2))
head(SRscore)
tibble(SRscore)

## -----------------------------------------------------------------------------
plot_SRscore_distr(SRscore)
plot_SRscore_distr(SRscore, log = TRUE)

## -----------------------------------------------------------------------------
plot_SRscore_rank(SRscore)

## -----------------------------------------------------------------------------
plot_SRscore_top(SRscore)
plot_SRscore_top(SRscore, top_n = 30)

## ----All in one---------------------------------------------------------------
res <- directly_calcSRscore(
  .data1 = MetadataABA,
  grp = grp,
  var1 = var1,
  var2 = var2,
  .data2 = TranscriptomeABA,
  is.log2 = TRUE,
  threshold = c(-2, 2)
)

names(res)
tibble(res$SRscore)

## ----find_diffexp1------------------------------------------------------------
set.seed(1)
res <- find_diffexp(
  sample(SRratio$ensembl_gene_id, 1),
  SRratio,
  SRscore,
  MetadataABA
)

tibble(res$result)
tibble(res$SRscore)

## ----find_diffexp2------------------------------------------------------------
set.seed(1)
res2 <- find_diffexp(
  sample(SRratio$ensembl_gene_id, 10),
  SRratio,
  SRscore,
  MetadataABA
)

tibble(res2$result)
tibble(res2$SRscore)

## ----ea, eval=FALSE-----------------------------------------------------------
# library(clusterProfiler)
# library(ggplot2)
# 
# ego <- enrichGO(
#   gene = SRscore$ensembl_gene_id[SRscore$score >= 1],
#   OrgDb = "org.At.tair.db",
#   keyType = "TAIR",
#   ont = "BP",
#   maxGSSize = 2000
# )
# 
# dotplot(ego, showCategory = 5, font.size = 14) +
#   theme(text = element_text(size = 14))

## ----heatmap------------------------------------------------------------------
library(ComplexHeatmap)
library(RColorBrewer)

cor_breaks <- seq(-2, 2, length.out = 51)
cor_color <- colorRampPalette(c("blue", "white", "red"))(51)

annotation_row <- res2$result[, c("treatment", "tissue")]
pal_treatment <- brewer.pal(length(unique(annotation_row$treatment)), "Set1")
pal_tissue <- brewer.pal(length(unique(annotation_row$tissue)), "Set2")

names(pal_treatment) <- unique(annotation_row$treatment)
names(pal_tissue) <- unique(annotation_row$tissue)

ComplexHeatmap::pheatmap(
  as.matrix(res2$result[, sapply(res2$result, is.numeric)]),
  breaks = cor_breaks,
  color = cor_color,
  cluster_rows = FALSE,
  name = "SRratio",
  annotation_row = annotation_row,
  annotation_colors = list(
    treatment = pal_treatment,
    tissue = pal_tissue
    )
)

## ----tm-----------------------------------------------------------------------
library(genefilter)
library(DT)

cl <- colnames(Filter(is.numeric, SRGA))
df <- as.matrix(column_to_rownames(SRGA, var = "ensembl_gene_id")[cl])

template <- "AT1G09350"

close_genes <- genefinder(
  df,
  ilist = template,
  numResults = 5,
  method = "euclidean"
)

datatable(
  SRGA[SRGA$ensembl_gene_id == template, ],
  options = list(dom = "lrtBip"),
  rownames = FALSE
)

datatable(
  SRGA[close_genes[[1]]$indices, ],
  options = list(dom = "lrtBip"),
  rownames = FALSE,
)

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