### R code from vignette source 'vignettes/TCC/inst/doc/TCC.Rnw'

###################################################
### code chunk number 1: quick_start_case_1
###################################################
library(TCC)
data(hypoData)
samplesize <- 100
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "bayseq",
                       iteration = 1, samplesize = samplesize) 
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)


###################################################
### code chunk number 2: quick_start_case_2
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)


###################################################
### code chunk number 3: quick_start_case_3
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
tcc <- new("TCC", hypoData[,c(1,4)], group)
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
tcc <- estimateDE(tcc, test.method = "deseq", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)


###################################################
### code chunk number 4: load_hypoData
###################################################
library(TCC)
data(hypoData)
head(hypoData)
dim(hypoData)
group <- c(1, 1, 1, 2, 2, 2)


###################################################
### code chunk number 5: group_argument_setting
###################################################
group <- c(1, 1, 1, 1, 2, 2, 2, 2, 2)


###################################################
### code chunk number 6: construct_tcc
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc


###################################################
### code chunk number 7: summary_tcc
###################################################
head(tcc$count)
tcc$group


###################################################
### code chunk number 8: subset_tcc
###################################################
dim(tcc$count)
tcc.sub1 <- subset(tcc, c(rep(TRUE, 20), rep(FALSE, 980)))
dim(tcc.sub1$count)
tcc.sub2 <- tcc[1:20]
dim(tcc.sub2$count)


###################################################
### code chunk number 9: filter_low_count_genes_tcc
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- filterLowCountGenes(tcc)
dim(tcc$count)


###################################################
### code chunk number 10: filter_low_count_genes
###################################################
filter <- as.logical(rowSums(hypoData) > 0)
dim(hypoData[filter, ])
tcc <- new("TCC", hypoData[filter, ], group)
dim(tcc$count)


###################################################
### code chunk number 11: run_tbt_tcc
###################################################
set.seed(1000)
library(TCC)
data(hypoData)
samplesize <- 100
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "bayseq",
                       iteration = 1, samplesize = samplesize)


###################################################
### code chunk number 12: check_tbt_execution_time
###################################################
tcc$norm.factors
tcc$DEGES$execution.time


