Hello
I'm currently doing 30 samples of RNAseq anaylsis with DESeq, but now I have problem with model matrix again.
I posted my sample colData below.
donor disease batch donor.n
A Control B1 1
A Control B2 1
A Control B3 1
B Control B1 2
B Control B2 2
B Control B3 2
C Control B1 3
C Control B2 3
C Control B3 3
D Control B1 4
D Control B2 4
D Control B3 4
E muta B1 1
E muta B2 1
E muta B3 1
F muta B1 2
F muta B2 2
F muta B3 2
G muta B1 3
G muta B2 3
G muta B3 3
H mutb B1 1
H mutb B2 1
H mutb B3 1
I mutb B1 2
I mutb B2 2
I mutb B3 2
J mutb B1 3
J mutb B2 3
J mutb B3 3
I want to find disease samples' transcriptomic change without batch variation or donor variation, so I added donor.n to distinguish donor's effect.
Problem was that I had 4 donors samples in control group so I have to remove levels without samples as vignettes of DESeq says.
I tried with this code to make matrix model just like vignettes says.
m1<-model.matrix(~ disease+ disease:donor.n+ disease:batch, samples)
colnames(m1)
unname(m1)
all.zero <- apply(m1, 2, function(x) all(x==0))
all.zero
idx <- which(all.zero)
m1 <- m1[,-idx]
unname(m1)
It worked with DESeqDataSetFromTximport(txi,
colData = samples,
design = m1)
but when I tired to make the results, I couldn't find disease_muta_vs_Control
nor disease_mutb_vs_Control
in resultNames()
.
This is the results of the resultsNames()
resultsNames(dds)
[1] "Intercept" "diseasemuta" "diseaseControl" "diseasemutb.donor.n2"
[5] "diseasemuta.donor.n2" "diseaseControl.donor.n2" "diseasemutb.donor.n3" "diseasemuta.donor.n3"
[9] "diseaseControl.donor.n3" "conditionControl.donor.n4" "diseasemutb.batchB2" "diseasemuta.batchB2"
[13] "diseaseControl.batchB2" "diseasemutb.batchB3" "diseasemuta.batchB3" "diseaseControl.batchB3"
I can make results with res <- results(dds, contrast=list("diseasemuta","diseaseControl"))
to get muta versus Control but I cannot compare mutb versus Control.
Can anyone explain what I missed and how to make model matix of this situation?
Always thanks to everyone's kind response.
Thanks for the kind response.
Q1. Though, I did relevel already.
dds$disease<- relevel(dds$disease, ref = "Control")
and results was
Error in relevel.default(dds$disease, ref = "Control") : 'relevel' only for (unordered) factors
which means I set Control as reference level already.
Therefore, I can not get
diseasemutb
..Q2. Moreover I got
54 rows did not converge in beta, labelled in mcols(object)$betaConv. Use larger maxit argument with nbinomWaldTest
erroers so I did
dds <- nbinomWaldTest(dds, maxit = 1000)
,but it end up with the error
Error in terms.formula(design(object)) : argument is not a valid model
I dont know what to do at this point...
Q3. I wanted to perform likelihood ratio test, so did
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ disease:donor.n + disease:batch)
as you expected it shows the error
Error in DESeq(dds, test = "LRT", reduced = ~ disease:donor.n + disease:batch) : if one of 'full' and 'reduced' is a matrix, the other must be also a matrix
I know I am quiet bad at model.matrix, but can you explain what model should be used to reduce the donor and batch effect in this LRT codes?
I dont have any colleague who is familliar with DESeq2 or RNA seq to ask, so always thanks to any advice.
Thanks to your fabulous codes, packages and Ideas too.
Make sure that you are defining
disease
as an un-ordered factor."can you explain what model should be used to reduce the donor and batch effect in this LRT codes"
Unfortunately, I don't have time to provide this level of statistical advice on the support site but I have to limit myself to DESeq2 software related issues. I'd recommend collaborating with someone who can help design your full and reduced models.
Ahh now I see what I did wrong...
Thank you so much and sorry for wasting your time.
I think it is time to study full and reduced models.