This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to control for factors like sex in an LRT in DESeq2

Hi all, I'm working on using DESeq2 to analyze a more complicated experiment than I've done before. I have six biological replicates of four conditions, and I want to identify differentially expressed genes between the four conditions.

I did not expect sex to have much of an impact, but when I created a PCA the samples were clearly segregating according to sex, so I would like to control for it.

Here's what my metadata looks like: 

Sample1 ConditionA  Male
Sample2 ConditionA  Female
Sample3 ConditionA  Male
Sample4 ConditionA  Male
Sample5 ConditionA  Female
Sample6 ConditionA  Female
Sample7 ConditionB  Male
Sample8 ConditionB  Female
Sample9 ConditionB  Male
Sample10    ConditionB  Male
Sample11    ConditionB  Female
Sample12    ConditionB  Female
Sample13    ConditionC  Male
Sample14    ConditionC  Female
Sample15    ConditionC  Male
Sample16    ConditionC  Male
Sample17    ConditionC  Female
Sample18    ConditionC  Female
Sample19    ConditionD  Male
Sample20    ConditionD  Female
Sample21    ConditionD  Male
Sample22    ConditionD  Male
Sample23    ConditionD  Female
Sample24    ConditionD  Female

Here's what I've tried in DESeq, but it doesn't seem to do what I'm aiming for. 

design <- formula(~sex + condition)
reduced <- formula(~1)
rsem.in <- DESeqDataSetFromMatrix(countData = ct, colData = sampleMetaData, design = design, tidy = T)
rsem.in <- rsem.in[ rowSums( counts( rsem.in)) >= 1, ]
rsem.de <- DESeq( rsem.in, test = "LRT", reduced = reduced)
rsem.de.res <- results( rsem.de, alpha = 0.05)
summary(rsem.de.res)

Any insight would be greatly appreciated!

rna-seq deseq2
• 309 views • link
written by Ashley S. • Ashley S.

it doesn't seem to do what I'm aiming for.

  1. What are you aiming for?
  2. How is what it's doing not matching your expectations?
• 0 views • link
written by Ram Ram

I'm aiming to find genes that are differentially expressed between conditions, not genes that are differentially expressed between sexes. I'd like to find the effect of the four conditions while controlling for sex. When I run it as is, the top genes with the lowest p-values seem to have those low p-values be driven by differences in sex rather than differences in condition.

• 0 views • link
written by Ashley S. • Ashley S.

If sex is not a factor, why include it in the design formula? Why not just design = ~ condition?

• 0 views • link
written by Ram Ram

Because I can tell from PCAs that the samples are clustering largely due to sex, so I'm inclined to believe that sex is playing a role in gene expression here. I'm concerned that since it seems to be driving so much variation (it is PC2 in a PCA) I need to control for it in my design. Is that incorrect?

• 0 views • link
written by Ashley S. • Ashley S.

Your process is correct. If sex is PC2, then I'm afraid that doesn't allow for cross-sex comparison. You may have to work with ~condition per sex.

I'm not familiar with the LRT test or the usage of the reduced model, so someone with knowledge on that might be able to help you better.

• 0 views • link
written by Ram Ram

1 answer

Can you show the PCA and associated code?

• 0 views • link
written by ATpoint ATpoint

Sure. Thanks!!

My PCA

rownames(sampleMetaData) <- sampleNames
design <- formula(~sex + tissue)
reduced <- formula(~1)
rsem.in <- DESeqDataSetFromMatrix(countData = ct, colData = sampleMetaData, design = design, tidy = T)
rsem.in <- rsem.in[ rowSumscountsrsem.in)) >= 1, ]
rownamesrsem.in) <-  sapplystrsplitcrownamesrsem.in)), split = '_', fixed = TRUE), function(x) (x[1]))
rsem.de <- DESeqrsem.in, test = "LRT", reduced = reduced)
rsem.de.res <- resultsrsem.de, alpha = 0.05)

