inst/doc/TCC.R

In TCC: TCC: Differential expression analysis for tag count data with robust normalization strategies

### R code from vignette source 'TCC.Rnw'

###################################################
### code chunk number 1: 3408021901
###################################################
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 2: 4301933845
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
tcc <- new("TCC", hypoData[,c(1,4)], 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: 4378919285
###################################################
library(TCC)
data(hypoData)
head(hypoData)
dim(hypoData)
group <- c(1, 1, 1, 2, 2, 2)


###################################################
### code chunk number 4: 1209835923
###################################################
group <- c(1, 1, 1, 1, 2, 2, 2, 2, 2)


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


###################################################
### code chunk number 6: 5048219812
###################################################
head(tcc$count)
tcc$group


###################################################
### code chunk number 7: 4085249378
###################################################
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 8: 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 9: 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 10: check_tbt_execution_time
###################################################
tcc$norm.factors
tcc$DEGES$execution.time


###################################################
### code chunk number 11: run_tbt_org
###################################################
set.seed(1000)
library(TCC)
data(hypoData)
samplesize <- 100
group <- c(1, 1, 1, 2, 2, 2)
floorPDEG <- 0.05
FDR <- 0.1
### 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@annotation <- data.frame(rowID = 1:nrow(hypoData))
cD <- getPriors.NB(cD, samplesize = samplesize, estimation = "QL", cl = NULL)
cD <- getLikelihoods(cD, pET = "BIC", cl = NULL)
result <- topCounts(cD, group = "DE", number = nrow(hypoData))
result <- result[order(result$rowID), ]
pval <- result$FWER.DE
qval <- result$FDR.DE
if (sum(qval < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(qval < FDR)
} else {
  is.DEG <- as.logical(rank(pval, ties.method = "min") <= 
                       nrow(hypoData) * floorPDEG)
}
### 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 12: 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 13: 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)
pval <- result$table$PValue
qval <- p.adjust(pval, method = "BH")
if (sum(qval < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(qval < FDR)
} else {
  is.DEG <- as.logical(rank(pval, 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 14: 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 15: 4390811048
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "deseq2", test.method = "deseq2",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors
tcc$DEGES$execution.time


###################################################
### code chunk number 16: 0869034832
###################################################
library(TCC)
data(hypoData)
FDR <- 0.1
floorPDEG <- 0.05
group <- factor(c(1, 1, 1, 2, 2, 2))
cl <- data.frame(group = group)
design <- formula(~ group)
dds <- DESeqDataSetFromMatrix(countData = hypoData, colData = cl,
                              design = design)
### STEP 1 ###
dds <- estimateSizeFactors(dds)
sizeFactors(dds) <- sizeFactors(dds) / mean(sizeFactors(dds))

### STEP 2 ###
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
result <- results(dds)
result$pvalue[is.na(result$pvalue)] <- 1
pval <- result$pvalue
qval <- p.adjust(pval, method = "BH")
if (sum(qval < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(qval < FDR)
} else {
  is.DEG <- as.logical(rank(pval, ties.method = "min") <=
                       nrow(hypoData) * floorPDEG)
}
### STEP 3 ###
dds <- DESeqDataSetFromMatrix(countData = hypoData[!is.DEG, ], colData = cl,
                              design = design)
dds <- estimateSizeFactors(dds)
norm.factors <- sizeFactors(dds) / colSums(hypoData)
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


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


###################################################
### code chunk number 18: 4352419012
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
tcc <- new("TCC", hypoData[, c(1, 4)], group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
iteration = 1, FDR = 0.1, floorPDEG = 0.05)
tcc$norm.factors


###################################################
### code chunk number 19: 5670112812
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
FDR <- 0.1
floorPDEG <- 0.05
d <- DGEList(counts = hypoData[, c(1, 4)], group = group)
### STEP 1 ###
d <- calcNormFactors(d)
### STEP 2 ###
d <- estimateGLMCommonDisp(d, method = "deviance", robust = TRUE, subset = NULL)
result <- exactTest(d)
pval <- result$table$PValue
qval <- p.adjust(pval, method = "BH")
if (sum(qval < FDR) > (floorPDEG * nrow(hypoData))) {
  is.DEG <- as.logical(qval < FDR)
} else {
  is.DEG <- as.logical(rank(pval, ties.method = "min") <=
                       nrow(hypoData) * floorPDEG)
}
### STEP 3 ###
d <- DGEList(counts = hypoData[!is.DEG, c(1, 4)], group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors * 
                colSums(hypoData[!is.DEG, c(1, 4)]) / colSums(hypoData[, c(1, 4)])
norm.factors <- norm.factors / mean(norm.factors)
norm.factors


###################################################
### code chunk number 20: 9347533928
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
head(hypoData)


###################################################
### code chunk number 21: 2394782901
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
group <- c(1, 1, 1, 2, 2, 2)
pair <- c(1, 2, 3, 1, 2, 3)
cl <- data.frame(group = group, pair = pair)
tcc <- new("TCC", hypoData, cl)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05, paired = TRUE)
tcc$norm.factors


###################################################
### code chunk number 22: 8939487592
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
group <- factor(c(1, 1, 1, 2, 2, 2))
pair <- factor(c(1, 2, 3, 1, 2, 3))
design <- model.matrix(~ group + pair)
coef <- 2
FDR <- 0.1
floorPDEG <- 0.05
d <- DGEList(counts = hypoData, 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)
pval <- lrt$table$PValue
qval <- p.adjust(pval, method = "BH")
if (sum(qval < FDR) > (floorPDEG * nrow(hypoData))){
  is.DEG <- as.logical(qval < FDR)
} else {
  is.DEG <- as.logical(rank(pval, 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 23: 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 24: 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 25: 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 26: degesEdger_multigroup_tcc_nf
###################################################
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1, design = design, coef = coef)
tcc$norm.factors


###################################################
### code chunk number 27: 3498028981
###################################################
library(TCC)
data(hypoData_mg)
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 = "edger",
                       iteration = 1)
tcc$norm.factors


###################################################
### code chunk number 28: 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 <- as.data.frame(topTags(lrt, n = nrow(hypoData_mg)))
result <- result$table[rownames(hypoData_mg), ]
pval <- lrt$table$PValue
qval <- p.adjust(pval, method = "BH")
if (sum(qval < FDR) > (floorPDEG * nrow(hypoData_mg))) {
  is.DEG <- as.logical(qval < FDR)
} else  {
  is.DEG <- as.logical(rank(pval, 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 29: hypoData_deg_label
###################################################
library(TCC)
data(hypoData)
nonDEG <- 201:1000
summary(hypoData[nonDEG, ])


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


###################################################
### code chunk number 31: calc_median_of_nondeg_hypoData_hide
###################################################
hypoData.median <- apply(hypoData[nonDEG, ], 2, median)
hypoData.14.median <- apply(hypoData[nonDEG, c(1, 4)], 2, median)


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


###################################################
### code chunk number 33: 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)
range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 34: 4390011292
###################################################
buff_1 <- range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 35: 2391104042
###################################################
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
effective.libsizes <- colSums(hypoData) * norm.factors
normalized.count <- sweep(hypoData, 2,
                    mean(effective.libsizes) / effective.libsizes, "*")
apply(normalized.count[nonDEG, ], 2, median)
range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 36: 4330011292
###################################################
buff_2 <- range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 37: 5647309237
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
nonDEG <- 201:1000
group <- c(1, 1, 1, 2, 2, 2)
pair <- c(1, 2, 3, 1, 2, 3)
cl <- data.frame(group = group, pair = pair)
tcc <- new("TCC", hypoData, cl)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 1, FDR = 0.1, floorPDEG = 0.05, paired = TRUE)
normalized.count <- getNormalizedData(tcc)
head(normalized.count, n = 4)
apply(normalized.count[nonDEG, ], 2, median)
range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 38: 1674110292
###################################################
buff_1 <- range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 39: 4387803948
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
nonDEG <- 201:1000
group <- factor(c(1, 1, 1, 2, 2, 2))
d <- DGEList(counts = hypoData, group = group)
d <- calcNormFactors(d)
norm.factors <- d$samples$norm.factors
effective.libsizes <- colSums(hypoData) * norm.factors
normalized.count <- sweep(hypoData, 2,
                          mean(effective.libsizes) / effective.libsizes, "*")
apply(normalized.count[nonDEG, ], 2, median)
range(apply(normalized.count[nonDEG, ], 2, median))


###################################################
### code chunk number 40: 1674110992
###################################################
buff_2 <- range(apply(normalized.count[nonDEG, ], 2, median))


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


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


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


###################################################
### code chunk number 44: 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 45: degesEdger_tcc_multigroup_normed_data_hdie
###################################################
normByiDEGES <- range(apply(normalized.count[nonDEG, ], 2, median))


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


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


###################################################
### code chunk number 48: 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 49: get_estimated_results
###################################################
result <- getResult(tcc, sort = TRUE)
head(result)


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


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


###################################################
### code chunk number 52: 3498757592
###################################################
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


###################################################
### code chunk number 53: 5901287562
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
d <- DGEList(count = hypoData, group = group)
d <- calcNormFactors(d)
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
out <- exactTest(d)
result <- as.data.frame(topTags(out, n = nrow(hypoData), sort.by = "PValue"))
head(result)


###################################################
### code chunk number 54: 1027518347
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", iteration = 0)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


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


###################################################
### code chunk number 56: 7512913498
###################################################
table(tcc$estimatedDEG)


###################################################
### code chunk number 57: 7512013498
###################################################
tb <- table(tcc$estimatedDEG)


###################################################
### code chunk number 58: 3489400103
###################################################
plot(tcc)


###################################################
### code chunk number 59: 4012399910
###################################################
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


###################################################
### code chunk number 60: 4930011190
###################################################
library(TCC)
data(hypoData)
group <- factor(c(1, 1, 1, 2, 2, 2))
cl <- data.frame(group = group)
design <- formula(~ group)
dds <- DESeqDataSetFromMatrix(countData = hypoData, colData = cl,
                              design = design)
dds <- estimateSizeFactors(dds)
sizeFactors(dds) <- sizeFactors(dds) / mean(sizeFactors(dds))
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
result <- results(dds)
head(result[order(result$pvalue),])
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
result$pvalue[is.na(result$pvalue)] <- 1
AUC(rocdemo.sca(truth = truth, data = -rank(result$pvalue)))


###################################################
### code chunk number 61: 4390023901
###################################################
library(TCC)
data(hypoData)
group <- c(1, 1, 1, 2, 2, 2)
tcc <- new("TCC", hypoData, group)
tcc <- calcNormFactors(tcc, norm.method = "deseq2", iteration = 0)
tcc <- estimateDE(tcc, test.method = "deseq2", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))
tcc$norm.factors


###################################################
### code chunk number 62: 4090911011
###################################################
buff_1 <- AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


###################################################
### code chunk number 63: 0309001992
###################################################
library(TCC)
data(hypoData)
group <- factor(c(1, 1, 1, 2, 2, 2))
cl <- data.frame(group = group)
design <- formula(~ group)
dds <- DESeqDataSetFromMatrix(countData = hypoData, colData = cl,
                              design = design)
dds <- estimateSizeFactors(dds)
size.factors <- sizeFactors(dds)
size.factors
hoge <- size.factors / colSums(hypoData)
norm.factors <- hoge / mean(hoge)
norm.factors
ef.libsizes <- norm.factors * colSums(hypoData)
ef.libsizes


###################################################
### code chunk number 64: 4328593702
###################################################
library(TCC)
data(hypoData)
group <- factor(c(1, 1, 1, 2, 2, 2))
cl <- data.frame(group = group)
design <- formula(~ group)
dds <- DESeqDataSetFromMatrix(countData = hypoData, colData = cl,
                              design = design)
dds <- estimateSizeFactors(dds)
sizeFactors(dds) <- sizeFactors(dds) / mean(sizeFactors(dds))
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
result <- results(dds)
head(result[order(result$pvalue),])
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
result$pvalue[is.na(result$pvalue)] <- 1
AUC(rocdemo.sca(truth = truth, data = -rank(result$pvalue)))


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


###################################################
### code chunk number 66: 8439100943
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
tcc <- new("TCC", hypoData[, c(1, 4)], 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 67: 4390087612
###################################################
result <- getResult(tcc, sort = TRUE)
head(result)


###################################################
### code chunk number 68: 0238735032
###################################################
table(tcc$estimatedDEG)


###################################################
### code chunk number 69: 0233731032
###################################################
tb <- table(tcc$estimatedDEG)


###################################################
### code chunk number 70: 3489320103
###################################################
plot(tcc)


###################################################
### code chunk number 71: 5702784221
###################################################
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


###################################################
### code chunk number 72: 4890357892
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
d <- DGEList(counts = hypoData[, c(1, 4)], group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, method = "deviance", robust = TRUE, subset = NULL)
out <- exactTest(d)
result <- as.data.frame(topTags(out, n = nrow(hypoData), sort.by = "PValue"))
head(result)


###################################################
### code chunk number 73: 4741957232
###################################################
library(TCC)
data(hypoData)
group <- c(1, 2)
tcc <- new("TCC", hypoData[, c(1, 4)], group)
tcc <- calcNormFactors(tcc, norm.method = "tmm", iteration = 0)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1)
result <- getResult(tcc, sort = TRUE)
head(result)
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


###################################################
### code chunk number 74: 0285012872
###################################################
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
head(hypoData)


###################################################
### code chunk number 75: 4891209302
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
group <- c(1, 1, 1, 2, 2, 2)
pair <- c(1, 2, 3, 1, 2, 3)
cl <- data.frame(group = group, pair = pair)
tcc <- new("TCC", hypoData, cl)
tcc <- calcNormFactors(tcc, norm.method = "tmm", test.method = "edger",
                       iteration = 3, FDR = 0.1, floorPDEG = 0.05, paired = TRUE)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1, paired = TRUE)


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


###################################################
### code chunk number 77: 3909884931
###################################################
table(tcc$estimatedDEG)


###################################################
### code chunk number 78: 3934884932
###################################################
tb <- table(tcc$estimatedDEG)


###################################################
### code chunk number 79: 3482340103
###################################################
plot(tcc)


###################################################
### code chunk number 80: 3490930192
###################################################
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


###################################################
### code chunk number 81: 8402288128
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
group <- factor(c(1, 1, 1, 2, 2, 2))
pair <- factor(c(1, 2, 3, 1, 2, 3))
design <- model.matrix(~ group + pair)
coef <- 2
d <- DGEList(counts = hypoData, group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)
fit <- glmFit(d, design)
lrt <- glmLRT(fit, coef = coef)
topTags(lrt, n = 6)
p.value <- lrt$table$PValue
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -rank(p.value)))


###################################################
### code chunk number 82: 4209576561
###################################################
library(TCC)
data(hypoData)
colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC",
                        "N_dogA", "N_dogB", "N_dogC")
group <- c(1, 1, 1, 2, 2, 2)
pair <- c(1, 2, 3, 1, 2, 3)
cl <- data.frame(group = group, pair = pair)
tcc <- new("TCC", hypoData, cl)
tcc <- calcNormFactors(tcc, norm.method = "tmm", iteration = 0, paired = TRUE)
tcc <- estimateDE(tcc, test.method = "edger", FDR = 0.1, paired = TRUE)
result <- getResult(tcc, sort = TRUE)
head(result)
library(ROC)
truth <- c(rep(1, 200), rep(0, 800))
AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank))


###################################################
### code chunk number 83: 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 84: num_of_fdr_less_than_2
###################################################
sum(result$q.value < 0.2)


###################################################
### code chunk number 85: 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 86: 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 87: str_of_simulation_field
###################################################
str(tcc$simulation)


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


###################################################
### code chunk number 89: 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 90: calc_AUC_of_simulation_data
###################################################
calcAUCValue(tcc)


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


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


###################################################
### code chunk number 93: 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 94: 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 95: 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 96: 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 97: 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 98: analysis_without_replicate_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")
calcAUCValue(tcc)


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


###################################################
### code chunk number 100: 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 101: plot_fc_pseudocolor_multigroup
###################################################
plotFCPseudocolor(tcc)


###################################################
### code chunk number 102: 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 103: 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 104: 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 105: 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 106: 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 107: 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 108: plot_fc_pseudocolor_multifactor
###################################################
head(tcc$count)
tcc$group
plotFCPseudocolor(tcc)


###################################################
### code chunk number 109: 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 110: maplot_simulate_data_1000
###################################################
plot(tcc)


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


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


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


###################################################
### code chunk number 114: 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 115: 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 116: 2390118599
###################################################
sessionInfo()