###################################################
### code chunk number 13: run_tbt_org
###################################################
set.seed(1000)
library(TCC)
data(hypoData)
samplesize <- 100
group <- c(1, 1, 1, 2, 2, 2)
### STEP 1 ###
d <- DGEList(count = hypoData, group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors
norm.factors <- norm.factors / mean(norm.factors)
### STEP 2 ###
cD <- new("countData", data = hypoData, replicates = group,
         groups = list(NDE = rep(1, length = length(group)), DE = group),
         libsizes = colSums(hypoData) * norm.factors)
cD <- getPriors.NB(cD, samplesize = samplesize, estimation = "QL", cl = NULL)
cD <- getLikelihoods.NB(cD, pET = "BIC", cl = NULL)
is.DEG <- as.logical(rank(-cD@posteriors[, "DE"]) < 
                     (nrow(hypoData) * cD@estProps[2]))
### STEP 3 ###
d <- DGEList(count = hypoData[!is.DEG, ], group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors * colSums(hypoData[!is.DEG, ]) / 
                colSums(hypoData)
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


###################################################
### code chunk number 14: degesEdger_tcc
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors
tcc$DEGES$execution.time


###################################################
### code chunk number 15: degesEdger_edger
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
FDR <- 0.1
floorPDEG <- 0.05
d <- DGEList(counts = hypoData, group = group)
### STEP 1 ###
d <- calcNormFactors(d)
### STEP 2 ###
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
result <- exactTest(d)
q.value <- p.adjust(result$table$PValue, method = "BH")
if (sum(q.value < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(q.value < FDR)
} else {
  is.DEG <- as.logical(rank(result$table$PValue, ties.method = "min") <=
                       nrow(hypoData) * floorPDEG)
}
### STEP 3 ###
d <- DGEList(counts = hypoData[!is.DEG, ], group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors * colSums(hypoData[!is.DEG, ]) /
                  colSums(hypoData)
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


###################################################
### code chunk number 16: idegesEdger_tcc
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors
tcc$DEGES$execution.time


###################################################
### code chunk number 17: degesDeseq_tcc
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors
tcc$DEGES$execution.time


###################################################
### code chunk number 18: degesDeseq_deseq
###################################################
library(TCC)
data(hypoData)  
group <- c(1, 1, 1, 2, 2, 2)
FDR <- 0.1
floorPDEG <- 0.05
cds <- newCountDataSet(hypoData, group)
### STEP 1 ###
cds <- estimateSizeFactors(cds)
### STEP 2 ###
cds <- estimateDispersions(cds)
result <- nbinomTest(cds, 1, 2)
result$pval[is.na(result$pval)] <- 1
result$padj[is.na(result$padj)] <- 1
q.value <- result$padj
if (sum(q.value < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(q.value < FDR)
} else {
  is.DEG <- as.logical(rank(result$pval, ties.method = "min") <=
                       nrow(hypoData) * floorPDEG)
}
### STEP 3 ###
cds <- newCountDataSet(hypoData[!is.DEG, ], group)
cds <- estimateSizeFactors(cds)
norm.factors <- sizeFactors(cds) / colSums(hypoData)
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


###################################################
### code chunk number 19: hypoData_without_replicates
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
tcc <- new("TCC", hypoData[, c(1, 4)], group)
head(tcc$count)
tcc$group


###################################################
### code chunk number 20: degesDeseq_without_replicates_tcc
###################################################
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors


###################################################
### code chunk number 21: degesDeseq_without_replicates_deseq
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
FDR <- 0.1
floorPDEG <- 0.05
cds <- newCountDataSet(hypoData[, c(1, 4)], group)
### STEP 1 ###
cds <- estimateSizeFactors(cds)
### STEP 2 ###
cds <- estimateDispersions(cds, method = "blind", sharingMode = "fit-only")
result <- nbinomTest(cds, 1, 2)
result$pval[is.na(result$pval)] <- 1
result$padj[is.na(result$padj)] <- 1
q.value <- result$padj
if (sum(q.value < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(q.value < FDR)
} else {
  is.DEG <- as.logical(rank(result$pval, ties.method = "min") <=
                       nrow(hypoData) * floorPDEG)
}
### STEP 3 ###
cds <- newCountDataSet(hypoData[!is.DEG, c(1, 4)], group)
cds <- estimateSizeFactors(cds)
norm.factors <- sizeFactors(cds) / colSums(hypoData[, c(1, 4)])
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


###################################################
### code chunk number 22: load_hypoData_mg
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
tcc
dim(tcc$count)


###################################################
### code chunk number 23: degesTbt_multigroup_tcc
###################################################
set.seed(1000)
library(TCC)
data(hypoData_mg)
samplesize <- 100
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "bayseq",
                       iteration = 1, samplesize = samplesize)
tcc$norm.factors


###################################################
### code chunk number 24: degesEdger_multigroup_tcc
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
design <- model.matrix(~ as.factor(group))
coef <- 2:length(unique(group))


###################################################
### code chunk number 25: degesEdger_multigroup_tcc_nf
###################################################
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
tcc$norm.factors


###################################################
### code chunk number 26: degesEdger_multigroup_edger
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
FDR <- 0.1
floorPDEG <- 0.05
design <- model.matrix(~ as.factor(group))
coef <- 2:length(unique(group))
d <- DGEList(counts = hypoData_mg, group = group)
### STEP 1 ###
d <- calcNormFactors(d)
### STEP 2 ###
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)
fit <- glmFit(d, design)
lrt <- glmLRT(fit, coef = coef)
result <- topTags(lrt, n = nrow(hypoData_mg))
result <- result$table[rownames(hypoData_mg), ]
if (sum(result$FDR < FDR) > (floorPDEG * nrow(hypoData_mg))) {
  is.DEG <- as.logical(result$FDR < FDR)
} else  {
  is.DEG <- as.logical(rank(result$PValue, ties.method = "min") <=
                       nrow(hypoData_mg) * floorPDEG)
}
### STEP 3 ###
d <- DGEList(counts = hypoData_mg[!is.DEG, ], group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors * colSums(hypoData_mg[!is.DEG, ]) /
                colSums(hypoData_mg)
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


###################################################
### code chunk number 27: degesDeseq_multigroup_tcc
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
fit1 <- count ~ condition
fit0 <- count ~ 1


###################################################
### code chunk number 28: degesDeseq_multigroup_tcc_nf
###################################################
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 1)
tcc$norm.factors


###################################################
### code chunk number 29: degesDeseq_multigroup_deseq
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
FDR <- 0.1
floorPDEG <- 0.05
tcc <- new("TCC", hypoData_mg, group)
fit1 <- count ~ condition
fit0 <- count ~ 1
cds <- newCountDataSet(hypoData_mg, group)
### STEP 1 ###
cds <- estimateSizeFactors(cds)
### STEP 2 ###
cds <- estimateDispersions(cds)
reduced.model <- fitNbinomGLMs(cds, fit0)
full.model <- fitNbinomGLMs(cds, fit1)
p.value <- nbinomGLMTest(full.model, reduced.model)
p.value[is.na(p.value)] <- 1
q.value <- p.adjust(p.value, method = "BH")
if (sum(q.value < FDR) > (floorPDEG * nrow(hypoData_mg))) {
  is.DEG <- as.logical(q.value < FDR)
} else {
  is.DEG <- as.logical(rank(p.value, ties.method = "min") <=
                       nrow(hypoData_mg) * floorPDEG)
}
### STEP 3 ###
cds <- newCountDataSet(hypoData_mg[!is.DEG, ], group)
cds <- estimateSizeFactors(cds)
norm.factors <- sizeFactors(cds) / colSums(hypoData_mg)
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


###################################################
### code chunk number 30: hypoData_deg_label
###################################################
library(TCC)
data(hypoData)
nonDEG <- 201:1000
summary(hypoData[nonDEG, ])


###################################################
### code chunk number 31: calc_median_of_nondeg_hypoData
###################################################
apply(hypoData[nonDEG, ], 2, median)


###################################################
### code chunk number 32: calc_median_of_nondeg_hypoData_hide
###################################################
hypoData.median <- apply(hypoData[nonDEG, ], 2, median)


###################################################
### code chunk number 33: get_normalized_data
###################################################
normalized.count <- getNormalizedData(tcc)


###################################################
### code chunk number 34: degesEdger_tcc_normed_data
###################################################
library(TCC)
data(hypoData)  
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
normalized.count <- getNormalizedData(tcc)
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 35: degesEdger_edger_normed_data
###################################################
library(TCC)
data(hypoData)
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2)
FDR <- 0.1
floorPDEG <- 0.05
d <- DGEList(counts = hypoData, group = group)
###  Step 1  ###
d <- calcNormFactors(d)
###  Step 2  ###
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
result <- exactTest(d)
q.value <- p.adjust(result$table$PValue, method = "BH")
if (sum(q.value < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(q.value < FDR)
} else {
  is.DEG <- as.logical(order(rank(result$table$PValue)) <= 
                       nrow(hypoData) * floorPDEG)
}
###  Step 3  ###
d <- DGEList(counts = hypoData[!is.DEG, ], group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors * colSums(hypoData[!is.DEG, ]) /
                  colSums(hypoData)
norm.factors <- norm.factors / mean(norm.factors)
effective.libsizes <- colSums(hypoData) * norm.factors
normalized.count <- sweep(hypoData, 2,
                          mean(effective.libsizes) / effective.libsizes, "*")
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 36: tmm_edger_normed_data
###################################################
library(TCC)
data(hypoData)
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2)
d <- DGEList(count = hypoData, group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors
norm.factors <- norm.factors / mean(norm.factors)
effective.libsizes <- colSums(hypoData) * norm.factors
normalized.count <- sweep(hypoData, 2, 
                          mean(effective.libsizes) / effective.libsizes, "*")
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 37: tmm_tcc_normed_data
###################################################
library(TCC)
data(hypoData)
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", iteration = 0)
normalized.count <- getNormalizedData(tcc)
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 38: degesDeseq_tcc_normed_data
###################################################
library(TCC)
data(hypoData)  
group <- c(1, 1, 1, 2, 2, 2)
nonDEG <- 201:1000
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
normalized.count <- getNormalizedData(tcc)
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 39: degesDeseq_deseq_normed_data
###################################################
library(TCC)
data(hypoData)  
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2)
FDR <- 0.1
floorPDEG <- 0.05
cds <- newCountDataSet(hypoData, group)
###  Step 1  ###
cds <- estimateSizeFactors(cds)
###  Step 2  ###
cds <- estimateDispersions(cds)
result <- nbinomTest(cds, 1, 2)
result$pval[is.na(result$pval)] <- 1
result$padj[is.na(result$padj)] <- 1
q.value <- result$padj
if (sum(q.value < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(q.value < FDR)
} else {
  is.DEG <- as.logical(rank(result$pval, ties.method = "min") <=
                       nrow(hypoData) * floorPDEG)
}
###  Step 3  ###
cds <- newCountDataSet(hypoData[!is.DEG, ], group)
cds <- estimateSizeFactors(cds)
norm.factors <- sizeFactors(cds) / colSums(hypoData)
norm.factors <- norm.factors / mean(norm.factors)
effective.libsizes <- colSums(hypoData) * norm.factors
normalized.count <- sweep(hypoData, 2,
                          mean(effective.libsizes) / effective.libsizes, "*")
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 40: degesDeseq_tcc_without_replicates_normed_data
###################################################
library(TCC)
data(hypoData)  
nonDEG <- 201:1000
group <- c(1, 2)
tcc <- new("TCC", hypoData[, c(1, 4)], group)
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
normalized.count <- getNormalizedData(tcc)
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 41: degesDeseq_deseq_without_replicates_normed_data
###################################################
library(TCC)
data(hypoData) 
nonDEG <- 201:1000
group <- c(1, 2)
FDR <- 0.1
floorPDEG <- 0.05
cds <- newCountDataSet(hypoData[,c(1, 4)], group)
###  Step 1  ###
cds <- estimateSizeFactors(cds)
###  Step 2  ###
cds <- estimateDispersions(cds, method = "blind", sharingMode = "fit-only")
result <- nbinomTest(cds, 1, 2)
result$pval[is.na(result$pval)] <- 1
result$padj[is.na(result$padj)] <- 1
q.value <- result$padj
if (sum(q.value < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(q.value < FDR)
} else {
  is.DEG <- as.logical(rank(result$pval, ties.method = "min") <=
                       nrow(hypoData) * floorPDEG)
}
###  Step 3  ###
cds <- newCountDataSet(hypoData[!is.DEG, c(1, 4)], group)
cds <- estimateSizeFactors(cds)
norm.factors <- sizeFactors(cds) / colSums(hypoData[, c(1, 4)])
norm.factors <- norm.factors / mean(norm.factors)
effective.libsizes <- colSums(hypoData[, c(1, 4)]) * norm.factors
normalized.count <- sweep(hypoData[, c(1, 4)], 2,
                          mean(effective.libsizes) / effective.libsizes, "*")
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 42: deseq_deseq_without_replicates_normed_data
###################################################
library(TCC)
data(hypoData)
nonDEG <- 201:1000
group <- c(1, 2)
cds <- newCountDataSet(hypoData[, c(1, 4)], group)
cds <- estimateSizeFactors(cds)
normalized.count <- counts(cds, normalized = TRUE)
apply(normalized.count[nonDEG, ], 2, median)


###################################################
### code chunk number 43: hypoData_mg_deg_label
###################################################
library(TCC)
data(hypoData_mg)
nonDEG <- 201:1000
summary(hypoData_mg[nonDEG, ])


###################################################
### code chunk number 44: calc_median_of_nondeg_hypoData_mg
###################################################
apply(hypoData_mg[nonDEG, ], 2, median)


###################################################
### code chunk number 45: NormHypoDataMultiGGet
###################################################
hypoData_mg.median <- apply(hypoData_mg[nonDEG, ], 2, median)


###################################################
### code chunk number 46: degesEdger_tcc_multigroup_normed_data
###################################################
library(TCC)
data(hypoData_mg)
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
design <- model.matrix(~ as.factor(group))
coef <- 2:length(unique(group))
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 3)
normalized.count <- getNormalizedData(tcc)
apply(normalized.count[nonDEG, ], 2, median)
range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 47: degesEdger_tcc_multigroup_normed_data_hdie
###################################################
normByiDEGES <- range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 48: tmm_tcc_multigroup_normed_data
###################################################
library(TCC)
data(hypoData_mg)
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", iteration = 0)
normalized.count <- getNormalizedData(tcc)
apply(normalized.count[nonDEG, ], 2, median)
range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 49: tmm_tcc_multigroup_normed_data_hide
###################################################
normByTMM <- range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 50: idegesEdger_tcc
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)


###################################################
### code chunk number 51: get_estimated_results
###################################################
result <- getResult(tcc, sort = TRUE)
head(result)


###################################################
### code chunk number 52: summary_of_estimated_deg
###################################################
table(tcc$estimatedDEG) 


###################################################
### code chunk number 53: maplot_of_estimated_result
###################################################
plot(tcc)


###################################################
### code chunk number 54: idegerEdger_tcc_bayseq_summary
###################################################
set.seed(1000)
library(TCC)
data(hypoData)
samplesize <- 100
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
tcc <- estimateDE(tcc, test.method = "bayseq", 
                  FDR = 0.1, samplesize = samplesize)
result <- getResult(tcc, sort = TRUE)
head(result)
table(tcc$estimatedDEG)


###################################################
### code chunk number 55: idegerEdger_tcc_bayseq_org_summary
###################################################
set.seed(1000)
library(TCC)
data(hypoData)
samplesize <- 100
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
effective.libsizes <- colSums(tcc$count) * tcc$norm.factors
groups <- list(NDE = rep(1, length(group)), DE = group)
cD <- new("countData", data = tcc$count, replicates = group,
          libsizes = effective.libsizes, groups = groups)
cD <- getPriors.NB(cD, samplesize = samplesize, 
                   estimation = "QL", cl = NULL)
cD <- getLikelihoods.NB(cD, pET = "BIC", cl = NULL)
tmp <- topCounts(cD, group = "DE", number = nrow(tcc$count))
tmp <- tmp[rownames(tcc$count), ]
p.value <- 1 - tmp$Likelihood
q.value <- tmp$FDR
result <- cbind(p.value, q.value)
rownames(result) <- rownames(tmp)
head(result)


###################################################
### code chunk number 56: idegerDeseq_tcc_without_replicaets_summary
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
tcc <- new("TCC", hypoData[, c(1, 4)], group)
head(tcc$count)
tcc$group
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors
tcc <- estimateDE(tcc, test.method = "deseq", 
                  FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)
table(tcc$estimatedDEG)


###################################################
### code chunk number 57: degerEdegr_tcc_multigroup_summary
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
###  Normalization  ###
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
###  DE analysis  ###
set.seed(1000)
samplesize <- 100
tcc <- estimateDE(tcc, test.method = "bayseq", 
                  FDR = 0.1, samplesize = samplesize)
result <- getResult(tcc, sort = TRUE)
head(result)
table(tcc$estimatedDEG)


###################################################
### code chunk number 58: num_of_fdr_less_than_2
###################################################
sum(result$q.value < 0.2)


###################################################
### code chunk number 59: degerTbt_tcc_multigroup_summary
###################################################
set.seed(1000)
samplesize <- 100
effective.libsizes <- colSums(tcc$count) * tcc$norm.factors
groups <- list(NDE = rep(1, length(group)), DE = group)
cD <- new("countData", data = tcc$count, replicates = group,
          libsizes = effective.libsizes, groups = groups)
cD <- getPriors.NB(cD, samplesize = samplesize, 
                   estimation = "QL", cl = NULL)
cD <- getLikelihoods.NB(cD, pET = "BIC", cl = NULL)
tmp <- topCounts(cD, group = "DE", number = nrow(tcc$count))
tmp <- tmp[rownames(tcc$count), ]
p.value <- 1 - tmp$Likelihood
q.value <- tmp$FDR
result <- cbind(p.value, q.value)
rownames(result) <- rownames(tmp)
head(result)
sum(q.value < 0.1)
sum(q.value < 0.2)


###################################################
### code chunk number 60: degerTbt_tcc_edger_multigroup_summary
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
###  Normalization  ###
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
###  DE analysis  ###
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)
table(tcc$estimatedDEG)


###################################################
### code chunk number 61: degerTbt_tcc_edger_org_multigroup_summary
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
design <- model.matrix(~ as.factor(group))
coef <- 2:length(unique(group))
tcc <- new("TCC", hypoData_mg, group)
###  Normalization  ###
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
###  DE analysis  ###
d <- DGEList(tcc$count, group = group)
d$samples$norm.factors <- tcc$norm.factors
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)
fit <- glmFit(d, design)
lrt <- glmLRT(fit, coef = coef)
tmp <- topTags(lrt, n = nrow(tcc$count))
p.value <- tmp$table$PValue
q.value <- tmp$table$FDR
result <- cbind(p.value, q.value)
rownames(result) <- rownames(tmp)
head(result)
sum(q.value < 0.1)
sum(q.value < 0.2)


###################################################
### code chunk number 62: degerEdger_tcc_edger_org_multigroup_summary_2
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
###  Normalization  ###
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
###  DE analysis  ###
coef <- 3
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1, coef = coef)
result <- getResult(tcc, sort = TRUE)
head(result)
table(tcc$estimatedDEG)


###################################################
### code chunk number 63: degerEdger_tcc_deseq_tcc_multigroup_summary
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
###  Normalization  ###
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
###  DE analysis  ###
fit1 <- count ~ condition
fit0 <- count ~ 1
tcc <- estimateDE(tcc, test.method = "deseq", 
                  FDR = 0.1, fit0 = fit0, fit1 = fit1)
result <- getResult(tcc, sort = TRUE)
head(result)
table(tcc$estimatedDEG)


###################################################
### code chunk number 64: degerEdger_tcc_deseq_org_multigroup_summary
###################################################
library(TCC)
data(hypoData_mg)
group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
tcc <- new("TCC", hypoData_mg, group)
###  Normalization  ###
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
###  DE analysis  ###
fit1 <- count ~ condition
fit0 <- count ~ 1
cds <- newCountDataSet(tcc$count, group)
sizeFactors(cds) <- tcc$norm.factors * colSums(tcc$count)
cds <- estimateDispersions(cds)
reduced.model <- fitNbinomGLMs(cds, fit0)
full.model <- fitNbinomGLMs(cds, fit1)
p.value <- nbinomGLMTest(full.model, reduced.model)
p.value[is.na(p.value)] <- 1
q.value <- p.adjust(p.value, method = "BH")
tmp <- cbind(p.value, q.value)
rownames(tmp) <- tcc$gene_id
result <- tmp[order(p.value), ]
head(result)
sum(q.value < 0.1)
sum(q.value < 0.2)


###################################################
### code chunk number 65: generate_simulation_count_data
###################################################
set.seed(1000)
library(TCC)
tcc <- simulateReadCounts(Ngene = 1000, PDEG = 0.2, 
                              DEG.assign = c(0.9, 0.1),
                              DEG.foldchange = c(4, 4), 
                              replicates = c(3, 3))
dim(tcc$count)
head(tcc$count)
tcc$group


###################################################
### code chunk number 66: str_of_simulation_field
###################################################
str(tcc$simulation)


###################################################
### code chunk number 67: tale_of_simulation_trueDEG
###################################################
table(tcc$simulation$trueDEG)


###################################################
### code chunk number 68: analysis_simulation_data
###################################################
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)


