tests/testthat/test_05cbcbseq.R

In kokrah/cbcbSEQ: RNAseq analysis pipeline for UMD CBCB collaborators with batch effect assessment and correction

library(cbcbSEQ)
library(testthat)
##library(hpgltools)
context("Test hpgltools and cbcbSEQ")
library(pasilla)
data(pasillaGenes)

## WTF why does travis give me this:
## Error in loadNamespace(name) : there is no package called 'cbcbSEQ'
## The damn library call is just above!

## Load the pasilla data set
message("Loading pasilla, setting up count tables.")
datafile = system.file("extdata/pasilla_gene_counts.tsv", package="pasilla")
## Load the counts and drop super-low counts genes
counts = read.table(datafile, header=TRUE, row.names=1)
counts = counts[rowSums(counts) > ncol(counts),]
## Set up a quick design to be used by cbcbSEQ and hpgltools
design = data.frame(row.names=colnames(counts),
    condition=c("untreated","untreated","untreated",
        "untreated","treated","treated","treated"),
    libType=c("single-end","single-end","paired-end",
        "paired-end","single-end","paired-end","paired-end"))
metadata = design
colnames(metadata) = c("condition", "batch")
metadata$sampleid = rownames(metadata)

## Make sure it is still possible to create an expt
message("Setting up an expt class to contain the pasilla data and metadata.")
pasilla_expt = create_expt(count_dataframe=counts, meta_dataframe=metadata)
cbcb_data = as.matrix(counts)
##pasilla_expt = create_expt(count_dataframe=counts, meta_dataframe=metadata)
##hpgl_data = Biobase::exprs(pasilla_expt$expressionset)
##test_that("Does data from an expt equal a raw dataframe?", {
##    expect_equal(cbcb_data, hpgl_data)
##})

## Check that normalization tools work similarly
message("Testing quantile raw normalization.")
cbcb_quantile = cbcbSEQ::qNorm(cbcb_data)
##hpgl_quantile_data = hpgl_norm(pasilla_expt, transform="raw", norm="quant", convert="raw", filter_low=FALSE, verbose=TRUE)
##hpgl_quantile = hpgl_quantile_data[["count_table"]]
##test_that("Are the quantile normalizations identical?", {
##    expect_equal(cbcb_quantile, hpgl_quantile)
##})

message("Testing quantile(cpm()) normalization using edgeR's cpm().")
cbcb_qcpm = cbcbSEQ::qNorm(cbcb_data)
##library(edgeR)
## I don't know how to call cpm without using a library call first.
##cbcb_qcpm = edgeR::cpm(cbcb_qcpm)

cbcb_quantile = cbcbSEQ::qNorm(cbcb_data)
##hpgl_quantile = hpgl_norm(pasilla_expt, norm="quant")
##hpgl_quantile = hpgl_quantile$count_table
##test_that("Are quantiles identical?", {
##    expect_equal(cbcb_quantile, hpgl_quantile)
##})

##hpgl_qcpm = hpgl_norm(pasilla_expt, norm="quant", convert="edgecpm", filter_low=FALSE, verbose=TRUE)
##hpgl_qcpm = hpgl_qcpm$final_counts$count_table
##hpgl_qcpm = hpgl_qcpm$count_table
##test_that("Are cpm conversions identical?", {
##    expect_equal(cbcb_qcpm, hpgl_qcpm)
##})

## log2/cpm that
message("Testing log2(quantile(cpm())) normalization using the cpm from voom()")
cbcb_l2qcpm_data = cbcbSEQ::log2CPM(cbcb_quantile)
cbcb_l2qcpm = cbcb_l2qcpm_data$y
##hpgl_l2qcpm_data = hpgl_norm(pasilla_expt, transform="log2", norm="quant", convert="cpm", filter_low=FALSE)
##hpgl_l2qcpm = hpgl_l2qcpm_data[["count_table"]]
##hpgl_l2qcpm_expt = normalize_expt(pasilla_expt, transform="log2", norm="quant", convert="cpm", filter_low=FALSE)
##hpgl_l2qcpm2 = Biobase::exprs(hpgl_l2qcpm_expt$expressionset)
##test_that("Are l2qcpm conversions/transformations identical using two codepaths?", {
##    expect_equal(cbcb_l2qcpm, hpgl_l2qcpm)
##    expect_equal(cbcb_l2qcpm, hpgl_l2qcpm2)
##})

