DESeq2: Model matrix problem with 30 samples
1
0
Entering edit mode
dagurim • 0
@dagurim-22721
Last seen 11 months ago

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.

DESeq2 RNAseq • 283 views
0
Entering edit mode
@mikelove
Last seen 7 hours ago
United States

You'll need to set Control to the reference level of disease, then you will have diseasemutb in the names of the model matrix, which is the disease_mutb_vs_Control term (DESeq2 when you provide design as a formula renames these coefficients to make them more human readable than model.matrix, but these are the same terms).

0
Entering edit mode

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.

0
Entering edit mode

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.

0
Entering edit mode

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.