I am having trouble choosing a model matrix to analyze certain RNA-Seq data set. I will first describe the experiment and afterwards show between which designs I am doubting.
Design of the experiment
There are six distinct animals and each animal received the control (untrt) and the treatment (trt). The control and the treatment were separated by a fixed amount of time. This makes it a paired design, where each animal serves as it own control.
After the treatment, sperm was extracted from each animal. With that sperm sample, five embryos were generated. This results in 60 samples (6 animals times two treatments times 5 embryos = 60 samples).
This is a part of my metadata.
> head(setup, 10) animal treatment attribute sample21 animal1 untrt good sample22 animal1 untrt good sample23 animal1 untrt good sample24 animal1 untrt good sample25 animal1 untrt good sample26 animal1 trt good sample27 animal1 trt good sample28 animal1 trt good sample29 animal1 trt good sample30 animal1 trt good
This shows that the experiment is completely balanced.
> sapply(setup, table) $animal animal1 animal2 animal3 animal4 animal5 animal6 10 10 10 10 10 10 $treatment trt untrt 30 30 $attribute bad good 30 30
The goal of the experiment is to discover if the gene expression of the embryos differs between the treatments. But I am confused what to do with the different embryos? I suppose the embryos can't be considered biological replicates, because one sperm sample was used to generate five embryos.
My first idea was to sum the counts for all embryos from a given animal and a given treatment (for example sum all samples from animal1 and untrt). But I am not sure if it OK to do this. Don't you loose information about the variability between the embryos when summing ?
Keep the metadata and samples as they are and fit a model with animal as blocking factor (
~ animal + treatment).
In both cases, I would test the last coefficient of model.
Thanks in advance.
P.S. This is my first post on this forum, so I hope I followed all conventions.