I am working with RNAseq analysis.. and I would like to get confirm for my model design before doing the analysis.
I have 8 samples from 4 patients (each patient has two samples before treatment and after treatment.. ) I would like to find DEG genes which are differentially expressed by treatment.. However, I suspect that it also be affected by sex factor. so I would like to consider the sex factor (but among 4 patient, I only have one male, three female)
Here is the exact description for data model.
sample1 patient1 pre_treatment female
sample2 patient1 post_treatment female
sample3 patient2 pre_treamment male
sample4 patient2 post_treament male
sample5 patient3 pre_treamment female
sample6 patient3 post_treament female
sample7 patient4 pre_treamment female
sample8 patient4 post_treament female
df_enzyme <- factor(rep(c("pre", "post"), len=8))
df_patient <- rep(1:4, each=2, len=8)
df_sex <- factor(c("F","F","M","M","F","F","F","F"))
df_enzyme <- relevel(df_enzyme, ref="pre")
df_meta <- data.frame(row.names=colnames(data), patient = df_patient, enzyme = df_enzyme, sex = df_sex)
# sincce I don't have replication, for the dispersion estimation, I used blind..
d_fabry <- estimateDispersions(d_fabry, method="blind", sharingMode="fit-only")
I made the following model.
d_fit0a = fitNbinomGLMs(d_fabry, count ~ patient)
d_fit0b = fitNbinomGLMs (d_fabry, count ~ enzyme)
d_fit0c = fitNbinomGLMs (d_fabry, count ~ sex)
d_fit1a = fitNbinomGLMs(d_fabry, count ~ patient + enzyme)
d_fit1b = fitNbinomGLMs(d_fabry, count ~ patient + sex)
d_fit2 = fitNbinomGLMs(d_fabry, count ~ patient + enzyme + sex)
d_fit3 = fitNbinomGLMs(d_fabry, count ~ patient + enzyme + sex + enzyme:sex
for the treatment(enzyme) effect, which one is right?? maybe the first one?? ( it assume that sex would have effect.. right?)
dh_pvalsGLM = nbinomGLMTest( dh_fit2, dh_fit1b )
dh_pvalsGLM = nbinomGLMTest( dh_fit1a, dh_fit0a)
for the treatment and sex interaction effct, (since I suspect that depending on different sex, effect of treatment might be different or not), the following modle is right??
dh_pvalsGLM = nbinomGLMTest( dh_fit3, dh_fit2 )
Again, I would like to find DEG which are responsive to treament(enzyme), but I suspect that depending on sex, the responsiveness might be different (maybe not) I also, test with this assumption. Does my model do the right thing for this??