Answers.Body: <p>I think the problem is in the design. There is no room for error or variation from a common ground.</p>
  
  <p>I'm not sure of the solution, but I think you'll need to delete a column to give the design a degree of freedom, and probably add the intercept. You'll need to adjust the contrasts accordingly. </p>
  
  <p>However, I would recommend to ask this at <a href="https://support.bioconductor.org" rel="nofollow noreferrer">support.bioconductor.org</a>. There are more experts on linear modelling and contrasts there than me. (If you ask there make it easier to copy and paste your code)</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-17T14:35:36.097Z
Answers.Text: Thanks for the replies. I'm not really sure how to do the modifications you recommended to my design or contrasts. I am writing a question on bioconductor, but how can I make it easier to copy and paste my code? This was my first post on a bioinformatics support forum.
Answers.UserDisplayName: Eddie
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-18T09:29:46.047Z
Answers.Text: Ideally you could use a package as [reprex](https://cran.r-project.org/package=reprex). But if you use code blocks and remove the `>` at the beginning of the code, comment the output (so that if I paste it on the console it doesn't error) you'll do good. Another thing to consider is that I cannot `load("KC_Raw.RData")` because I don't have it, so try to make the example self contained. But the code provided is enough to understand your problem so don't worry if you use the same code used here there.
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-16T21:24:03.403Z
Answers.LastActivityDate: 2020-06-16T21:24:03.403Z
Answers.UserDisplayName: llrs
Body: <p>I am working with a RNA-seq data set in maize that has a relatively complex design. There are two levels of treatment A (nitrogen fertilizer level in the field, high or low), two levels of treatment B (nitrogen nutrients in <em>in vitro</em> cultures, high and low) and two levels of treatment C (two time points of sampling), all with 3 reps. </p>
  
  <pre><code>&gt; library(edgeR)
  &gt; load("KC_Raw.RData")
  &gt; y &lt;- DGEList(counts = KCraw.data[,2:25])
  &gt; keep &lt;- rowSums(cpm(y) &gt; 10) &gt;= 3
  &gt; targets &lt;- data.frame(rownames=colnames(KCraw.data)[2:25] ,
  +                       Time=rep(c(rep("2DIC",12),rep("5DIC",12))) ,
  +                       FieldN=rep(c(rep("FH",6), rep("FL",6)),2) ,
  +                       CultureN=rep(c(rep("CL",3),rep("CH",3)),4))
  &gt; Group &lt;- factor(paste(targets$FieldN,targets$Time,targets$CultureN,sep="."))
  &gt; targets &lt;- cbind(targets,Group=Group)
  &gt; targets
     rownames Time FieldN CultureN      Group
  1   KC1_H2L 2DIC     FH       CL FH.2DIC.CL
  2   KC2_H2L 2DIC     FH       CL FH.2DIC.CL
  3   KC3_H2L 2DIC     FH       CL FH.2DIC.CL
  4   KC4_H2H 2DIC     FH       CH FH.2DIC.CH
  5   KC5_H2H 2DIC     FH       CH FH.2DIC.CH
  6   KC6_H2H 2DIC     FH       CH FH.2DIC.CH
  7   KC7_L2L 2DIC     FL       CL FL.2DIC.CL
  8   KC8_L2L 2DIC     FL       CL FL.2DIC.CL
  9   KC9_L2L 2DIC     FL       CL FL.2DIC.CL
  10 KC10_L2H 2DIC     FL       CH FL.2DIC.CH
  11 KC11_L2H 2DIC     FL       CH FL.2DIC.CH
  12 KC12_L2H 2DIC     FL       CH FL.2DIC.CH
  13 KC13_H5L 5DIC     FH       CL FH.5DIC.CL
  14 KC14_H5L 5DIC     FH       CL FH.5DIC.CL
  15 KC15_H5L 5DIC     FH       CL FH.5DIC.CL
  16 KC16_H5H 5DIC     FH       CH FH.5DIC.CH
  17 KC17_H5H 5DIC     FH       CH FH.5DIC.CH
  18 KC18_H5H 5DIC     FH       CH FH.5DIC.CH
  19 KC19_L5L 5DIC     FL       CL FL.5DIC.CL
  20 KC20_L5L 5DIC     FL       CL FL.5DIC.CL
  21 KC21_L5L 5DIC     FL       CL FL.5DIC.CL
  22 KC22_L5H 5DIC     FL       CH FL.5DIC.CH
  23 KC23_L5H 5DIC     FL       CH FL.5DIC.CH
  24 KC24_L5H 5DIC     FL       CH FL.5DIC.CH
  </code></pre>
  
  <p>I have used edgeR in R to calculate differential expression for contrasts involving 3 reps at one treatment combination to 3 reps at another treatment combination, for example</p>
  
  <pre><code>&gt; y &lt;- DGEList(counts = KCraw.data[keep,2:25], group = Group)
  &gt; y &lt;- calcNormFactors(y)
  &gt; 
  &gt; TMM &lt;- KCraw.data[keep,2:25]
  &gt; for (i in 1:24) {
  +   TMM[,i] &lt;- TMM[,i] / (y<span class="math-container">$samples$</span>lib.size[i] * y<span class="math-container">$samples$</span>norm.factors[i]) * 1e6
  + }
  &gt; 
  &gt; y &lt;- DGEList(counts = TMM,group = Group)
  &gt; 
  &gt; design &lt;- model.matrix(~0+Group)
  &gt; colnames(design) &lt;- levels(Group)
  &gt; y &lt;- calcNormFactors(y,method = "TMM")
  &gt; y &lt;- estimateDisp(y,design)
  &gt; fitQL &lt;- glmQLFit(y,design)
  &gt; fit &lt;- glmFit(y,design)
  &gt; myKC.contrasts &lt;- makeContrasts(
  +   H2H.H2L = FH.2DIC.CH - FH.2DIC.CL,
  +   L2H.L2L = FL.2DIC.CH - FL.2DIC.CL,
  +   H2H.L2H = FH.2DIC.CH - FL.2DIC.CH,
  +   H2L.L2L = FH.2DIC.CL - FL.2DIC.CL,
  +   H5H.H5L = FH.5DIC.CH - FH.5DIC.CL,
  +   L5H.L5L = FL.5DIC.CH - FL.5DIC.CL,
  +   H5H.L5H = FH.5DIC.CH - FL.5DIC.CH,
  +   H5L.L5L = FH.5DIC.CL - FL.5DIC.CL,
  +   H2H.L2L = FH.2DIC.CH - FL.2DIC.CL,
  +   H5H.L5L = FH.5DIC.CH - FL.5DIC.CL,
  +   H5L.H2L = FH.5DIC.CL - FH.2DIC.CL, 
  +   H5H.H2H = FH.5DIC.CH - FH.2DIC.CH,
  +   L5L.L2L = FL.5DIC.CL - FL.2DIC.CL,
  +   L5H.L2H = FL.5DIC.CH - FL.2DIC.CH,
  +   levels=design)
  &gt; design
     FH.2DIC.CH FH.2DIC.CL FH.5DIC.CH FH.5DIC.CL FL.2DIC.CH FL.2DIC.CL FL.5DIC.CH FL.5DIC.CL
  1           0          1          0          0          0          0          0          0
  2           0          1          0          0          0          0          0          0
  3           0          1          0          0          0          0          0          0
  4           1          0          0          0          0          0          0          0
  5           1          0          0          0          0          0          0          0
  6           1          0          0          0          0          0          0          0
  7           0          0          0          0          0          1          0          0
  8           0          0          0          0          0          1          0          0
  9           0          0          0          0          0          1          0          0
  10          0          0          0          0          1          0          0          0
  11          0          0          0          0          1          0          0          0
  12          0          0          0          0          1          0          0          0
  13          0          0          0          1          0          0          0          0
  14          0          0          0          1          0          0          0          0
  15          0          0          0          1          0          0          0          0
  16          0          0          1          0          0          0          0          0
  17          0          0          1          0          0          0          0          0
  18          0          0          1          0          0          0          0          0
  19          0          0          0          0          0          0          0          1
  20          0          0          0          0          0          0          0          1
  21          0          0          0          0          0          0          0          1
  22          0          0          0          0          0          0          1          0
  23          0          0          0          0          0          0          1          0
  24          0          0          0          0          0          0          1          0
  attr(,"assign")
  [1] 1 1 1 1 1 1 1 1
  attr(,"contrasts")
  attr(,"contrasts")$Group
  [1] "contr.treatment"
  
  &gt; myKC.contrasts
              Contrasts
  Levels       H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
    FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
    FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
    FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
    FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
    FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
    FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
    FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
    FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
              Contrasts
  Levels       L5H.L2H
    FH.2DIC.CH       0
    FH.2DIC.CL       0
    FH.5DIC.CH       0
    FH.5DIC.CL       0
    FL.2DIC.CH      -1
    FL.2DIC.CL       0
    FL.5DIC.CH       1
    FL.5DIC.CL       0
  </code></pre>
  
  <p>After analyzing these contrasts, I wanted to estimate some sort of simple effect, such as the culture media nitrogen level. To do this, I ran the following code.</p>
  
  <pre><code>&gt; myKC.contrasts &lt;- cbind(myKC.contrasts,
  +                         Development = c(1,1,-1,-1,1,1,-1,-1),
  +                         FieldN = c(1,1,1,1,-1,-1,-1,-1),
  +                         CultureN = c(1,-1,1,-1,1,-1,1,-1)
  + )
  &gt; myKC.contrasts
             H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
  FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
  FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
  FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
  FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
  FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
  FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
  FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
  FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
             L5H.L2H Development FieldN CultureN
  FH.2DIC.CH       0           1      1        1
  FH.2DIC.CL       0           1      1       -1
  FH.5DIC.CH       0          -1      1        1
  FH.5DIC.CL       0          -1      1       -1
  FL.2DIC.CH      -1           1     -1        1
  FL.2DIC.CL       0           1     -1       -1
  FL.5DIC.CH       1          -1     -1        1
  FL.5DIC.CL       0          -1     -1       -1
  </code></pre>
  
  <p>Once I rerun the analysis for the CultureN contrast and look at the result for a particular gene, I see that it's estimated log2FC is equal to the sum of every simple contrast. </p>
  
  <pre><code>&gt; lrt &lt;- glmQLFTest(fitQL, contrast=myKC.contrasts[,"CultureN"])
  &gt; topTags(lrt,n=nrow(y$counts))["GRMZM2G445575",]
  Coefficient:  1*FH.2DIC.CH -1*FH.2DIC.CL 1*FH.5DIC.CH -1*FH.5DIC.CL 1*FL.2DIC.CH -1*FL.2DIC.CL 1*FL.5DIC.CH -1*FL.5DIC.CL 
                   logFC   logCPM        F       PValue          FDR
  GRMZM2G445575 -6.63617 5.417106 151.5261 3.691525e-11 2.825777e-08
  # FC is a data frame of the logFC of each constrast in columns for each gene in rows
  &gt; sum(FC["GRMZM2G445575",c("H2H.H2L","L2H.L2L","H5H.H5L","L5H.L5L")])
  [1] -6.636197
  </code></pre>
  
  <p>My first question is if this analysis is a valid way of summarizing the simple effects of each treatment.
  I would like to be able to also include the effects of the H2H.L2L and H5H.L5L contrast in the FieldN and CultureN comparison, but I am not sure how to do this, or if this would be valid because each of these contrasts includes treatments that have different levels of two treatment factors. </p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-06-16T21:26:32.630Z
Comments.Text: I don't quite get why you need to manually change the counts with a for loop. It is used as a size factor in the linear model. Now, for the pairwise comparison, what you have in the first part is ok
Comments.UserDisplayName: StupidWolf
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-06-16T21:27:57.930Z
Comments.Text: For the general culture N etc, run a linear model and run the estimateGLM... etc stuff, it should give you the effects you want. You cannot get everything from one model
Comments.UserDisplayName: StupidWolf
Tags: <rna-seq><differential-expression><edger><design>
Title: Comparing multiple treatments to multiple other treatments in edgeR for simple effects in a complex experimental design
ThreadText: # Comparing multiple treatments to multiple other treatments in edgeR for simple effects in a complex experimental design
  
  (Asked by: username_0 on 2020-06-16)
  
  I am working with a RNA-seq data set in maize that has a relatively complex design. There are two levels of treatment A (nitrogen fertilizer level in the field, high or low), two levels of treatment B (nitrogen nutrients in *in vitro* cultures, high and low) and two levels of treatment C (two time points of sampling), all with 3 reps.
  
  ```
  > library(edgeR)
  > load("KC_Raw.RData")
  > y <- DGEList(counts = KCraw.data[,2:25])
  > keep <- rowSums(cpm(y) > 10) >= 3
  > targets <- data.frame(rownames=colnames(KCraw.data)[2:25] ,
  +                       Time=rep(c(rep("2DIC",12),rep("5DIC",12))) ,
  +                       FieldN=rep(c(rep("FH",6), rep("FL",6)),2) ,
  +                       CultureN=rep(c(rep("CL",3),rep("CH",3)),4))
  > Group <- factor(paste(targets$FieldN,targets$Time,targets$CultureN,sep="."))
  > targets <- cbind(targets,Group=Group)
  > targets
     rownames Time FieldN CultureN      Group
  1   KC1_H2L 2DIC     FH       CL FH.2DIC.CL
  2   KC2_H2L 2DIC     FH       CL FH.2DIC.CL
  3   KC3_H2L 2DIC     FH       CL FH.2DIC.CL
  4   KC4_H2H 2DIC     FH       CH FH.2DIC.CH
  5   KC5_H2H 2DIC     FH       CH FH.2DIC.CH
  6   KC6_H2H 2DIC     FH       CH FH.2DIC.CH
  7   KC7_L2L 2DIC     FL       CL FL.2DIC.CL
  8   KC8_L2L 2DIC     FL       CL FL.2DIC.CL
  9   KC9_L2L 2DIC     FL       CL FL.2DIC.CL
  10 KC10_L2H 2DIC     FL       CH FL.2DIC.CH
  11 KC11_L2H 2DIC     FL       CH FL.2DIC.CH
  12 KC12_L2H 2DIC     FL       CH FL.2DIC.CH
  13 KC13_H5L 5DIC     FH       CL FH.5DIC.CL
  14 KC14_H5L 5DIC     FH       CL FH.5DIC.CL
  15 KC15_H5L 5DIC     FH       CL FH.5DIC.CL
  16 KC16_H5H 5DIC     FH       CH FH.5DIC.CH
  17 KC17_H5H 5DIC     FH       CH FH.5DIC.CH
  18 KC18_H5H 5DIC     FH       CH FH.5DIC.CH
  19 KC19_L5L 5DIC     FL       CL FL.5DIC.CL
  20 KC20_L5L 5DIC     FL       CL FL.5DIC.CL
  21 KC21_L5L 5DIC     FL       CL FL.5DIC.CL
  22 KC22_L5H 5DIC     FL       CH FL.5DIC.CH
  23 KC23_L5H 5DIC     FL       CH FL.5DIC.CH
  24 KC24_L5H 5DIC     FL       CH FL.5DIC.CH
  
  ```
  
  I have used edgeR in R to calculate differential expression for contrasts involving 3 reps at one treatment combination to 3 reps at another treatment combination, for example
  
  ```
  > y <- DGEList(counts = KCraw.data[keep,2:25], group = Group)
  > y <- calcNormFactors(y)
  > 
  > TMM <- KCraw.data[keep,2:25]
  > for (i in 1:24) {
  +   TMM[,i] <- TMM[,i] / (y$samples$lib.size[i] * y$samples$norm.factors[i]) * 1e6
  + }
  > 
  > y <- DGEList(counts = TMM,group = Group)
  > 
  > design <- model.matrix(~0+Group)
  > colnames(design) <- levels(Group)
  > y <- calcNormFactors(y,method = "TMM")
  > y <- estimateDisp(y,design)
  > fitQL <- glmQLFit(y,design)
  > fit <- glmFit(y,design)
  > myKC.contrasts <- makeContrasts(
  +   H2H.H2L = FH.2DIC.CH - FH.2DIC.CL,
  +   L2H.L2L = FL.2DIC.CH - FL.2DIC.CL,
  +   H2H.L2H = FH.2DIC.CH - FL.2DIC.CH,
  +   H2L.L2L = FH.2DIC.CL - FL.2DIC.CL,
  +   H5H.H5L = FH.5DIC.CH - FH.5DIC.CL,
  +   L5H.L5L = FL.5DIC.CH - FL.5DIC.CL,
  +   H5H.L5H = FH.5DIC.CH - FL.5DIC.CH,
  +   H5L.L5L = FH.5DIC.CL - FL.5DIC.CL,
  +   H2H.L2L = FH.2DIC.CH - FL.2DIC.CL,
  +   H5H.L5L = FH.5DIC.CH - FL.5DIC.CL,
  +   H5L.H2L = FH.5DIC.CL - FH.2DIC.CL, 
  +   H5H.H2H = FH.5DIC.CH - FH.2DIC.CH,
  +   L5L.L2L = FL.5DIC.CL - FL.2DIC.CL,
  +   L5H.L2H = FL.5DIC.CH - FL.2DIC.CH,
  +   levels=design)
  > design
     FH.2DIC.CH FH.2DIC.CL FH.5DIC.CH FH.5DIC.CL FL.2DIC.CH FL.2DIC.CL FL.5DIC.CH FL.5DIC.CL
  1           0          1          0          0          0          0          0          0
  2           0          1          0          0          0          0          0          0
  3           0          1          0          0          0          0          0          0
  4           1          0          0          0          0          0          0          0
  5           1          0          0          0          0          0          0          0
  6           1          0          0          0          0          0          0          0
  7           0          0          0          0          0          1          0          0
  8           0          0          0          0          0          1          0          0
  9           0          0          0          0          0          1          0          0
  10          0          0          0          0          1          0          0          0
  11          0          0          0          0          1          0          0          0
  12          0          0          0          0          1          0          0          0
  13          0          0          0          1          0          0          0          0
  14          0          0          0          1          0          0          0          0
  15          0          0          0          1          0          0          0          0
  16          0          0          1          0          0          0          0          0
  17          0          0          1          0          0          0          0          0
  18          0          0          1          0          0          0          0          0
  19          0          0          0          0          0          0          0          1
  20          0          0          0          0          0          0          0          1
  21          0          0          0          0          0          0          0          1
  22          0          0          0          0          0          0          1          0
  23          0          0          0          0          0          0          1          0
  24          0          0          0          0          0          0          1          0
  attr(,"assign")
  [1] 1 1 1 1 1 1 1 1
  attr(,"contrasts")
  attr(,"contrasts")$Group
  [1] "contr.treatment"
  
  > myKC.contrasts
              Contrasts
  Levels       H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
    FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
    FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
    FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
    FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
    FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
    FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
    FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
    FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
              Contrasts
  Levels       L5H.L2H
    FH.2DIC.CH       0
    FH.2DIC.CL       0
    FH.5DIC.CH       0
    FH.5DIC.CL       0
    FL.2DIC.CH      -1
    FL.2DIC.CL       0
    FL.5DIC.CH       1
    FL.5DIC.CL       0
  
  ```
  
  After analyzing these contrasts, I wanted to estimate some sort of simple effect, such as the culture media nitrogen level. To do this, I ran the following code.
  
  ```
  > myKC.contrasts <- cbind(myKC.contrasts,
  +                         Development = c(1,1,-1,-1,1,1,-1,-1),
  +                         FieldN = c(1,1,1,1,-1,-1,-1,-1),
  +                         CultureN = c(1,-1,1,-1,1,-1,1,-1)
  + )
  > myKC.contrasts
             H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
  FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
  FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
  FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
  FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
  FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
  FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
  FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
  FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
             L5H.L2H Development FieldN CultureN
  FH.2DIC.CH       0           1      1        1
  FH.2DIC.CL       0           1      1       -1
  FH.5DIC.CH       0          -1      1        1
  FH.5DIC.CL       0          -1      1       -1
  FL.2DIC.CH      -1           1     -1        1
  FL.2DIC.CL       0           1     -1       -1
  FL.5DIC.CH       1          -1     -1        1
  FL.5DIC.CL       0          -1     -1       -1
  
  ```
  
  Once I rerun the analysis for the CultureN contrast and look at the result for a particular gene, I see that it's estimated log2FC is equal to the sum of every simple contrast.
  
  ```
  > lrt <- glmQLFTest(fitQL, contrast=myKC.contrasts[,"CultureN"])
  > topTags(lrt,n=nrow(y$counts))["GRMZM2G445575",]
  Coefficient:  1*FH.2DIC.CH -1*FH.2DIC.CL 1*FH.5DIC.CH -1*FH.5DIC.CL 1*FL.2DIC.CH -1*FL.2DIC.CL 1*FL.5DIC.CH -1*FL.5DIC.CL 
                   logFC   logCPM        F       PValue          FDR
  GRMZM2G445575 -6.63617 5.417106 151.5261 3.691525e-11 2.825777e-08
  # FC is a data frame of the logFC of each constrast in columns for each gene in rows
  > sum(FC["GRMZM2G445575",c("H2H.H2L","L2H.L2L","H5H.H5L","L5H.L5L")])
  [1] -6.636197
  
  ```
  
  My first question is if this analysis is a valid way of summarizing the simple effects of each treatment.
  I would like to be able to also include the effects of the H2H.L2L and H5H.L5L contrast in the FieldN and CultureN comparison, but I am not sure how to do this, or if this would be valid because each of these contrasts includes treatments that have different levels of two treatment factors.
  
  **Tags:** rna-seq, differential-expression, edger, design
  
  > **Comments:**
  > **username_1**: I don't quite get why you need to manually change the counts with a for loop. It is used as a size factor in the linear model. Now, for the pairwise comparison, what you have in the first part is ok  
  > **username_1**: For the general culture N etc, run a linear model and run the estimateGLM... etc stuff, it should give you the effects you want. You cannot get everything from one model  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2020-06-16. Score: 0
  
  I think the problem is in the design. There is no room for error or variation from a common ground.
  
  I'm not sure of the solution, but I think you'll need to delete a column to give the design a degree of freedom, and probably add the intercept. You'll need to adjust the contrasts accordingly.
  
  However, I would recommend to ask this at [support.bioconductor.org](https://support.bioconductor.org). There are more experts on linear modelling and contrasts there than me. (If you ask there make it easier to copy and paste your code)
  
  > **Comments:**
  > **username_0**: Thanks for the replies. I'm not really sure how to do the modifications you recommended to my design or contrasts. I am writing a question on bioconductor, but how can I make it easier to copy and paste my code? This was my first post on a bioinformatics support forum.  
  > **username_2**: Ideally you could use a package as [reprex](https://cran.r-project.org/package=reprex). But if you use code blocks and remove the `>` at the beginning of the code, comment the output (so that if I paste it on the console it doesn't error) you'll do good. Another thing to consider is that I cannot `load("KC_Raw.RData")` because I don't have it, so try to make the example self contained. But the code provided is enough to understand your problem so don't worry if you use the same code used here there.  
  
  ---
  
  