###################################################
### code chunk number 69: calc_AUC_of_simulation_data
###################################################
calcAUCValue(tcc)


###################################################
### code chunk number 70: calc_AUC_of_simulation_data_hide
###################################################
auc.degesedger <- calcAUCValue(tcc)


###################################################
### code chunk number 71: calc_AUC_of_simulation_data_ROC
###################################################
AUC(rocdemo.sca(truth = as.numeric(tcc$simulation$trueDEG != 0),
                data = -tcc$stat$rank))


###################################################
### code chunk number 72: calc_AUC_of_tmm_tcc
###################################################
tcc <- calcNormFactors(tcc, norm.method = "tmm", iteration = 0)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
calcAUCValue(tcc)


###################################################
### code chunk number 73: calc_AUC_of_tmm_org
###################################################
d <- DGEList(counts = tcc$count, group = tcc$group$group)
d <- calcNormFactors(d)
d$samples$norm.factors <- d$samples$norm.factors / 
                          mean(d$samples$norm.factors)
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
result <- exactTest(d)
result$table$PValue[is.na(result$table$PValue)] <- 1
AUC(rocdemo.sca(truth = as.numeric(tcc$simulation$trueDEG != 0),
                data = -rank(result$table$PValue)))


###################################################
### code chunk number 74: calc_AUC_of_degesTbt_tcc_hide
###################################################
set.seed(1000)
samplesize <- 100
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "bayseq",
                       iteration = 1, samplesize = samplesize)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
