Hello,

I've heard of two different ways of doing batch effect correction:

1. Explicit:

        design(myData_collapsed) <- ~ Sex + Genotype
        dds <- DESeq(myData_collapsed)
        deseqOutput <- results(dds, alpha = 0.05, contrasts = c('Genotype', 'A', 'B'))

2. Implicit:

        # make a full model matrix
        mod  <- model.matrix(~ Genotype, colData(myData_collapsed))
        # make a null model to compare it to
        mod0 <- model.matrix(~   1, colData(myData_collapsed))
        
        # calculate SVs
        svseq <- svaseq( assay(myData_collapsed), mod, mod0, n.sv = 4)
        myData_collapsed$SV1 <- svseq$sv[,1]
        myData_collapsed$SV2 <- svseq$sv[,2]
        myData_collapsed$SV3 <- svseq$sv[,3]
        myData_collapsed$SV4 <- svseq$sv[,4]
    
        # redesign matrix
        design(myData_collapsed) <- ~ SV1 + SV2 + SV3 + SV4 + Genotype
    
         # DESeq
        dds <- DESeq(myData_collapsed)
        deseqOutput <- results(dds, alpha = 0.05, contrasts = c('Genotype', 'A', 'B'))


I only ever see people doing one of the two strategies above. However, it seems to me that it should be possible to do both, just as long as you make sure to include the explicit batches in the design of your model matrix when calculating SVs, like so:


    # make a full model matrix
    mod  <- model.matrix(~ Sex + Genotype, colData(myData_collapsed))
    
    # make a null model to compare it to
    mod0 <- model.matrix(~   1, colData(myData_collapsed))
    
    # make SVs
    svseq <- svaseq( assay(myData_collapsed), mod, mod0, n.sv = 4)
    
    # new DESeq matrix design
    myData_collapsed$SV1 <- svseq$sv[,1]
    myData_collapsed$SV2 <- svseq$sv[,2]
    myData_collapsed$SV3 <- svseq$sv[,3]
    myData_collapsed$SV4 <- svseq$sv[,4]
    design(myData_collapsed) <- ~ SV1 + SV2 + SV3 + SV4 + Sex + Genotype
    
    # perform DESeq
    dds <- DESeq(myData_collapsed)
    deseqOutput <- results(dds, alpha = 0.05, contrasts = c('Genotype', 'A', 'B'))



Is that correct?

(Or - just realized this - should the null model include the explicit batch effect, e.g. mod0 <- model.matrix(~   Sex, colData(myData_collapsed)) ?)

Thank you for your help! I appreciate all the help I've been getting from BioStars recently, and I hope I can do a good job at documentation for others Googling this issue.





The starting point to any RNA-seq experiment should be a well-designed study that negates the need for any form of batch adjustment / correction. This is often not the case, as we know.

In the 'explicit' part, you are not actually doing any adjustment for batch. You've just got gender (sex) and genotype in the design model - are you implying that gender is your 'batch' effect for which you are aiming to adjust? Gender can be a confounding factor in mixed-gender / transgender RNA-seq studies, in particular for genes transcribed from chromosome X and / or those that have escaped chrX inactivation.

In your 'implicit' part, with your code, you are looking for any residual batch (or other) effects in your data without knowing *a priori* what they may be. This type of adjustment should not be performed without good reason. Performing it without justification could lead to unnecessary adjustment of your coefficients and the obtaining of unrealistic *P* values during differential expression. It could be performed in a situation like this: our laboratory receives samples of minimal information and we are told to process them. We realise, after an initial test, that there are likely unknown batch effects, so, we employ SVA to model these effects and adjust for them (as you do).


-----------------------------------------------
--------------------------------

# batch effect known

When we know what is the batch effect, we should include it in the design formula, e.g.:

    design = ~ batch + sex + genotype

This has the effect of modelling the known batch effect and making adjustments for it when it comes to differential expression analysis. No further adjustments from this are necessary unless you wish to use your normalised / transformed counts for other tools, downstream. In such a situation, there are 2 options:

 - use `removeBatchEffect()` (*limma*) on the transformed expression levels to subtract out the batch effect
 - set `blind=FALSE` when transforming by regularised log or variance
   stabilisation (*DESeq2* only) ( see here: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#blind-dispersion-estimation )


# batch effect(s) unknown

If the batch effect(s) is/are unknown, we can use SVA, as you have done, in which case we arrive at:

    design = ~ SV1 + SV2 + Sex + Genotype

*NB - include as minimal a number of surrogate variables in the model as possible; otherwise, you'll over-adjust*

You can then later subtract out the modeled batch effect as above.

----------------------------------------------------

Yet other methods aim to directly adjust raw or normalised counts for either of the situations mentioned above. I go over this here in this previous answer:
https://www.biostars.org/p/266507/#280157

Kevin
