I would be very grateful for advice regarding the analysis of my microarray data. I am pretty new to gene expression analysis, and limma seems to be a very powerful package, but I am still a bit overwhelmed by it because it offers so many different functions that I can't make sense of yet.

In my experiment, I thawed a vial of cells and cultured them. The cells were then split into four plates and further cultured. Three of the four plates were treated with different drugs for several days; the fourth plate was not treated and served as a control.

At different time points, I collected a portion of the cells for gene expression analysis.

In analyzing the microarray data, I am now interested in looking at the differentially expressed genes (DEGs) at the different time points to determine whether the prolonged treatment has an enhancing effect, i.e., whether the number of significant DEGs increases. Since the cell culture's duration will also affect the genome, I would like to normalize the treated cells to the control sample of each time point so that I can only see the effects of the treatment.

Unfortunately, we could not make replicates of each sample due to cost, so I only measured one sample per treatment and time point.

How would one proceed here? So far, I have been inspired by chapter 17.4 (Time course effects of corn oil on rat thymus) of the limma Handbook. Here, however, the data are not normalized to a control, so I cannot adopt this pipeline one-to-one. Furthermore, replicates were measured here, and I wonder if I have to use any special functions to interpret my single measurements correctly.

My idea was to analyze the three treatments separately, i.e., first load only the data from the first treatment and the control, normalize the data to the control and then perform linear regression modeling using a y ~ 0 + time function. When the data are normalized to the control, no genes should be differentially expressed at time t0, which is why a forced intercept at (0,0) makes biological sense. From the model, the effect of treatment over time would then be seen. I would subsequently run this pipeline for the remaining two treatments as well.

I am grateful to the community for any hints, code snippets, or advice.

Targets.tsv

FileName    Name        Treatment       Time    Cy3

file1       ctrl_t0     none            0       control
file2       treat1_t0   treat1          0       sample
file3       treat2_t0   treat2          0       sample
file4       treat3_t0   treat3          0       sample
file5       ctrl_t1     none            1       control
file6       treat1_t1   treat1          1       sample
file7       treat2_t1   treat2          1       sample
file8       treat3_t1   treat3          1       sample
file9       ctrl_t2     none            2       control
file10      treat1_t2   treat1          2       sample
file11      treat2_t2   treat2          2       sample
file12      treat3_t2   treat3          2       sample
file13      ctrl_t3     none            3       control
file14      treat1_t3   treat1          3       sample
file15      treat2_t3   treat2          3       sample
file16      treat3_t3   treat3          3       sample
file17      ctrl_t4     none            4       control
file18      treat1_t4   treat1          4       sample
file19      treat2_t4   treat2          4       sample
file20      treat3_t4   treat3          4       sample

Selection of the files of one treatment*

targetinfo <- readTargets(targetsFile, row.names = 'Name')
targetinfo <- subset(targetinfo, Treatment == "treat1" | Treatment == "none")

Resulting targets

FileName    Name        Treatment       Time    Cy3

file1       ctrl_t0     none            0       control
file2       treat1_t0   treat1          0       sample
file5       ctrl_t1     none            1       control
file6       treat1_t1   treat1          1       sample
file9       ctrl_t2     none            2       control
file10      treat1_t2   treat1          2       sample
file13      ctrl_t3     none            3       control
file14      treat1_t3   treat1          3       sample
file17      ctrl_t4     none            4       control
file18      treat1_t4   treat1          4       sample

Linear regression*

time_factor <- factor(targetinfo$Time)
design <- model.matrix(~ 0 + time_factor)
fit <- lmFit(eset, design = design, weights = array.weights)

With this approach, the data is NOT normalized to the control as I do not know how to set the control samples as a reference. I saw a modelMatrix function in limma, which seems to do what I'm seeking, but I don't understand how I should combine it with my linear regression model ~ 0 + time_factor.

It is not necessary to subset the data. It is easier as well as more powerful to analyse all the data together.

On the other hand, this experiment has some special features and needs a bespoke analysis approach.

I think that assuming simple linear regression over 5 time-points is too restrictive. I would allow for a quadratic trend at the very least. I suggest allowing for linear and quadratic baseline trends like this:

X <- poly(Targets$Time, degree=2)
Baseline.Lin  <- X[,1]
Baseline.Quad <- X[,2]

Then setup relative trends for the treatments:

Treat1.Lin  <- (Targets$Treatment=="treat1" & Targets$Time>0) * X[,1]
Treat1.Quad <- (Targets$Treatment=="treat1" & Targets$Time>0) * X[,2]
Treat2.Lin  <- (Targets$Treatment=="treat2" & Targets$Time>0) * X[,1]
Treat2.Quad <- (Targets$Treatment=="treat2" & Targets$Time>0) * X[,2]
Treat3.Lin  <- (Targets$Treatment=="treat3" & Targets$Time>0) * X[,1]
Treat3.Quad <- (Targets$Treatment=="treat3" & Targets$Time>0) * X[,2]

The design matrix will be

design <- model.matrix( ~ Baseline.Lin + Baseline.Quad
                        + Treat1.Lin + Treat1.Quad
                        + Treat2.Lin + Treat2.Quad
                        + Treat3.Lin + Treat3.Quad )

After fitting a linear model

fit <- lmFit(y, design)
fit <- eBayes(fit, trend=TRUE, robust=TRUE)

you can test whether Treatment 1 affects expression relative to the baseline (control) by

topTable(fit, coef=c("Treat1.Lin","Treat1.Quad"))

You can test whether the trend is increasing or decreasing by

topTable(fit, coef="Treat1.Lin")

Similarly for the other two treatments.

Note that you do not need to explicitly normalize the treated expression to the controls. That adjustment is done automatically by the linear model.

Note also that if the time-points are not equally spaced, you should define Time to be the real time-points in days.