auc.degestbt <- calcAUCValue(tcc)


###################################################
### code chunk number 75: calc_AUC_of_degesTbt_tcc
###################################################
set.seed(1000)
samplesize <- 100
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "bayseq",
                       iteration = 1, samplesize = samplesize)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
calcAUCValue(tcc)


###################################################
### code chunk number 76: generate_simulation_data_without_replicate
###################################################
set.seed(1000)
library(TCC)
tcc <- simulateReadCounts(Ngene = 1000, PDEG = 0.2, 
                              DEG.assign = c(0.9, 0.1),
                              DEG.foldchange = c(4, 4), 
                              replicates = c(1, 1))
dim(tcc$count)
head(tcc$count)
tcc$group


###################################################
### code chunk number 77: analysis_without_replicate_simulation_data
###################################################
tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
tcc <- estimateDE(tcc, test.method = "deseq")
calcAUCValue(tcc)


###################################################
### code chunk number 78: calc_AUC_of_deseq_tcc
###################################################
tcc <- calcNormFactors(tcc, norm.method = "deseq", iteration = 0)
tcc <- estimateDE(tcc, test.method = "deseq")
calcAUCValue(tcc)


###################################################
### code chunk number 79: calc_AUC_of_degesDeseq_org
###################################################
cds <- newCountDataSet(tcc$count, tcc$group$group)
cds <- estimateSizeFactors(cds)
norm.factors <- sizeFactors(cds) / colSums(tcc$count)
norm.factors <- norm.factors / mean(norm.factors)
sizeFactors(cds) <- colSums(tcc$count) * norm.factors
cds <- estimateDispersions(cds, method="blind", sharingMode="fit-only")
result <- nbinomTest(cds, 1, 2)
result$pval[is.na(result$pval)] <- 1
AUC(rocdemo.sca(truth = as.numeric(tcc$simulation$trueDEG != 0),
                data = -rank(result$pval))) 


