Continuous covariates

[BEFH]

Description of feature

Many users including us would like to be able to include continuous covariates for the differential expression tests.

They could be added to the contrasts sheet like so:

id,variable,reference,target,blocking,continuous
condition_control_treated,condition,control,treated,,age;rin
condition_control_treated_blockrep,condition,control,treated,replicate;batch,

Then in the modules/nf-core/deseq2/differential/templates/deseq_de.R:

At the beginning:

opt <- list(
    output_prefix = ifelse('$task.ext.prefix' == 'null', '$meta.id', '$task.ext.prefix'),
    count_file = '$counts',
    sample_file = '$samplesheet',
    contrast_variable = '$contrast_variable',
    reference_level = '$reference',
    target_level = '$target',
    blocking_variables = NULL,
    continuous_variables = NULL,
    control_genes_file = '$control_genes_file',
    transcript_lengths_file = '$transcript_lengths_file',
    sizefactors_from_controls = FALSE,
    gene_id_col = "gene_id",
    sample_id_col = "experiment_accession",
    subset_to_contrast_samples = FALSE,
    exclude_samples_col = NULL,
    exclude_samples_values = NULL,
    test = "Wald",
    fit_type = "parametric",
    sf_type = 'ratio',
    min_replicates_for_replace = 7,
    use_t = FALSE,
    lfc_threshold = 0,
    alt_hypothesis = 'greaterAbs',
    independent_filtering = TRUE,
    p_adjust_method = 'BH',
    alpha = 0.1,
    minmu = 0.5,
    vs_method = 'vst', # 'rlog', 'vst', or 'rlog,vst'
    shrink_lfc = TRUE,
    cores = 1,
    vs_blind = TRUE,
    vst_nsub = 1000
)

Then for building the variable list:

contrast_variable <- make.names(opt\$contrast_variable)
blocking.vars <- c()

if (!contrast_variable %in% colnames(sample.sheet)) {
    stop(
        paste0(
        'Chosen contrast variable \"',
        contrast_variable,
        '\" not in sample sheet'
        )
    )
} else if (any(!c(opt\$reflevel, opt\$treatlevel) %in% sample.sheet[[contrast_variable]])) {
    stop(
        paste(
        'Please choose reference and treatment levels that are present in the',
        contrast_variable,
        'column of the sample sheet'
        )
    )
} else if (is_valid_string(opt\$blocking_variables)) {
    blocking.vars = make.names(unlist(strsplit(opt\$blocking_variables, split = ';')))
    if (!all(blocking.vars %in% colnames(sample.sheet))) {
        missing_block <- paste(blocking.vars[! blocking.vars %in% colnames(sample.sheet)], collapse = ',')
        stop(
            paste(
                'Blocking variables', missing_block,
                'do not correspond to sample sheet columns.'
            )
        )
    }
}

if (is_valid_string(opt\$continuous_variables)) {
    continuous.vars = make.names(unlist(strsplit(opt\$continuous_variables, split = ';')))
    if (!all(continuous.vars %in% colnames(sample.sheet))) {
        missing_continuous <- paste(continuous.vars[! continuous.vars %in% colnames(sample.sheet)], collapse = ',')
        stop(
            paste(
                'Continuous variables', missing_continuous,
                'do not correspond to sample sheet columns.'
            )
        )
    }
}

And finally for building the model:

# Now specify the model. Use cell-means style so we can be explicit with the
# contrasts

model <- '~ 0'

if (is_valid_string(opt\$blocking_variables)) {
    model <- paste(model, paste(blocking.vars, collapse = ' + '), sep=' + ')
}

if (is_valid_string(opt\$continuous_variables)) {
    model <- paste(model, paste(continuous.vars, collapse = ' + '), sep=' + ')
}

# Make sure all the appropriate variables are factors

for (v in c(blocking.vars, contrast_variable)) {
    sample.sheet[[v]] <- as.factor(sample.sheet[[v]])
}

for (v in continuous.vars) {
    sample.sheet[[v]] <- as.numeric(sample.sheet[[v]])
}

# Variable of interest goes last, see
# https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs

model <- paste(model, contrast_variable, sep = ' + ')

However, there is still clearly work to be done in the Nextflow pipeline and documentation, and I don't know what to change there. If anyone could provide guidance on what to edit and what testing to do, I could do a pull request, but I don't know what to do at this point.

[WackerO]

Hey there @BEFH,
it would indeed be great if you could open some PRs for this!

The first thing would be to change the deseq2 module.

I think if possible, you should also provide a test for the new functionality. The current deseq2 tests are done with the mus_musculus dataset, you can find the links to the files here.
If the samplesheet does not contain any info that is suitable to use as a covariate, you could check the test-datasets repo; maybe there is some dataset that already works, otherwise you can possibly provide something suitable?

[BEFH]

I think there needs to be some discussion over how to do this. My suggestion of an additional contrast sheet column is one way to do it, or we could automatically detect if the column looks like a float (double in R) before turning it into a factor and instead make sure it's coerced to a double, or there could be a pipeline option to list continuous variables (probably preferable).

[WackerO]

Hmm, I'm not sure if adding another param would be the best solution.
As this concerns the contrasts, I think probably your other two suggestions make more sense. The check for numeric values is the easiest solution I think (but must make sure R understands negatives and decimals; yml sometimes struggles with that).
The additional column in the contrasts file is of course more specific, meaning that it would allow people to also have factor variables that are numeric...However, @pinin4fjords, I think for this one would also need to change the validator, no?

[pinin4fjords]

This really just needs to be done in more or less the same way as the blocking. So a column for continuous covariates in the contrasts file. Yes, we'd have to make sure the validator understood about the new column, and verified that the listed covariates were present etc.

There's no reason the existing test data can't be adapted, we can add a new fake continuous covariate column- the test data doesn't need to have scientific validity.

[BEFH]

My one concern with the design here is whether DESeq2 uses a Type I SS or a Type II SS. If it uses a Type II model, then linear covariates and blocking factors can be in different columns with no issues because order doesn't matter. If it's type I, I don't actually know.

[WackerO]

I'm not quite sure what you mean by Type I or II SS, but according to the DESeq2 vignette, the order of variables in the design does not matter (search for order of, second result).