Hi,

I’m trying to do some differential expression analysis which I believe that it should be trivial but I couldn’t find any reference for what I’m looking for. I’m pretty sure that the solution will be much shorter than the description below for what I’m looking for (sorry for that…):

I would like to find differentially expressed genes that indicate that a treatment A has a higher impact comparing a treatment B, where each of the treatments is compared to DIFFERENT reference levels as follows:

• I have overall of 24 samples.
• 6 first samples are of 6 different patients (group A) before a treatment ”A” was given. I consider the 6 patients as biological replicants, as I’m not interested to differentiate between patients.
• Next 6 samples are of the same 6 patients as above, but after the treatment “A” was given. Once again, I consider all these 6 samples as biological replicants.
• Next 6 samples are of another group (group B) of 6 patients, before the treatment “B” was given. All of these 6 samples are considered as biological replicants.
• Last 6 samples are of the 6 patients of group B, after the treatment “B” was given.
•  

I would like to find a differential expression for the following: [fold change of  ” after treatment A” vs ”before treatment A”], compared to [fold change of ”after treatment B” vs ” before treatment B”].

I managed easily to apply DESeq2 for a differential expression of  ”before treatment X” vs ”after treatment X”, but now I have 2 result matrices that include LOG2FC, pvalue and adjusted pvalue, and I’m not sure if and how I should proceed and compare between them. I prefer that DESeq2 will do the entire required analysis from the beginning, if possible.

I was looking in the DESeq2’s vignette and I tried the following code which works well, but I have no idea if it does what I mean:

d <- expand.grid(condition=c(rep("beforeT",6), rep("afterT",6)), treatment=c("A", "B"))

coldataBOTH <- data.frame(row.names=colnames(BOTH), d)

ddsBOTH <- DESeqDataSetFromMatrix(countData=BOTH, colData=coldataBOTH, design=~condition+treatment)

#(the matrix BOTH is the raw count matrix that includes 24 libraries)

mcols(ddsBOTH) <- geneID

ddsBOTH <- DESeq(ddsBOTH)

resB <- results(ddsBOTH)

I suspect that what I really need is described under either the section “Multi-factor designs”or probably under ”Interactions” in the following great, just recently updated, vignette:

https://bioconductor.org/packages/3.7/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions

However, I’m not sure that I could find there a good correspondence with what I’m trying to do.

I would really appreciate if someone can show me an explicit example, tailored for my case (I guess that it should be just an addition of 2-3 lines to my code)…

If DESEQ2 is not the right tool I would be happy to get any recommendation for a proper one. Or any other statistical suggestions.

p.s – just for the case that the someone in the future will encounter my question and it’s not what he will be looking for, I found some hints for what I am trying to do in previous topics in this forum:

DESEQ2--compare two log2-transformed fold changes

DESeq2: Paired Test

But once again I failed to relate and implement the help that was given in those topics into my case.

Hi Marina

Mike must have overlooked that you have two different sets of patients, and each patient only gets either treatment "A" or treatment "B". An interaction would be needed if you had the both treatment applied to the same patient, but in your case, it is actually straightforward.

I use the following slightly obscure R code to produce an example sample table:

sampleTable <- data.frame(
   sample = sprintf( "Sample%02d", 1:24 ),
   patient = rep( sprintf( "Pat%02d", 1:12 ), each=2 ),
   treatment = c( rep( c( "before", "afterA" ), 6 ), rep( c( "before", "afterB" ), 6 ) ) )

It produces the following data frame, with one row per sample:

     sample patient treatm
1  Sample01   Pat01 before
2  Sample02   Pat01 afterA
3  Sample03   Pat02 before
4  Sample04   Pat02 afterA
5  Sample05   Pat03 before
6  Sample06   Pat03 afterA
7  Sample07   Pat04 before
8  Sample08   Pat04 afterA
9  Sample09   Pat05 before
10 Sample10   Pat05 afterA
11 Sample11   Pat06 before
12 Sample12   Pat06 afterA
13 Sample13   Pat07 before
14 Sample14   Pat07 afterB
15 Sample15   Pat08 before
16 Sample16   Pat08 afterB
17 Sample17   Pat09 before
18 Sample18   Pat09 afterB
19 Sample19   Pat10 before
20 Sample20   Pat10 afterB
21 Sample21   Pat11 before
22 Sample22   Pat11 afterB
23 Sample23   Pat12 before
24 Sample24   Pat12 afterB

I assume that your sample table looks like this, too. Make your own sample table, with one row per sample, and two columns, labelled patient with the ID of the patient and treatment, with values before, afterA, and afterB.

It is helpful (though not really necessary) to use

mf$treatment <- relevel( mf$treatment, "before" )

to make the before level the default base of comparisons. (Note that we have only one "before" treatment level. After all, before treatment, there is no difference yet between A and B.)

Now, I get me some (random, non-sensical) example count data for 24 samples with

library( DESeq2 )
exampleCounts <- counts(makeExampleDESeqDataSet( m=24 ))

and run DESeq

dds <- DESeqDataSetFromMatrix( exampleCounts, mf, ~ patient + treatment )
dds <- DESeq(dds)

You can now use three result commands to get results for three comparisons of interest:

results( dds, c( "treatment", "afterA", "before" ) )
results( dds, c( "treatment", "afterB", "before" ) )
results( dds, c( "treatment", "afterB", "afterA" ) )

The first two show the genes affected by the two treatments, respectively. The last one compares the two treatments against each other and looks for significant differences. As "before" is the same for both, we do not need to subtract it or even mention it.

Note that our design formula above (~ patient + treatment) mentions patient to account for this so-called "blocking factor". This causes the exact opposite of looking at the differences between patients: namely, it causes the baseline expression values (the "before" expression) to be subtracted so that the treatment effects are calculated as the averages of the differences between after and before treatment, with the difference taken for each patient first, and averaged over paatients afterwards.

If you want to understand this last point better, remind yourself first from an elementary statistics textbook about the difference between an unpaired and a paired t test, and the read up on how this generalizes to blocking factors in ANOVA. What you have is often called a "split-plot design" in textbooks.