I have problems performing a time-series experiment with DESeq2 while accounting for age as a covariate. We have one group with each individual having three time points. We're interested in detecting differential effects between time points while accounting for within individual variability and age as a covariates. The following code runs just fine but once adding age into the equation DESeq2 gives an error that the model matrix is not full rank. I suspect the problem being with variables SUBJECT_ID and AGE being linear combination of each other based on the DESeq2 vignette.
Does anyone have solution how this time series experiment could be done so that age could be accounted for? Was wondering that would replacing SUBJECT_ID with AGE correct the problem as subject_id is an arbitrary number designated for each individual while age is specific. The study is less than a year so basically the age is constant at each three time point.
Thank you so much in advance!
######################## DIFFERENTIAL ANALYSIS with DESeq2 Case-only #################
########################################################################################
#Lets check that the order is right
coldata <- subset(coldata, coldata$CASE_CONTROL==1)
coldata <- coldata[order(coldata$CASE_CONTROL),]
sample_list <- as.vector(coldata$SAMPLE_ID)
countdata <- countdata[,sample_list]
str(coldata)
#Lets create the modelmatrixes for full and reduced model
#Full
x <- model.matrix(~ PRE_MID_POST + SUBJECT_ID, coldata)
all.zero <- apply(x, 2, function(x) all(x==0))
all.zero
idx <- which(all.zero)
full <- x[,-idx]
full1 <- unname(full)
#Reduced
x <- model.matrix(~ SUBJECT_ID, coldata)
all.zero <- apply(x, 2, function(x) all(x==0))
all.zero
idx <- which(all.zero)
reduced <- x[,-idx]
reduced1 <- unname(reduced)
##################################################################################
############# Lets do the DE with LRT and both modelmatrices #####################
##################################################################################
source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")
library(DESeq2)
ddsFullCountTable <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata, design = ~1)
ddsFullCountTable <- ddsFullCountTable[ rowSums(counts(ddsFullCountTable)) > 1, ]
keep <- rowSums(counts(ddsFullCountTable) >= 5) >= 5
ddsFullCountTable <- ddsFullCountTable[keep,]
dds <- DESeq(ddsFullCountTable, test="LRT", full=full, reduced=reduced)
