eEtting up DESeq 2 analysis

[andrewelamb]

I have some RNA-seq data I'm trying to do DESeq2 analysis on, and it's more complicated than what I've done before.

Patients came in two different times and had their blood drawn and sequenced. There is some missing data, so that some patients only have data for visit 1 or 2. In between visits patients were given either a placebo or one of two treatments.

The answer we're trying to answer is what expression levels changed in between visits in a different way between the treatment groups. Based on the vignettes and other answered questions it seems the design would be:

~ patient + visit + treatment + visit:treatment

the interaction visit:treatment would be the 'difference in differences' of each treatment over the time between visits. Is this correct?

[dpryan]

You can simply exclude the patients for whom you only have a single sample. They'll get ignored in the analysis anyway. Anyway, yes your model is correct and you do indeed care most about the "visit:treatment" term.

[andrewelamb]

Thanks for the answer!

I did get this error however:

Error in checkFullRank(modelMatrix) :
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.

my pheno file looks like:

sampleName visit condition patient
1 V2 control 1
2 V5 control 1
3 V2 treatment 2
4 V5 treatment 2
5 V5 treatment 3
6 V2 treatment 3
7 V5 treatment 4
8 V2 treatment 4
9 V2 control 5
10 V5 control 5

Removing patients from the experimental design worked. Is there any way, or value, to preserve the patient data?

[dpryan]

Indeed, I should have foreseen that :P

If you were to instead use "~patient+condition:visit+visit" and got rid of the "conditiontreatment:visitV2" column in the model matrix then the result would work. The original problem was that each condition is comprised of a set of patients, so you can't have patient coefficients and a "condition" coefficient (which is just the average of the patient coefficients!).

Sorry that that's so confusing.

[andrewelamb]

Thank you for the help!

I apologize, I'm not entirely clear on how to set up my model matrix based on your answer. It seems I would still need every column if I were to use "~patient+condition:visit+visit".

[dpryan]

I had a typo in my reply, I meant to remove the "conditiontreatment:visitV2" from the model matrix. That'll make it full rank,

[andrewelamb]

Ahh I see, I'm getting my sample table and the model matrix confused.

So is this the correct way to use my own model matrix?

design_string <- "~patient+condition:visit+visit"
sample_table <- read.table(input_file, row.names = NULL, header = T, sep = ",")
deseq_object <- DESeqDataSetFromHTSeqCount(sampleTable = sample_table,
design = ~condition, #have to have something here
directory = count_folder)
mm <- model.matrix(as.formula(design_string), sample_table)
mm <- mm[,-19] # gets rid of conditiontreatment:visitV2
deseq_object <- DESeq(deseq_object, full=mm, betaPrior=FALSE)

[dpryan]

Something along those lines at least.