Thanks a lot and sorry for the delay response. It is really a big hand for me. Actually, this is a big concern of one reviewer of our work. I have read all the materials you provided, below is my table for construct model matrix. I would like to know whether it is right. Another question is how should I code if I want to know genes differently expressed between two Physiological states. e.g. use DESeq2
ddsMat_compare <- DESeqDataSetFromMatrix(countData = countdata_compare,
colData = design_compare,
design = ~ PhysiologicalState+PhysiologicalState:individual)
ddsruv_deseq = DESeq(ddsMat_compare)
???
or edgeR
y <- DGEList(counts=countdata_compare, group=conditions)
y <- calcNormFactors(y,method="TMM")
y <- estimateDisp(y, design_edge_used, robust=TRUE)
fit <- glmFit(y, design_edge_used,robust=TRUE)
glmLRT(fit, ???)
samples individual PhysiologicalState
4H_S1_P1 1 4H
4H_S1_P2 1 4H
4H_S2_P1 2 4H
4H_S2_P2 2 4H
4H_S3_P1 3 4H
4H_S3_P2 3 4H
48H_S1_P1 1 48H
48H_S1_P2 1 48H
48H_S2_P1 2 48H
48H_S2_P2 2 48H
48H_S3_P1 3 48H
48H_S3_P2 3 48H
AT_S1_P1 1 ACT
AT_S1_P2 1 ACT
AT_S2_P1 2 ACT
AT_S2_P2 2 ACT
AT_S3_P1 3 ACT
AT_S3_P2 3 ACT
TP_2011_S1_P1 1 TOP
TP_2011_S1_P2 1 TOP
TP_2011_S2_P1 2 TOP
TP_2011_S2_P2 2 TOP
TP_2011_S3_P1 3 TOP
TP_2011_S3_P2 3 TOP
TP_2015_S1_P1 4 TOP
TP_2015_S1_P2 4 TOP
TP_2015_S2_P1 5 TOP
TP_2015_S2_P2 5 TOP
TP_2015_S3_P1 6 TOP
TP_2015_S3_P2 6 TOP