[dmosk]

I see you figured it out while I was still typing my response

Thanks so much! I was actually also going to ask about whether betaPrior should be set to false in this case. I do have one more question, though: Is it valid to average over T1:G1 and T1:G2 to claim a main effect for T1, or is this inestimable given the coefficients in the model?

ETA: To clarify, I mean if I wanted to compare T1 to T2, could I average over the interactions similarly to how you showed in your example, or can I only compare T1 to T2 within G1 or G2?

[Michael Love]

Yes, this is a valid way to generate the T1 effect, and the fold change should end up being the average of the group2 common scale counts over the group1 counts for tissue T1. However, the dispersion estimates will be more accurate than with a model that doesn't include the patient effect, because in that case, patient specific differences would end up looking like higher dispersion. Another way this could be modeled is with mixed effects models, with the patient as a random effect. But this is how you should set up the contrast for modeling patient as fixed effect.

[lmolokin]

design not of full rank

Hi,
I'm attempting to use the same formula as above (~treatment + treatment:subject + treatment:time) to get within treatment time effects while accounting for subject differences. The problem I've run into is that the formula creates an effect that doesn't exist due to one less subject in the treated group vs control. In edgeR, I can remove the non-existent effect by simply removing the column from the design matrix but I'm not sure how to do that in DESeq2.

Here's the code I use for creating the dds:

dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = designframe,
                              formula(~tx+tx:subj+tx:time)

Is there a way to use model.matrix with the DESeqDataSetFromMatrix function so I can specify the truncated, full rank design?

Can I use ignoreRank for this situation?

Thanks.

[Michael Love]

Don't use ignoreRank, that'll lead to trouble.

I answered this in another thread today, should work for you:

http://seqanswers.com/forums/showthread.php?t=47881

In this development cycle, I will bump up the modelMatrix argument to higher level functions to make things simpler. I didn't realize until after we'd designed the interface that model.matrix() will create a matrix with a column of 0's.

[lmolokin]

Thanks for the quick reply.

I see how you use modelMatrix after already having created the dds object. But how can I initially create the dds object if I can't specify the model within it?

When I input:

dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = designframe,
                              formula(~tx+tx:subj+tx:time))

I get the error:
"Error in DESeqDataSet(se, design = design, ignoreRank) :
the model matrix is not full rank, so the model cannot be fit as specified.
one or more variables or interaction terms in the design formula
are linear combinations of the others and must be removed"

Thanks!

[Michael Love]

You can use any valid design except ~ 1, this design won't be used because you provide the model matrix at each step.

[lmolokin]

Thanks!

Just ran through the workflow without a hitch.

My next question is how do I specify different contrasts since resultsNames(dds) now only gives me just [1] Intercept?

[Michael Love]

hi,

I just tried this out and found a bug, associated with using a design of ~ 1. You can use any valid design other than this one (~ tx), and it should work. I plan to support this functionality at a higher level in the next release, and in the process will add unit tests for situations like this!

[lmolokin]

Thanks for the help thus far, Michael.

Is there a way to insert the modified model.matrix into the dds object in order to test contrasts from modelMatrix=m?

In other words, so that when I call resultsNames(dds) I will see all of the interactions terms from ~tx+tx:subj+tx:time instead of just what the formula placeholder (~tx) gives?

[Michael Love]

Unfortunately, this trick won't work beyond the simple LRT results table for the current release. Contrasts on user provided model matrices aren't possible with the current release code. Update 2015.01.05 below

It's relatively easy for me to implement, but this has to be done in the development branch. I will likely implement this in Nov or Dec.

[lmolokin]

I see.

Is there no way to specify a model matrix within the DESeqDataSetFromMatrix function in lieu of a formula?

[Michael Love]

I tried an example, and you can use numeric or list style contrasts after nbinomLRT with a user-supplied model matrix in version 1.6. Like so:

res2 <- results(dds, test="Wald", contrast=c(0,0,1))
res2 <- results(dds, test="Wald", contrast=list("conditionB"))

[lmolokin]

I'm not sure I understand where the "user specified model matrix" is applied to dds.

I create the dds object with a place-holder formula here:

dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = designframe,
                              formula(~tx))

Then create the actual model and supply it for dispersion estimations here:

m <- model.matrix(~tx+tx:subj+tx:time)
colnames(m)
m <- m[,-24] # get rid of the 0's column
dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds, modelMatrix=m)
dds <- estimateDispersionsFit(dds)
dds <- estimateDispersionsMAP(dds, modelMatrix=m)
dds <- nbinomWaldTest(dds)

But then inputting:

res_ct2.ct1 <- results(dds, test="Wald", contrast=list("txc.time2","txc.time1"))

...isn't valid because resultsNames(dds) doesn't contain interactions from modelMatrix=m, it only has the elements given by the place-holder formula ~tx.

[Michael Love]

Ah, now I see the problem. You can't use nbinomWaldTest(dd) in the last step. You need to use all the lines of code from the thread I pointed you to, which was Update 2015.01.05 below

dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds, modelMatrix=m)
dds <- estimateDispersionsFit(dds)
dds <- estimateDispersionsMAP(dds, modelMatrix=m)
dds <- nbinomLRT(dds, full=m, reduced=m[,-idx])

where idx is a numeric vector of columns to remove from the model matrix for the LRT. There is no support with nbinomWaldTest in the current version to supply model matrices, but I've promised to work on this in devel.

[lmolokin]

You're right! It appears to have worked now.

I received this message after running the LRT: 32 rows did not converge in beta, labelled in mcols(object)$fullBetaConv. Use larger maxit argument with nbinomLRT.

Does this suggest that I need to rerun it with different nbinomLRT parameters? Not sure what to make of the message.

Thanks.