###################################################
### code chunk number 80: calc_AUC_of_deseq_org
###################################################
cds <- newCountDataSet(tcc$count, tcc$group$group)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds, method="blind", sharingMode="fit-only")
result <- nbinomTest(cds, 1, 2)
result$pval[is.na(result$pval)] <- 1
AUC(rocdemo.sca(truth = as.numeric(tcc$simulation$trueDEG != 0),
                data = -rank(result$pval))) 


###################################################
### code chunk number 81: generate_simulation_data_multigroup
###################################################
set.seed(1000)
library(TCC)
tcc <- simulateReadCounts(Ngene = 1000, PDEG = 0.3,
                         DEG.assign = c(0.7, 0.2, 0.1),
                         DEG.foldchange = c(3, 10, 6),
                         replicates = c(2, 4, 3))
dim(tcc$count)
tcc$group
head(tcc$count)


###################################################
### code chunk number 82: plot_fc_pseudocolor_multigroup
###################################################
plotFCPseudocolor(tcc)


###################################################
### code chunk number 83: degesEdger_edger_tcc_multigroup
###################################################
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
calcAUCValue(tcc)


###################################################
### code chunk number 84: tmm_edger_tcc_multigroup
###################################################
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 0)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
calcAUCValue(tcc)


###################################################
### code chunk number 85: plot_fc_pseudocolor_multigroup_2
###################################################
set.seed(1000)
library(TCC)
tcc <- simulateReadCounts(Ngene = 10000, PDEG = 0.34,
               DEG.assign = c(0.1, 0.3, 0.05, 0.1, 0.05, 0.21, 0.09, 0.1),
               DEG.foldchange = c(3.1, 13, 2, 1.5, 9, 5.6, 4, 2),
               replicates = c(1, 1, 2, 1, 1, 1, 1, 1))
