6.1 Two-group comparison

This section compares normpatch and TCC for a two-group comparison. The example uses hypoData, a small count matrix distributed with TCC. The same count matrix and group labels are supplied to both packages.

data(hypoData)

group <- factor(c(1, 1, 1, 2, 2, 2))
exp_design <- data.frame(group = group)

6.1.1 edgeR screening

Both normpatch and TCC support a DEGES workflow that estimates TMM normalization factors and screens candidate DEGs with an edgeR quasi-likelihood test. With matching DEGES settings, the two packages return identical normalization factors in this example.

tcc_tmm <- new("TCC", hypoData, group)
tcc_tmm <- calcNormFactors(
    tcc_tmm,
    norm.method = "tmm",
    test.method = "edger",
    iteration = 3,
    FDR = 0.1,
    floorPDEG = 0.05
)

tcc_tmm$norm.factors

##   G1_rep1   G1_rep2   G1_rep3   G2_rep1   G2_rep2   G2_rep3 
## 0.8706680 0.8501287 0.8389209 1.0844163 1.1491726 1.2066935

x_tmm <- newSeqCountData(hypoData, exp_design = exp_design)
x_tmm <- calc_nf(
    x_tmm,
    norm_func = norm_tmm,
    test_func = test_edger,
    iter = 3,
    p_cutoff = 0,
    q_cutoff = 0.1,
    min_rm_genes = 0.05
)
x_tmm@meta$nf

##   G1_rep1   G1_rep2   G1_rep3   G2_rep1   G2_rep2   G2_rep3 
## 0.8706680 0.8501287 0.8389209 1.0844163 1.1491726 1.2066935

The normalization factors are identical in this example.

tcc_tmm$norm.factors - x_tmm@meta$nf

## G1_rep1 G1_rep2 G1_rep3 G2_rep1 G2_rep2 G2_rep3 
##       0       0       0       0       0       0

6.1.2 DESeq2 screening

Both normpatch and TCC also support a DESeq2-based DEGES workflow, but the default screening test is not identical for the two-group case. normpatch uses RLE normalization with a DESeq2 LRT during DEGES, whereas the original TCC DESeq2 workflow uses RLE normalization with a Wald test. Therefore, the default DESeq2 workflows are expected to give different normalization factors for this example.

tcc_deseq2 <- new("TCC", hypoData, group)
tcc_deseq2 <- calcNormFactors(
    tcc_deseq2,
    norm.method = "deseq2",
    test.method = "deseq2",
    iteration = 3,
    FDR = 0.1,
    floorPDEG = 0.05
)
tcc_deseq2$norm.factors

##   G1_rep1   G1_rep2   G1_rep3   G2_rep1   G2_rep2   G2_rep3 
## 0.8741251 0.8628172 0.8169570 1.0872078 1.1621538 1.1967390

x_deseq2_lrt <- newSeqCountData(hypoData, exp_design = exp_design)
x_deseq2_lrt <- calc_nf(
    x_deseq2_lrt,
    norm_func = norm_rle,
    test_func = test_deseq2,
    iter = 3,
    p_cutoff = 0,
    q_cutoff = 0.1,
    min_rm_genes = 0.05
)
x_deseq2_lrt@meta$nf

##   G1_rep1   G1_rep2   G1_rep3   G2_rep1   G2_rep2   G2_rep3 
## 0.8749343 0.8632637 0.8171796 1.0860253 1.1618237 1.1967733

To reproduce the normalization factors from the original TCC DESeq2 workflow with normpatch, define a custom DEG screening function based on the DESeq2 Wald test. The function follows the calc_nf() testing interface: it receives a count matrix, an experimental design data frame, and normalization factors, and then returns a data frame containing p_value and q_value.

test_deseq2_wald <- function(x, exp_design, nf = NULL, ...) {
    if (is.null(nf)) nf <- rep(1, ncol(x))

    x <- as.matrix(x)
    exp_design <- as.data.frame(exp_design)
    rownames(exp_design) <- colnames(x)

    effective_lib_sizes <- nf * colSums(x)

    dds <- DESeqDataSetFromMatrix(
        countData = x,
        colData = exp_design,
        design = ~ group
    )
    sizeFactors(dds) <- effective_lib_sizes / mean(effective_lib_sizes)
    dds <- estimateDispersions(dds)
    dds <- nbinomWaldTest(dds)

    test_output <- as.data.frame(results(dds))
    p <- test_output$pvalue
    p[is.na(p)] <- 1
    q <- stats::p.adjust(p, method = "BH")

    data.frame(
        p_value = p,
        q_value = q,
        row.names = NULL,
        check.names = FALSE
    )
}

x_deseq2_wald <- newSeqCountData(hypoData, exp_design = exp_design)
x_deseq2_wald <- calc_nf(
    x_deseq2_wald,
    norm_func = norm_rle,
    test_func = test_deseq2_wald,
    iter = 3,
    p_cutoff = 0,
    q_cutoff = 0.1,
    min_rm_genes = 0.05
)
x_deseq2_wald@meta$nf

##   G1_rep1   G1_rep2   G1_rep3   G2_rep1   G2_rep2   G2_rep3 
## 0.8741251 0.8628172 0.8169570 1.0872078 1.1621538 1.1967390

With the custom Wald screening function, the normalization factors estimated by normpatch match those produced by the original TCC DESeq2 workflow in this example.

tcc_deseq2$norm.factors - x_deseq2_wald@meta$nf

## G1_rep1 G1_rep2 G1_rep3 G2_rep1 G2_rep2 G2_rep3 
##       0       0       0       0       0       0