Hello all,

I'm having difficulty in creating a design formula for my rna-seq experiment. I have two main groups each group has its own internal control. How can I first normalize them according to their internal control and then compare the two groups with each other?

My coldata:

How can I design my coldata matrix and do the comparison with deseq2?

Thank you.

A few questions:

• What is the nature of the 2 controls in relation to drug-continued/withdrawed?
• Why is every row duplicated? (replicates?)
• What is the difference between "first" and "second"?

Assuming first/second refer to batches (i.e. there is a technical or biological difference unrelated to treatment):

• Then you have 2 factors (make sure they are actually stored as factors): "group" with 2 levels, and "treatment" with 4 levels. The reference levels in their current order will be "first" and "control1", and will form the intercept. So the design looks like it should simply be "~ group + treatment".

• Important: If "group" can affect "treatment" then it needs to be included as an interaction term, so then your design would be "~ group + treatment + group:treatment". For example: if instead of "group", it was "sex". Males may respond differently to treatment than females. So the design would be "~ group + treatment + sex:treatment". I hope that's clear.

If first/second are not batches (i.e just an arbitrary grouping made by you to keep control2 with drug-withdrawed):

• It might be better to concatenate the names as "first_control1" etc. Then you have one factor (treatment) with 4 levels and the design would just be "~treatment"

This blog: A guide to designs and contrasts in DESeq2 has helped me a lot with my contrasts. Take a look at section 4 - Multiple factors multiple levels. When you get to extracting results is when you will need to take care with "deseq2::results(dds, contrast = ???)". The author also builds a contraster() function that is more intuitive, and helped me realise what I am actually comparing when using "contrast =". There is also an explanation of use of "~ 0" (as mentioned by James above) with another worked example if you choose to go that way.

If you are more mathematically minded ( I am not ) ... this might help: What is a multifactorial design with interactions and negative control experiments.

How can I first normalize them

Don't. Deseq2 will do it for you once you figure out the correct format for your metadata, design and contrasts.

You are asking about an interaction term, which is (drug-continued - control1) - (drug-withdrawed - control2), which assumes that there are actual differences between the two control samples. If the two controls are essentially the same thing, they cancel out and then it's just drug-continued - drug-withdrawed.

There's probably a more automatic way of doing this, but I don't totally get how DESeq2 does the contrasts, so I do it the dumb direct way. As a fake example

> dds <- makeExampleDESeqDataSet()
> design(dds)
~condition
> colData(dds)
DataFrame with 12 rows and 1 column
         condition
          <factor>
sample1          A
sample2          A
sample3          A
sample4          A
sample5          A
...            ...
sample8          B
sample9          B
sample10         B
sample11         B
sample12         B

## change this to have four levels

> colData(dds) <- DataFrame(condition = factor(rep(LETTERS[1:4], each = 3)))
> colData(dds)
DataFrame with 12 rows and 1 column
    condition
     <factor>
1           A
2           A
3           A
4           B
5           B
...       ...
8           C
9           C
10          D
11          D
12          D

## and change to a cell means model, which isn't totally necessary, but maybe easier to understand

> design(dds) <- ~ 0 + condition
> design(dds)
~0 + condition
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

## the signs in the interaction term I described above, after 
## distributing into the parentheses are 1, -1, -1, 1
## so we use that as the contrast

> results(dds, c(1, -1, -1, 1))
log2 fold change (MLE): +1,-1,-1,+1 
Wald test p-value: +1,-1,-1,+1 
DataFrame with 1000 rows and 6 columns
          baseMean log2FoldChange
         <numeric>      <numeric>
gene1     88.17007       0.895912
gene2     11.37401       0.267528
gene3     10.31525      -3.993132
gene4      2.00921       0.069651
gene5     18.40443      -0.322246
...            ...            ...
gene996    13.0700       1.448528
gene997    14.3798       0.200951
gene998    79.2857      -0.550139
gene999    16.5750      -0.014954
gene1000   14.8213      -2.309181
             lfcSE       stat
         <numeric>  <numeric>
gene1     0.592058  1.5132152
gene2     1.463057  0.1828555
gene3     1.442082 -2.7690057
gene4     2.609663  0.0266897
gene5     1.069997 -0.3011655
...            ...        ...
gene996   0.992003  1.4602058
gene997   1.149728  0.1747813
gene998   0.633383 -0.8685729
gene999   1.124772 -0.0132952
gene1000  1.049954 -2.1993176
             pvalue      padj
          <numeric> <numeric>
gene1    0.13022500  0.852226
gene2    0.85491143  0.998859
gene3    0.00562277  0.763831
gene4    0.97870726  0.998859
gene5    0.76328833  0.998859
...             ...       ...
gene996   0.1442335  0.852768
gene997   0.8612515  0.998859
gene998   0.3850808  0.979411
gene999   0.9893923  0.998859
gene1000  0.0278553  0.763831

How can I first normalize them according to their internal control

That's generally not how it's done. You compare the two groups you want, and see what genes are different. If you want to know differences between drug_continued and drug_stopped, compare them to each other. Don't add the complication of variability and uncertainty from a third set of samples.