dim(tcc$count)
tcc$group
head(tcc$count)
plotFCPseudocolor(tcc)


###################################################
### code chunk number 86: generate_simulation_data_multifactor_group
###################################################
group <- data.frame(
    GROUP = c("WT", "WT", "WT", "WT", "KO", "KO", "KO", "KO"),
    TIME = c("1d", "1d", "2d", "2d", "1d", "1d", "2d", "2d")
)


###################################################
### code chunk number 87: generate_simulation_data_multifactor_fc
###################################################
DEG.foldchange <- data.frame(
    FACTOR1.1 = c(2, 2, 2, 2, 1, 1, 1, 1),
    FACTOR1.2 = c(1, 1, 1, 1, 3, 3, 3, 3),
    FACTOR2 = c(1, 1, 0.5, 0.5, 1, 1, 4, 4)
)


###################################################
### code chunk number 88: generate_simulation_data_multifactor
###################################################
set.seed(1000)
tcc <- simulateReadCounts(Ngene = 10000, PDEG = 0.2,
               DEG.assign = c(0.5, 0.2, 0.3),
               DEG.foldchange = DEG.foldchange,
               group = group)


###################################################
### code chunk number 89: plot_fc_pseudocolor_multifactor
###################################################
head(tcc$count)
tcc$group
plotFCPseudocolor(tcc)