#PCA
rsem.rlog <- rlog( rsem.de )
plotPCA(rsem.rlog, intgroup=c("sex", "tissue")) + theme_minimal()
• 0 views • link
updated by Ram Ram • written by Ashley S. • Ashley S.

Hey Ashley, this is actually common enough, i.e., a difference based on sex is discovered via PCA. In some cases, the effect may be small and negligible, while, in others [possibly like yours] it can be inferred that sex ought to indeed be included in the design formula.

I see no problem leaving it in the design formula in this case, and any test statistics that you derive for tissue will be adjusting for sex. One thing, I think that, for your 'reduced' model, you may want to be using (and see example HERE):

design <- formula(~ sex + tissue)
... ...
rsem.de <- DESeqrsem.in, test = "LRT", reduced = '~ sex')

If you want to eliminate the effect of sex for your rlog (or variance-stabilised) data for downstream analyses, then you could simply use:

limma::removeBatchEffect(assay(rsem.rlog), colDatarsem.in)[,'sex'])

By the way, I think that, if you do the following, the effect of sex will disappear on that PCA bi-plot:

rsem.rlog <- rlog( rsem.de, blind = FALSE )
plotPCA(rsem.rlog, intgroup=c("sex", "tissue")) + theme_minimal()

[the effect of sex would also disappear after you use limma::removeBatchEffect()).

• 0 views • link
written by Kevin Blighe Kevin Blighe

Thank you so much!! Using that in the reduced model was exactly what I needed, so that's a huge help. Thank you!!

I just had a couple questions about the way to remove these effects from the rlog data. I tried both methods and I think I'm missing something with both--I'm so sorry that these are probably pretty basic issues.

When I tried to plot a PCA using the blind = FALSE and then the same line of plotPCA, I got this PCA.

https://ibb.co/7GDqXWf

It's very similar to the previous version of my PCA, except the cluster placement is switched, so it doesn't seem as though the effect of sex is diminished? Is there any insight you have on this?

Also, I'm not sure I quite understand how to use the limma function of removeBatchEffect, or at least how to use the output in downstream analyses.

When I try this as suggested: 

rsem.BatEffRem <- limma::removeBatchEffect(assay(rsem.log), colDatarsem.in)[,'sex'])

rsem.BatEffRem is a matrix that I then don't know how to handle. Functions like plotPCA don't seem to accept it, and I can't take the log because there are negative values within it. Do you have any guidance about how to continue using this matrix to generate PCAs/heatmaps/etc?

Thank you so much!!

• 0 views • link
updated by Kevin Blighe Kevin Blighe • written by Ashley S. • Ashley S.

Hey again, hmmm, sometimes the blind = FALSE approach works, while not so in others. If the effect of sex is inconsistent, then it will not work. As I think about it, differences in sex will affect gene expression disproportionately; whereas, for example, 'consistent' differences induced by, e.g., single vs. paired end sequencing technologies would affect all genes in a consistent fashion and be easier to eliminate by this soft approach.

Another thing about which to be aware is that the DESeq2 PCA function automatically removes a large chunk of your dataset based on gene-wise variance, i.e., before it even generates the bi-plot.

With the output of, removeBatchEffect(), for PCA, you can probably use my own Bioconductor package: PCAtools: everything Principal Component Analysis

You probably just need to do:

require(PCAtools)
rsem.BatEffRem.pca <- pca(rsem.BatEffRem, metadata = colData(dds))
biplot(rsem.BatEffRem.pca,
  colby = 'sex',
  shape = 'tissue')

You may have to re-set the row and colnames of rsem.BatEffRem, though.

This removes no genes by default. You can remove based on variance, like in DESeq2's function, via the removeVar parameter that is passed to pca()

That matrix should also be ready for clustering analysis. ComplexHeatmap is the supreme clustering / heatmap package, but has learning curve.

• 0 views • link
written by Kevin Blighe Kevin Blighe

Great, thank you!! I will try your package!

• 0 views • link
written by Ashley S. • Ashley S.

Log in to answer this question.

More posts like this