Hi All,

I am writing you because I need help in designing the design matrix for analyzing the results of my RNA-seq experiment.

I have sequenced 57 samples. Each sample is described by 6 factors ( Timepoint, Sex, Treatment, Family, sequencing Lane and Run of extraction).

I am struggling to find an effective way to write a matrix design for such a high number of factor. My first idea was to create a design matrix as follows:

f <- paste(data$TimePoint, data$Sex, data$Treatment, data$Sibgroup, data$Run, data$Lane, sep=".")
f <- as.factor(f)
design <- model.matrix(~f)
colnames(design) <- levels(f)

I wanted to create it this way because it would have made simple to make the contrast matrix (for ex. comparing male vs females would have required just to make the sum of the combination of factor including "male" and substract it to the sum for the combination of factors including the "female" term). 

This doesn't work, because the observed combination of the 6 factors are 57, which means no sample share the exact same combination of factors and voom needs replication of at least one combination to run. 

Can you suggest me a way to write the design and make the constrast matrix (for ex. to compare male vs females).

Thank you in advance,

let me know if you need more informations,

best regards

Oliver

First, check on a MDS plot and see if your samples are segregating according to your factors. I can't imagine the sequencing lane and extraction run are biologically interesting parameters, so if they're not causing major differences in the MDS plot, then you might as well ignore them.

Of the remaining factors, you can stitch them together in an additive model:

design <- model.matrix(~ Sex + Treatment + Sibgroup + ..., data)

Whether this is sensible or not depends on how the factors interact with each other. The pure additive model above assumes that every factor behaves independently; however, you can imagine, e.g., that the response to treatment might vary between sexes. In that case, you might prefer to do something like:

combined <- paste0(data$Sex, ".", data$Treatment)
design <- model.matrix(~ combined + Sibgroup + ..., data)

... in order to account for the interaction between the two factors. The same argument applies with any combination of factors, so you'll need to make some decisions about what to put in combined and what to leave as additive terms.

Also, if TimePoint is a continuous covariate, that will need to be handled differently. Just putting TimePoint as an additive factor assumes that it has a linear effect on gene (log-)expression for all samples. If you do something like:

design <- model.matrix(~ Sex + Sex:TimePoint + ..., data)

... you can get coefficients for sex-specific effects of time (i.e., the effect of time for each sex is different). You can also do this with splines to account for non-linear responses; check out the ns function and examples in the limma user's guide.

P.S. You should tag limma and voom separately, otherwise people who follow either tag won't see this question.

Hi Aaron, thanks for your answer.

Sequencing Lane and extraction Run seem to explain part of the variance on the MDS plot, this is why I want to include them in the model. Time point is a discrete factor with 3 possible values.

Ok so let's say I am writing the model as follows:

combined <- paste0(data$TimePoint, ".", data$Sex, ".", data$Treatment)
design <- model.matrix(~ combined + data$Sibgroup + data$Lane + data$Run)

The colnames of the design matrix are these:

> colnames(design)
 [1] "tp1.Female.control" "tp1.Female.EE2"     "tp1.Male.control"   "tp1.Male.EE2"       "tp2.Female.control" "tp2.Female.EE2"    
 [7] "tp2.Male.control"   "tp2.Male.EE2"       "tp3.Female.control" "tp3.Female.EE2"     "tp3.Male.control"   "tp3.Male.EE2"      
[13] "data$Sibgroup94"    "data$Sibgroup95"    "data$Sibgroup96"    "data$Sibgroup97"    "data$Runr2"         "data$Runr3"        
[19] "data$LaneB"         "data$LaneC"         "data$LaneD"         "data$LaneE"         "data$LaneF"         "data$LaneG"        
[25] "data$LaneH"         "data$LaneI"         "data$LaneL" 

Now, the combined factors have 12 possible combinations, the sibgroup has 5 levels, the Lane 10 and the Run 3. As you can see the design matrix doesn't explicitly show the first levels for the 3 "single" factors (sibgroup, lane and run). My question now is how should write the contrast for, for example: tp1.Female.control vs tp1. Male.control?

Thank you in advance,

Oliver

Ok thanks for your answer. Now let's say that I want to include also the sibgroup inside the combined group. I write the design as:

combined <- paste0(data$TimePoint, ".", data$Sex, ".", data$Treatment, ".", data$Sibgroup)
design <- model.matrix(~ combined + data$Lane + data$Run)

Which returns me this error while I am running the voom command:

Coefficients not estimable: Runr2 Runr3 LaneC LaneD LaneE LaneF LaneH LaneI
Error in approxfun(l, rule = 2) :
  need at least two non-NA values to interpolate
Inoltre: Warning message:
Partial NA coefficients for 35348 probe(s)
Called from: top level 

I don't see why for these factors (for which there is replication on the values that they can assume) it is impossible to calculate the coefficients.

Another thing: let's say I don't want to check any interaction and I write a model as:

design <- model.matrix(~0+TimePoint+Sex+Treatment+Run+Lane, data)

The design matrix will look like:

> colnames(design)
 [1] "TimePointtp1" "TimePointtp2" "TimePointtp3" "SexMale"      "TreatmentEE2" "Runr2"        "Runr3"        "LaneB"       
 [9] "LaneC"        "LaneD"        "LaneE"        "LaneF"        "LaneG"        "LaneH"        "LaneI"        "LaneL"   

In this case, how should I write the contrast to compare male vs female?