Hi BioC community,

I work on a PCR data set designed like that :

- 21 patients are included in the study

- For each patient, his gene expression profile was measured (by PCR) at the inclusion of the study on a sample of his back which was not treated. We call this gene expression measure the “baseline” measure, the time of this measure is called “T1”.

- After this baseline measurement, two areas are defined in the back of the patient : a left area (also called Z2) and a right area (also called Z1). A treatment (called T) and a control (= no treatment, called NT) are randomized between the two areas. So, each patient has one of the two following sequences :

           - Treatment in the left side of his back and no treatment in the right side (sequence called “T.NT” in the following)

           - No treatment in the left side of his back and treatment in the right side (sequence called “NT.T” in the following)

- His gene expression was then measured 8 days (called T8) and 22 days (called T22) after the baseline measurement on the two areas treated and no treated.

So, for each of the 21 patients, we have 5 samples (T1, T8Z1, T8Z2, T22Z1, T22Z2) on which we measured gene expression.

The gene expression was measured by PCR (TLDA technology), which gives the Ct of 64 genes, each gene measured in triplicate.

The variable to explain (=the response) is the evolution of gene expression relatively to T1.

So, briefly, the following steps were performed before fitting the model :

• Normalize the dataset to obtain DeltaCt for 63 genes (1 gene was used to normalize the others genes)
• Compute, the gene expression evolution as : DeltaCt_TimeX-DeltaCtT1 (where TimeX represents time8 or time22 according to the sample which is considered)

The first two questions are inter-product analysis, i.e comparison between the two products treatment and no treatment :

1/ Is the product effect significantly associated with the evolution (all times together) ?

2/ For each time (8 and 22) separately : Is the product effect significantly associated with the evolution ?

The last two questions are intra-product analysis, i.e analysis for each product separately :

3/ For each product (treatment and no treatment) separately : Is the time effect significantly associated with the evolution ?

4/For each product and each time (NT.8, NT.22, T.8, T.22) separately : is there a significant association with the evolution ?

These kind of design is very used in our field and an ancova model is routinely used to answer to the 4 questions, when we have only one response variable, for example a clinical variable.

Here is an example of the dataset we use in this clinical context:

> head(clinical)

  patientid.f time.f  product.f sequence.f side.f         T1           evolution

         117      8           T           T.NT       left       -0.1780001  0.2070001

         117      8           NT        T.NT       right     -0.1780001  0.5546665

         110      8           T           NT.T       right     -2.2918336  1.2368339

         124     22          T           T.NT       left       -1.7133338  -0.2925002

         102     22          NT        T.NT       right     -0.4601666  1.4905005

         101      8           T           NT.T       right     -1.0029996 -0.1920007

To be more precise, we fit the following model, 

> options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))
> library(lme4)
> library(lmerTest)
> model = lmer(evolution ~ T1 + product.f + side.f + time.f + product.f:time.f + sequence.f + (1|patientid.f), data=clinical)

where :

•  the evolution is the response
• T1 is the baseline measurement of the clinical variable (continous parameter allowing for taking into account the baseline level of the clinical parameter in the evolution).
• product is the product factor
• side is a blocking variable indicating the area where the sample was measured
• time is the time of the measure
• sequence is another blocking variable indicating one of the two sequence (NT.T and T.NT) presented above
• A random effect is set on the patient.

The answers of the 4 questions are respectively given by :

1/ anova(model, type=3) : product effect studied

2/ difflsmeans(model, test.effs="product.f:time.f")

3/ lsmeans(model, test.effs="product.f")

4/ lsmeans(model, test.effs="product.f:time.f")

So my question is mainly : how to transpose this analysis in the context of gene expression with limma, taking into account triplicates in measurement, covariate T1 depending on gene and missing values ?

Here is an extract of the eset used, with 2 genes in triplicates (each line represent a technical replicate, eset is ordered by gene and triplicates), and 5 samples. A cell is the evolution of gene expression relatively to T1. There are missing values.

> eset[1:6,1:5]

                P101T08Z01      P101T08Z02    P101T22Z01 P101T22Z02   P102T08Z01

EGFR.1    1.3973325         1.922667          -0.9410000    1.34133212     3.066335

EGFR.2    0.2103322         1.816666          -1.0289993    0.86233203     2.921334

EGFR.3    0.4623330         0.870667           0.8929996    0.78133265     1.499333

HER2.1   -1.9818335         0.820501           -1.6911662   -0.04483382    1.203333

HER2.2   -1.8958330         1.277500           -1.2571662    0.26116594    1.263334

HER2.3   -2.0888338         1.008501            NA               0.19516595            NA

And here the code I propose to answer the four questions :