###################################################
### code chunk number 90: simulate_data_1000
###################################################
set.seed(1000)
library(TCC)
tcc <- simulateReadCounts(Ngene = 20000, PDEG = 0.30, 
                              DEG.assign = c(0.85, 0.15),
                              DEG.foldchange = c(8, 16), 
                              replicates = c(2, 2))
head(tcc$count)


###################################################
### code chunk number 91: maplot_simulate_data_1000
###################################################
plot(tcc)


###################################################
### code chunk number 92: norm_simulate_data_1000
###################################################
normalized.count <- getNormalizedData(tcc)
colSums(normalized.count)
colSums(tcc$count)
mean(colSums(tcc$count))


###################################################
### code chunk number 93: maplot_normed_simulate_data_1000_hide
###################################################
xy <- plot(tcc)
isnot.na <- as.logical(xy[, 1] != min(xy[, 1]))
median.G1 <- median(xy[(tcc$simulation$trueDEG == 1) & isnot.na, 2])
median.G2 <- median(xy[(tcc$simulation$trueDEG == 2) & isnot.na, 2])
median.nonDEG <- median(xy[(tcc$simulation$trueDEG == 0) & isnot.na, 2])


###################################################
### code chunk number 94: maplot_normed_simulate_data_1000
###################################################
plot(tcc, median.lines = TRUE) 


###################################################
### code chunk number 95: maplot_normed_simulate_data_1000_2_hide
###################################################
tcc <- calcNormFactors(tcc, "tmm", "edger", iteration = 3, 
                       FDR = 0.1, floorPDEG = 0.05)
xy <- plot(tcc)
isnot.na <- as.logical(xy[, 1] != min(xy[, 1]))
median.nonDEG <- median(xy[(tcc$simulation$trueDEG == 0) & isnot.na, 2])


###################################################
### code chunk number 96: maplot_normed_simulate_data_1000_2
###################################################
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger", 
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05)
plot(tcc, median.line = TRUE)


###################################################
### code chunk number 97: session_info
###################################################
sessionInfo()