Hello,

I am facing a tricky problem of experimental design that I would like to share with you. 

We are working with several groups of samples, each of them containing sorted cells of the same cell type. In groups A and B, we have a certain percentage of the cells which have been modified in different ways (=continuous value). Group C is a control group without any modification.

The goal of this experiment is to assess the impact of the modification on the cells by comparing A vs C and B vs C while controlling for the number of modified cells in each group.

In addition, we don't know if this modification has a linear impact on the cells: i.e. a continuous value of 1% could have 10 times less impact than a continuous value of 10% which could itself have 1000 times less impact than a 20% of cells.

The design table:

sample    group    continuous_value(%)

sample_a    A    35
sample_b    A    10
sample_c    B    1
sample_d    B    4
sample_e    C     0
sample_f    C    0

I have already tried two different approaches with DESeq2 to work on this dataset:

a. Use the continuous value as a numeric value in the design (with or without log2 transformation). I obtained some differentially expressed genes but have no idea how to say that it was the right thing to do. 

b. Transform the value into small bins. Unfortunately this did not work ("Error in DESeqDataSet(se, design = design, ignoreRank) : the model matrix is not full rank, so the model cannot be fit as specified.one or more variables or interaction terms in the design formula are linear combinations of the others and must be removed"). As you can see my controls are always at 0, my group B is quite low and the % in the group A is alway higher than the rest. Moreover, we currently don't have enough biological evidences to say: from this % to this one the cells can be fitted in the same box.

Since I don't understand enough all the mechanics behind linear model I would like to have some advices on this design. How can I compare my groups while controlling for this variable number of cells? (assuming or not that the % has a linear impact). 

Up to know I have used DESeq2 but I could also use other packages if they are more suitable for this kind of messy design.

Sorry for this long post,

Thanks in advance for your answers!

Hmm. The main problem with your set-up is that your ABC groupings are correlated with the continuous variable. In particular, group A has the highest variable values, while group B has lower values and group C has values of zero. Now, if you observe an increasing trend in expression from A to B to C, there's no way to know whether that's due to genuine DE between groups or whether it's being driven by the continuous variable (especially if you don't know how expression responds to changes in the continuous variable, i.e., the shape of the trend).

As a result, just putting in the continuous variable as a coefficient in the design matrix will probably be liberal. Significant deviations from the fitted line will be considered to represent DE between groups, but they may well be due to an uninteresting non-linear response in expression with respect to the continuous variable. If you want to be more conservative, you could do something like this:

group <- rep(c("A", "B", "C"), each=2)
cont.var <- c(35, 10, 4, 1, 0, 0)
X <- splines::ns(cont.var, df=2)
design <- model.matrix(~0 + group + X)

This fits a cubic spline to the expression values against changes in the continuous variables, which is a bit more flexible than the linear trend you would get from using cont.var directly in design. As a result of this flexibility, both the full and null models (when comparing groups) will fit equally well. This will reduce significance and avoid detecting spurious differences between groups that are driven by cont.var, at the cost of losing power to detect genuine DE.

Note that the choice of df=2 is simply because there's not enough points for more spline coefficients. Increasing df would saturate the model and preclude dispersion estimation. Ryan's suggestion is relevant here; if you have more samples, you should show them, because a more flexible setting of df may be appropriate.

You could also log-transform cont.var prior to using ns; this might provide better performance by making the spacing between points more similar, such that the outlier at 35 doesn't dominate the fit. Make sure you add a prior to avoid undefined values at zero, though.

Thanks a lot to all of you! I will give a try to what you advised.

Hello!

I finished to test your advices. Unfortunately the MDS plots I have seen are really not great and I could not detect any gene differentially expressed after the applying some contrasts. Could you tell me if the code provided makes sens to you?

Moreover, I was a bit surprise by the results of the spline function who added not one but five columns to the design table (df=5). Could you tell me if this is expected? 

Code:

library(limma)
library(edgeR)
# Load the data (column information in colData and matrix of count in matrix)
colData <- read.table(..., sep="\t", header=T, row.names=1)
Matrix <- read.table(..., sep="\t", header=T, row.names=1)
Matrix <- DGEList(Matrix)

# Design matrix
group <- factor(colData$group)
cont.var <- splines::ns(as.vector(colData$cont.var), df=5)
design <- model.matrix(~0+group+cont.var)
colnames(design) <- c("A", "B", "C", "D", "df1", "df2", "df3", "df4", "df5")

design[1:5,]

# Normalization of the reads with TMM
​Matrix <- calcNormFactors(Matrix, method="TMM")

# Remove the entries with not enough information
isexpr <- rowSums(cpm(Matrix) > 10) >= 2
Matrix <- Matrix[isexpr,]

# Voom 
Matrix <- voom(Matrix, design, plot=TRUE)
par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plotMDS(Matrix, xlim=c(-2.5,2.5), col= c(1,3,2,4,3,3,2,2,2,1,4,4,1,2,1), pch=c(16,18,17,19,18,18,17,17,17,16,19,19,16,17,16), cex=1.2, xlab="Leading logFC dim 1", ylab= "Leading logFC dim 2")
legend("topright", legend=c("D", "C", "A", "B"), pch=c(16,17,18,19), cex=1.2, title= "Groups", col= c(1,2,3,4), inset=c(-0.25,-0,2,0))

# Fit the GLM model, the contrasts and get the results
fit <- lmFit(Matrix,design)
contrast.matrix <- makeContrasts(A-B, A-C, ... , levels = design)
fit <- contrasts.fit(fit, contrast.matrix)
fit.de <- eBayes(fit)
topTable(fit.de, adjust="BH",number=50, coef=1) 

Thanks in advance