<p>It looks good to me. One think I could suggest is to see the distribution of your data/model (with limma/voom).</p>
<p>After you get the list of your genes with:</p>
<pre><code>top.table &lt;- topTable(fit.eBayes, sort.by = &quot;P&quot;, n = Inf)
</code></pre>
<p>Maybe you can check your multidimensional scaling (MDS):</p>
<pre><code>plotMDS(your_data, col = as.numeric(group))
</code></pre>
<p>Then you can also check with voom</p>
<pre><code>y &lt;- voom(data, model.matrix, plot = T)
</code></pre>
<p>I hope this helps</p>

CC BY-SA 4.0
2021-09-23T21:22:46.660Z
Thank you for your answer. Checking the distro might be very useful, indeed!
Lukas
CC BY-SA 4.0
2021-09-23T11:38:28.033Z
2021-09-23T11:38:28.033Z
flavinsky
<p>I'm starting to use Limma to find differential expressed sites within my data and would like to ask if my approach was correct.</p>
<p>The data has a control and a treatment set, each containing data of 48 patients, which were measured before treatment (t0) and after treatment (t1). So I have Control-t0 and Control-t1 as well as Treatment-t0 and Treatment-t1, whereas Control-t0 and Treatment-t0 should not show any effect.
Now I would like to know if there is a significant change in Treatment compared to Control. Therefore, I think, I have to adjust for the baseline (t0), right?</p>
<p>The code below is probably the most verbose option Limma provides, but I wanted to use the long method first before I use the shortcuts.
So would the Interaction-coefficient I defined within my cont-matrix be the difference between Control and Treatment adjusted for the baseline?
Also, I had to turn of adjusting for multiple testing here as my example data is random</p>
<pre><code>

#Data.frame with columns representing samples and rows different sites/genes
data_df &lt;- data.frame(replicate(8, sample(100,1000, rep=TRUE)))
colnames(data_df) &lt;- c(&quot;Patient_01_t0&quot;, &quot;Patient_01_t1&quot;,
                       &quot;Patient_02_t0&quot;, &quot;Patient_02_t1&quot;,
                       &quot;Patient_03_t0&quot;, &quot;Patient_03_t1&quot;,
                       &quot;Patient_04_t0&quot;, &quot;Patient_04_t1&quot;)
                     
#Metadata for data_df
meta &lt;- data.frame(Sample= colnames(data_df),
                   ID= c(rep(&quot;Patient_01&quot;,2),
                         rep(&quot;Patient_02&quot;,2),
                         rep(&quot;Patient_03&quot;,2),
                         rep(&quot;Patient_04&quot;,2)),
                   Time= c(rep(c(&quot;t0&quot;,&quot;t1&quot;),4)),
                   Group=c(rep(&quot;Control&quot;,4),
                           rep(&quot;Treatment&quot;,4))
)
                         
#Create levels for model-matrix
time &lt;- factor(meta<span class="math-container">$Time, labels=c('t0','t1') )
group&lt;- factor(meta$</span>Group, labels=c('Control','Treatment') )
lev=paste(group,&quot;.&quot;,time,sep=&quot;&quot;)
lev


design &lt;- model.matrix(~0+lev)
design

#Create contrast matrix with all factors of interest
cont_matrix &lt;- makeContrasts( t1_VS_t0_Ctrl = levControl.t1 - levControl.t0,
                              t1_VS_t0_Treat = levTreatment.t1 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t0= levControl.t0 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t1= levControl.t1 - levTreatment.t1,
                              Interaction = (levControl.t1-levControl.t0)-(levTreatment.t1-levTreatment.t0), 
                              levels=design)

#Calculate different expressed sites
fit &lt;- lmFit(data_df, design) 
fit &lt;- contrasts.fit(fit, cont_matrix) 
fit.eBayes &lt;- eBayes(fit, robust=TRUE)

summary(decideTests(fit.eBayes, adjust.method=&quot;none&quot;))
</code></pre>

<r><differential-expression><limma>
Help with Limma-model
# Help with Limma-model

(Asked by: username_0 on 2021-09-22)

I'm starting to use Limma to find differential expressed sites within my data and would like to ask if my approach was correct.
The data has a control and a treatment set, each containing data of 48 patients, which were measured before treatment (t0) and after treatment (t1). So I have Control-t0 and Control-t1 as well as Treatment-t0 and Treatment-t1, whereas Control-t0 and Treatment-t0 should not show any effect.
Now I would like to know if there is a significant change in Treatment compared to Control. Therefore, I think, I have to adjust for the baseline (t0), right?
The code below is probably the most verbose option Limma provides, but I wanted to use the long method first before I use the shortcuts.
So would the Interaction-coefficient I defined within my cont-matrix be the difference between Control and Treatment adjusted for the baseline?
Also, I had to turn of adjusting for multiple testing here as my example data is random


#Data.frame with columns representing samples and rows different sites/genes
data_df <- data.frame(replicate(8, sample(100,1000, rep=TRUE)))
colnames(data_df) <- c("Patient_01_t0", "Patient_01_t1",
                       "Patient_02_t0", "Patient_02_t1",
                       "Patient_03_t0", "Patient_03_t1",
                       "Patient_04_t0", "Patient_04_t1")
                     
#Metadata for data_df
meta <- data.frame(Sample= colnames(data_df),
                   ID= c(rep("Patient_01",2),
                         rep("Patient_02",2),
                         rep("Patient_03",2),
                         rep("Patient_04",2)),
                   Time= c(rep(c("t0","t1"),4)),
                   Group=c(rep("Control",4),
                           rep("Treatment",4))
)
                         
#Create levels for model-matrix
time <- factor(meta$Time, labels=c('t0','t1') )
group<- factor(meta$Group, labels=c('Control','Treatment') )
lev=paste(group,".",time,sep="")
lev


design <- model.matrix(~0+lev)
design

#Create contrast matrix with all factors of interest
cont_matrix <- makeContrasts( t1_VS_t0_Ctrl = levControl.t1 - levControl.t0,
                              t1_VS_t0_Treat = levTreatment.t1 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t0= levControl.t0 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t1= levControl.t1 - levTreatment.t1,
                              Interaction = (levControl.t1-levControl.t0)-(levTreatment.t1-levTreatment.t0), 
                              levels=design)

#Calculate different expressed sites
fit <- lmFit(data_df, design) 
fit <- contrasts.fit(fit, cont_matrix) 
fit.eBayes <- eBayes(fit, robust=TRUE)

summary(decideTests(fit.eBayes, adjust.method="none"))

**Tags:** r, differential-expression, limma

---

## 1 answers

### Answer by username_1 on 2021-09-23. Score: 3 | Accepted answer

It looks good to me. One think I could suggest is to see the distribution of your data/model (with limma/voom).
After you get the list of your genes with:
top.table <- topTable(fit.eBayes, sort.by = "P", n = Inf)

Maybe you can check your multidimensional scaling (MDS):
plotMDS(your_data, col = as.numeric(group))

Then you can also check with voom
y <- voom(data, model.matrix, plot = T)

I hope this helps

> **Comments:**
> **username_0**: Thank you for your answer. Checking the distro might be very useful, indeed!  

---