> product.f[1:10]
[1] T  NT T  NT NT T  NT T  NT T
Levels: NT T
> side.f[1:10]
[1] right left  right left  right left  right left  right left
Levels: left right
> patientid.f[1:10]
[1] 101 101 101 101 102 102 102 102 103 103
21 Levels: 101 102 103 104 105 106 107 108 109 110 111 112 113 115 116 ... 124
> time.f[1:10]
[1] 8  8  22 22 8  8  22 22 8  8
Levels: 8 22
> sequence.f[1:10]
[1] NT.T NT.T NT.T NT.T T.NT T.NT T.NT T.NT T.NT T.NT
Levels: NT.T T.NT

> gp <- factor(paste(product.f,time.f,sep="."))
> gp<-factor(gp,levels=c("NT.8","NT.22","T.8","T.22"))

# Design matrix taking into account the two block variables and gp
> design <- model.matrix(~  side.f + sequence.f + gp)
> colnames(design)
#[1] "(Intercept)"      "side.fleft"     "sequence.fT.NT" "gpNT.22"        
#[5] "gpT.8"            "gpT.22"         
> colnames(design)[1]<-"Intercept"
> colnames(design)[4]<-"gpNT.22vsgpNT.8"
> colnames(design)[5]<-"gpT.8vsgpNT.8"
> colnames(design)[6]<-"gpT.22vsgpNT.8"

> library(limma)
> corfit <- duplicateCorrelation(eset,design=design,block=patientid.f)
> corfit$consensus
# [1] 0.5792043
> fit <- lmFit(eset, design,block=patientid.f,correlation=corfit$consensus)
> cm <- makeContrasts(

                        productEffect = (gpT.8vsgpNT.8)/2- (gpT.22vsgpNT.8-(gpNT.22vsgpNT.8))/2,
                        productEffectTime8=gpT.8vsgpNT.8,
                        productEffectTime22=gpT.22vsgpNT.8-gpNT.22vsgpNT.8,
                        timeEffectproductT=gpT.22vsgpNT.8-gpT.8vsgpNT.8,
                        timeEffectproductNT=gpNT.22vsgpNT.8,
                        levels=design)
> fit <- contrasts.fit(fit, cm)

# As recommended in the post DE analysis of PCR array [was: dataset dim for siggenes], I use robust=T and  trend=T in eBayes for PCR data
> fit <- eBayes(fit,robust=T,trend=T)

Here are finaly the answers to my questions :

1/ > topTable(fit, coef="productEffect")

2/

> topTable(fit, coef="productEffectTime8")

> topTable(fit, coef="productEffectTime22")

3/

> topTable(fit, coef="timeEffectproductNT")

> topTable(fit, coef="timeEffectproductT")

4/ ????

And so my questions are :

A/ In this model, I do not take into account for the covariate T1 because it depends not only on the patient but also on the gene studied. So how can I do to take into account the baseline expression in the model?

B/ I use duplicateCorrelation for the patientid, to transpose what is done in the clinical setting. In this case, how can I take into account the triplicates ? I tried with “patientid.f” as a block variable to take into account pairing, but sequence and patientid.f are counfounded (all samples from the same patient have the same sequence) and so it is not possible to keep the two informations. Taking patientid as a block variable, there is also the problem of missing values which make the design unbalanced, and in this case it is recommended to use duplicateCorrelation.

C/ How to obtain the answers to question 4, that is to say the coefficients: gpNT.8, gpNT.22, gpT.8 and gpT.22 ?

Many thanks for having read this long post, and thanks in advance for your answers,

Eléonore

Your design seems unnecessarily complicated. Why don't you just block on the patient ID directly in the design? Each patient has T1 and treated/untreated samples at both time points, so the comparisons can be performed within each patient. The design matrix would be something like:

design <- model.matrix(~ 0 + patientid.f + gp)

... where gp has all 5 levels (T1, NT.T8, T.T8, NT.T22, T.T22), releveled to have T1 as the baseline. This should give you one coefficient for each patient, plus an extra four coefficients representing the fold changes of the non-T1 samples over T1. If you're paranoid about the left/right sequence, you can split those 4 coefficients into two sets; one for the left/right patients, and another for the right/left patients. In any case, you can drop any number of these coefficients (or compare them to one another) to test for your product effect over time.

Now, onto your questions. For (A), the model I've suggested above will necessarily include T1. For (B), I'd suggest you just average the triplicate measurements and fit the model to the averages - from what I understand, these triplicates are technical replicates, and the true replication in this setting is that between patients. For (C), you can just compare each time/product combination to the T1 baseline within each patient, again using the above model. Not sure how useful that is, though, as you won't be able to tell whether the time or product effect is driving the increase over the baseline.