## Check that PCA invocations are similar
message("Testing PCA invocations.")
cbcb_svd = cbcbSEQ::makeSVD(cbcb_l2qcpm)
##hpgl_pca_info = hpgl_pca(hpgl_l2qcpm_expt)
##hpgl_svd = hpgl_pca_info$pca
cbcb_res = cbcbSEQ::pcRes(cbcb_svd$v, cbcb_svd$d, design$condition, design$libType)
##hpgl_res = hpgl_pca_info$res
##test_that("Do the PCA invocations provide the same results?", {
##    expect_equal(cbcb_svd$v, hpgl_svd$v)
##    expect_equal(cbcb_svd$d, hpgl_svd$d)
##    expect_equal(cbcb_res, hpgl_res)
##})

## invocation of batch correction
message("Testing batch correction results using modified combat.")
## sva no longer has a function 'design.mat'
##cbcb_batch = cbcbSEQ::combatMod(cbcb_l2qcpm, batch=design$libType, mod=design$condition, noScale=TRUE)
##hpgl_batch_expt = normalize_expt(pasilla_expt, transform="log2", norm="quant", convert="cpm", batch="combatmod", filter_low=FALSE)
##hpgl_batch = exprs(hpgl_batch_expt$expressionset)
##test_that("Does combat batch correction end in the same dataframe?", {
##    expect_equal(cbcb_batch, hpgl_batch)
##})

## voom invocation
## This is a peculiar thing, the cbcbSEQ log2CPM() returns the normalized libsizes rather than those following log2() transformation.
## It then goes on to call voom() with this libsize rather than the log2().
## As a result, my normalization function is now keeping a copy of the libsizes and count tables from beginning to end
## to ensure that this is still accessible when required.\
message("Testing that the libsize input to voom is identical.")
cbcb_libsize = cbcb_l2qcpm_data$lib.size
##hpgl_libsize = hpgl_l2qcpm_data[["intermediate_counts"]][["normalization"]][["libsize"]]
##test_that("In preparing for voom(), are the library sizes the same?", {
##    expect_equal(cbcb_libsize, hpgl_libsize)
##})

## Given that, lets try a voom() invocation and see what happens
message("Testing different methods of invoking voom, with/without cpm/log2 converted/transformed data.")
condition = design$condition
test_model = model.matrix(~condition)
cbcb_voom = cbcbSEQ::voomMod(x=as.matrix(cbcb_l2qcpm), design=test_model, lib.size=cbcb_libsize)
##hpgl_voom = cbcbSEQ::voomMod(x=as.matrix(hpgl_l2qcpm), design=test_model, lib.size=hpgl_libsize, plot=TRUE)
##hpgl_voom2 = hpgltools::hpgl_voom(as.matrix(hpgl_l2qcpm), model=test_model, libsize=hpgl_libsize, logged=TRUE, converted=TRUE)
##hpgl_voom3 = hpgltools::hpgl_voom(as.matrix(hpgl_quantile), test_model, libsize=hpgl_libsize, logged=FALSE, converted=FALSE)

##message("I can't test cbcb_voom vs. hpgl_voom3 because I set the row/colnames.")
##test_that("Do different voom() invocations end with the same data?", {
##    expect_equal(cbcb_voom, hpgl_voom)
##    expect_equal(cbcb_voom$E, hpgl_voom2$E)
##    expect_equal(cbcb_voom$E, hpgl_voom3$E)
##})

## my hpgl_voom() sets row/column names and causes a test of the weights to fail.
## But checking manually shows them the same.
## expect_equal(cbcb_voom$weights, hpgl_voom$weights)

message("Test the final results from limma by invoking it through different codepaths.")
cbcb_fit = limma::lmFit(cbcb_voom)
cbcb_eb = limma::eBayes(cbcb_fit)
cbcb_top = limma::topTable(cbcb_eb, number=nrow(cbcb_eb))
##hpgl_toptables = hpgltools::limma_pairwise(hpgl_l2qcpm_expt, model_intercept=TRUE, model_batch=FALSE, libsize=hpgl_libsize)
##hpgl_top = hpgl_toptables$all_tables
##test_that("Limma results.", {
##    expect_equal(cbcb_top, hpgl_top)
##})