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

## ----setup, message=F, warning=F----------------------------------------------
#remotes::install_github("Divo-Lee/SIEVEseq")

## -----------------------------------------------------------------------------
library(SIEVEseq)

## -----------------------------------------------------------------------------
data("clrCounts3") 
#500 genes, 200 samples per group, differential variability for the first 50 genes,
#CLR-transformed counts table
dim(clrCounts3)
clrCounts3[1:5, c(1:3, 201:203)]

## -----------------------------------------------------------------------------
SN.plot(clrCounts3[2, 1:200]) # gene 2 in control group
SN.plot(clrCounts3[2, 201:400]) # gene 2 in case group 

## -----------------------------------------------------------------------------
clr.SN.fit(clrCounts3[2, 1:200]) # gene 2 in control group
clr.SN.fit(clrCounts3[3:4, 201:400]) # gene 3 and gene 4 in case group

## -----------------------------------------------------------------------------
data("clrCounts3")
 #CLR-transformed counts table, 500 genes, 200 samples per group, 
 #differential expression for the first 50 genes
data("clrCounts2")
 #CLR-transformed counts table, 500 genes, 200 samples per group, 
 #differential variability for the first 50 genes,
dim(clrCounts3); dim(clrCounts2)
 #CLR-transformed counts table, 500 genes, 200 samples per group, 
 #differential expression for the first 50 genes,
groups <- c(rep(0,200), rep(1,200))
 # 200 control samples, and 200 case samples
clrseq_result1 <-  clrSeq(clrCounts3, group = groups) # DE dataset
clrseq_result2 <-  clrSeq(clrCounts2, group = groups) # DV dataset

## -----------------------------------------------------------------------------
head(clrseq_result1, 3) # MLE, DE genes
 #
tail(clrseq_result1, 3) # MLE, non-DE genes
 #

## -----------------------------------------------------------------------------
head(clrseq_result2, 3) # MLE, DV genes
 #
tail(clrseq_result2, 3) # MLE, non-DV genes
 #

## -----------------------------------------------------------------------------
sieve_try1 <- clrSIEVE(clrSeq_result = clrseq_result1,
                       alpha_level = 0.05,
                       order_DE = FALSE,
                       order_LFC = FALSE,
                       order_DS = FALSE,
                       order_sieve = FALSE)
names(sieve_try1)

## -----------------------------------------------------------------------------
DE_test_result1 <- sieve_try1$clrDE_test # results of DE tests
head(DE_test_result1, 3) # DE genes
tail(DE_test_result1, 3) # non-DE genes

## -----------------------------------------------------------------------------
sieve_try2 <- clrSIEVE(clrSeq_result = clrseq_result2,
                       alpha_level = 0.05,
                       order_DE = FALSE,
                       order_LFC = FALSE,
                       order_DS = FALSE,
                       order_sieve = FALSE)
names(sieve_try1)

DV_test_result2 <- sieve_try2$clrDV_test
head(DV_test_result2, 3)
tail(DV_test_result2, 3)

## -----------------------------------------------------------------------------
DS_test_result3 <- sieve_try2$clrDS_test
head(DS_test_result3, 3)

## -----------------------------------------------------------------------------
SIEVE_results <- sieve_try1$clrSIEVE_tests
head(SIEVE_results, 3)

## -----------------------------------------------------------------------------
violin.plot.SIEVE(data = clrCounts3, 
                  "gene1",
                  group = groups,
                  group.names = c("control", "case")) # DE gene (non-DV and non-DS)
clrseq_result1[1,] # MLE, gene1 of clrCounts3. group 1: control; group 2: case

## ----message=F, warning=F-----------------------------------------------------
#install.packages("compositions")
#library(compositions) # a package for compositional data analysis
# clr-transformation
clr.transform <- function(data = NULL){
  # data: count table, genes in rows and samples in columns
  data[data == 0] <- 1/2 
  # A pseudo count 0.5 is added if the count is zero
  clr.count <- t(clr(t(data)))
  clr.count <- matrix(as.numeric(clr.count),
                      nrow = dim(data)[1],
                      ncol = dim(data)[2])
  row.names(clr.count) <- row.names(data)
  return(clr.count)
}