### 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 = "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 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: 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 16: 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 pval <- result$pval qval <- result$padj 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 ### 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 17: 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 18: 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 19: 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 20: 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 21: 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 22: 0408278414 ################################################### 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 = 1, FDR = 0.1, floorPDEG = 0.05) tcc$norm.factors ################################################### ### code chunk number 23: 4389572193 ################################################### 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 pval <- result$pval qval <- result$padj 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 ### 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 24: 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 25: 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 26: 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 27: 7573490899 ################################################### 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 = "deseq", test.method = "deseq", iteration = 1, FDR = 0.1, floorPDEG = 0.05, paired = TRUE) tcc$norm.factors ################################################### ### code chunk number 28: 3489438904 ################################################### 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)) FDR <- 0.1 floorPDEG <- 0.05 cl <- data.frame(group = group, pair = pair) full <- count ~ group + pair reduced <- count ~ pair ### STEP 1 ### cds <- newCountDataSet(hypoData, cl) ### STEP 2 ### cds <- estimateSizeFactors(cds) cds <- estimateDispersions(cds, method = "blind", sharingMode = "fit-only") full.model <- fitNbinomGLMs(cds, full) reduced.model <- fitNbinomGLMs(cds, reduced) pval <- nbinomGLMTest(full.model, reduced.model) pval[is.na(pval)] <- 1 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 ### 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 29: 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 30: 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 31: 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 32: 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 33: 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 34: 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 35: 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) full <- count ~ condition reduced <- count ~ 1 ################################################### ### code chunk number 36: degesDeseq_multigroup_tcc_nf ################################################### tcc <- calcNormFactors(tcc, norm.method = "deseq", test.method = "deseq", iteration = 1, full = full, reduced = reduced) tcc$norm.factors ################################################### ### code chunk number 37: 3894410011 ################################################### 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 = "deseq", test.method = "deseq", iteration = 1) tcc$norm.factors ################################################### ### code chunk number 38: 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) full <- count ~ condition reduced <- count ~ 1 cds <- newCountDataSet(hypoData_mg, group) ### STEP 1 ### cds <- estimateSizeFactors(cds) ### STEP 2 ### cds <- estimateDispersions(cds) full.model <- fitNbinomGLMs(cds, full) reduced.model <- fitNbinomGLMs(cds, reduced) pval <- nbinomGLMTest(full.model, reduced.model) pval[is.na(pval)] <- 1 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 ### 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 39: hypoData_deg_label ################################################### library(TCC) data(hypoData) nonDEG <- 201:1000 summary(hypoData[nonDEG, ]) ################################################### ### code chunk number 40: calc_median_of_nondeg_hypoData ################################################### apply(hypoData[nonDEG, ], 2, median) ################################################### ### code chunk number 41: 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 42: get_normalized_data ################################################### normalized.count <- getNormalizedData(tcc) ################################################### ### code chunk number 43: 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 44: 4390011292 ################################################### buff_1 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 45: 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 46: 4330011292 ################################################### buff_2 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 47: 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) range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 48: 4330110292 ################################################### buff_1 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 49: 3259021231 ################################################### library(TCC) data(hypoData) nonDEG <- 201:1000 group <- c(1, 1, 1, 2, 2, 2) cds <- newCountDataSet(hypoData, group) cds <- estimateSizeFactors(cds) normalized.count <- counts(cds, normalized = TRUE) apply(normalized.count[nonDEG, ], 2, median) range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 50: 1330110292 ################################################### buff_2 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 51: 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 52: 1664110292 ################################################### buff_1 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 53: 238967389 ################################################### 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 54: 1664110292 ################################################### buff_2 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 55: 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 56: 1674110292 ################################################### buff_1 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 57: 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 58: 1674110992 ################################################### buff_2 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 59: 5904580011 ################################################### 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 = "deseq", test.method = "deseq", iteration = 1, FDR = 0.1, floorPDEG = 0.05, paired = TRUE) normalized.count <- getNormalizedData(tcc) apply(normalized.count[nonDEG, ], 2, median) range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 60: 5490200292 ################################################### buff_1 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 61: 4390055561 ################################################### 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)) pair <- factor(c(1, 2, 3, 1, 2, 3)) cl <- data.frame(group = group, pair = pair) cds <- newCountDataSet(hypoData, cl) cds <- estimateSizeFactors(cds) normalized.count <- counts(cds, normalized = TRUE) apply(normalized.count[nonDEG, ], 2, median) range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 62: 9600320992 ################################################### buff_2 <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 63: hypoData_mg_deg_label ################################################### library(TCC) data(hypoData_mg) nonDEG <- 201:1000 summary(hypoData_mg[nonDEG, ]) ################################################### ### code chunk number 64: calc_median_of_nondeg_hypoData_mg ################################################### apply(hypoData_mg[nonDEG, ], 2, median) ################################################### ### code chunk number 65: NormHypoDataMultiGGet ################################################### hypoData_mg.median <- apply(hypoData_mg[nonDEG, ], 2, median) ################################################### ### code chunk number 66: 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 67: degesEdger_tcc_multigroup_normed_data_hdie ################################################### normByiDEGES <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 68: 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 69: tmm_tcc_multigroup_normed_data_hide ################################################### normByTMM <- range(apply(normalized.count[nonDEG, ], 2, median)) ################################################### ### code chunk number 70: 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 71: get_estimated_results ################################################### result <- getResult(tcc, sort = TRUE) head(result) ################################################### ### code chunk number 72: summary_of_estimated_deg ################################################### table(tcc$estimatedDEG) ################################################### ### code chunk number 73: maplot_of_estimated_result ################################################### plot(tcc) ################################################### ### code chunk number 74: 3498757592 ################################################### library(ROC) truth <- c(rep(1, 200), rep(0, 800)) AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank)) ################################################### ### code chunk number 75: 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 76: 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 77: 5408289011 ################################################### 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 = 3, FDR = 0.1, floorPDEG = 0.05) tcc <- estimateDE(tcc, test.method = "deseq2", FDR = 0.1) ################################################### ### code chunk number 78: 5029042481 ################################################### result <- getResult(tcc, sort = TRUE) head(result) ################################################### ### code chunk number 79: 7512913498 ################################################### table(tcc$estimatedDEG) ################################################### ### code chunk number 80: 7512013498 ################################################### tb <- table(tcc$estimatedDEG) ################################################### ### code chunk number 81: 3489400103 ################################################### plot(tcc) ################################################### ### code chunk number 82: 4012399910 ################################################### library(ROC) truth <- c(rep(1, 200), rep(0, 800)) AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank)) ################################################### ### code chunk number 83: 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 84: 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 85: 4090911011 ################################################### buff_1 <- AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank)) ################################################### ### code chunk number 86: 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 87: 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 88: 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 89: 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 90: 4390087612 ################################################### result <- getResult(tcc, sort = TRUE) head(result) ################################################### ### code chunk number 91: 0238735032 ################################################### table(tcc$estimatedDEG) ################################################### ### code chunk number 92: 0233731032 ################################################### tb <- table(tcc$estimatedDEG) ################################################### ### code chunk number 93: 3489320103 ################################################### plot(tcc) ################################################### ### code chunk number 94: 5702784221 ################################################### library(ROC) truth <- c(rep(1, 200), rep(0, 800)) AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank)) ################################################### ### code chunk number 95: 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 96: 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 97: 7400039021 ################################################### 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) ################################################### ### code chunk number 98: 4890184102 ################################################### result <- getResult(tcc, sort = TRUE) head(result) ################################################### ### code chunk number 99: 4901822901 ################################################### table(tcc$estimatedDEG) ################################################### ### code chunk number 100: 4911822001 ################################################### tb <- table(tcc$estimatedDEG) ################################################### ### code chunk number 101: 3387482103 ################################################### plot(tcc) ################################################### ### code chunk number 102: 4919187780 ################################################### library(ROC) truth <- c(rep(1, 200), rep(0, 800)) AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank)) ################################################### ### code chunk number 103: 4929100782 ################################################### library(TCC) data(hypoData) group <- factor(c(1, 2)) cl <- data.frame(group = group) dds <- newCountDataSet(countData = hypoData[, c(1, 4)], condition = group) dds <- estimateSizeFactors(dds) sizeFactors(dds) <- sizeFactors(dds) / mean(sizeFactors(dds)) dds <- estimateDispersions(dds, method = "blind", sharingMode = "fit-only") dds <- nbinomTest(dds, 1, 2) pvalue <- dds$pval library(ROC) truth <- c(rep(1, 200), rep(0, 800)) pvalue[is.na(pvalue)] <- 1 AUC(rocdemo.sca(truth = truth, data = -rank(pvalue))) ################################################### ### code chunk number 104: 8787372620 ################################################### library(TCC) data(hypoData) group <- c(1, 2) tcc <- new("TCC", hypoData[, c(1, 4)], group) tcc <- calcNormFactors(tcc, norm.method = "deseq", iteration = 0) tcc <- estimateDE(tcc, test.method = "deseq", 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 105: 0285012872 ################################################### data(hypoData) colnames(hypoData) <- c("T_dogA", "T_dogB", "T_dogC", "N_dogA", "N_dogB", "N_dogC") head(hypoData) ################################################### ### code chunk number 106: 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 107: 4390911012 ################################################### result <- getResult(tcc, sort = TRUE) head(result) ################################################### ### code chunk number 108: 3909884931 ################################################### table(tcc$estimatedDEG) ################################################### ### code chunk number 109: 3934884932 ################################################### tb <- table(tcc$estimatedDEG) ################################################### ### code chunk number 110: 3482340103 ################################################### plot(tcc) ################################################### ### code chunk number 111: 3490930192 ################################################### library(ROC) truth <- c(rep(1, 200), rep(0, 800)) AUC(rocdemo.sca(truth = truth, data = -tcc$stat$rank)) ################################################### ### code chunk number 112: 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 113: 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 114: 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 115: num_of_fdr_less_than_2 ################################################### sum(result$q.value < 0.2) ################################################### ### code chunk number 116: 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(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 117: 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 118: 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 <- as.data.frame(topTags(lrt, n = nrow(tcc$count))) p.value <- tmp$table$PValue q.value <- tmp$table$FDR result <- cbind(p.value, q.value) head(result) sum(q.value < 0.1) sum(q.value < 0.2) ################################################### ### code chunk number 119: 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 120: 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 ### full <- count ~ condition redueced <- count ~ 1 tcc <- estimateDE(tcc, test.method = "deseq", FDR = 0.1, full = full, reduced = reduced) result <- getResult(tcc, sort = TRUE) head(result) table(tcc$estimatedDEG) ################################################### ### code chunk number 121: 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 ### full <- count ~ condition reduced <- count ~ 1 cds <- newCountDataSet(tcc$count, group) cds <- estimateSizeFactors(cds) sizeFactors(cds) <- sizeFactors(cds) / mean(sizeFactors(cds)) cds <- estimateDispersions(cds) full.model <- fitNbinomGLMs(cds, full) reduced.model <- fitNbinomGLMs(cds, reduced) 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 122: 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 123: str_of_simulation_field ################################################### str(tcc$simulation) ################################################### ### code chunk number 124: tale_of_simulation_trueDEG ################################################### table(tcc$simulation$trueDEG) ################################################### ### code chunk number 125: 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 126: calc_AUC_of_simulation_data ################################################### calcAUCValue(tcc) ################################################### ### code chunk number 127: calc_AUC_of_simulation_data_hide ################################################### auc.degesedger <- calcAUCValue(tcc) ################################################### ### code chunk number 128: calc_AUC_of_simulation_data_ROC ################################################### AUC(rocdemo.sca(truth = as.numeric(tcc$simulation$trueDEG != 0), data = -tcc$stat$rank)) ################################################### ### code chunk number 129: 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 130: 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 131: 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 132: 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 133: 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 134: 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 135: calc_AUC_of_deseq_tcc ################################################### tcc <- calcNormFactors(tcc, norm.method = "deseq", iteration = 0) tcc <- estimateDE(tcc, test.method = "deseq") calcAUCValue(tcc) ################################################### ### code chunk number 136: calc_AUC_of_degesDeseq_org ################################################### cds <- newCountDataSet(tcc$count, tcc$group$group) cds <- estimateSizeFactors(cds) sizeFactors(cds) <- sizeFactors(cds) / mean(sizeFactors(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 137: 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 138: 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 139: plot_fc_pseudocolor_multigroup ################################################### plotFCPseudocolor(tcc) ################################################### ### code chunk number 140: 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 141: 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 142: 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 143: 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 144: 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 145: 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 146: plot_fc_pseudocolor_multifactor ################################################### head(tcc$count) tcc$group plotFCPseudocolor(tcc) ################################################### ### code chunk number 147: 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 148: maplot_simulate_data_1000 ################################################### plot(tcc) ################################################### ### code chunk number 149: norm_simulate_data_1000 ################################################### normalized.count <- getNormalizedData(tcc) colSums(normalized.count) colSums(tcc$count) mean(colSums(tcc$count)) ################################################### ### code chunk number 150: 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 151: maplot_normed_simulate_data_1000 ################################################### plot(tcc, median.lines = TRUE) ################################################### ### code chunk number 152: 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 153: 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 154: 2390118599 ################################################### sessionInfo()