Unbalanced block design analysis with DESeq2
Hello,
I have RNAseq samples from 11 Patients. For each patient the samples were taken from two difference locations (Ileum and Caecum) and had two different inflammation status at those locations. However, the design is not balanced as I do not have 4 samples for each patient. Below is a description of the dataset Code:
1 Caecum NI I am interested in testing if there is a difference between the Location in general, the Status in general, the Status given the Location and the interaction between status and location. Therefore, I used DESeq2 with the following design Code:
~ Patient + Location + Status + Location:Status the effect of Location, Status and the interaction between the two. To check the difference for Location in general, I extracted the results of the DESeq analysis using contrast c("Location", "Ileum", "Caecum"). For Status, I used contrast = c("Status", "NI", "I") To test for the effect of Status given the Location, I used contrast = list("LocationCaecum.StatusI", "LocationCaecum.StatusNI"). To check for the effect of Location on Status (i.e. inflamed caecum has additional effect than Caecum alone and inflammation alone) I used contrast = list(c("StatusI", "LocationCaecum.StatusI"), c("StatusNI", "LocationCaecum.StatusNI")) Considering that the high amount of missing samples (or incomplete blocks), does the design formula I am using make sense? For instance, when comparing on Location, does the pairing have any value as it will only work for the following 8 samples (Or am I wrong?) Code:
Sample Location Status then do a DESeq2 analysis for each case separately. For example select all NI samples, then compare them between Ileum and Caecum with the design below Code:
design = formula(~ Patient + Status) Thank you very much in advance, Youssef 
Your model (~ Patient + Location + Status + Location:Status) has full rank, so it should work in general. You can drop at least patient #4, since it adds nothing.
Edit: BTW, make sure that "patient" is a factor (and reset the levels if you remove patient #4). 
Hi,
I have analyzed the samples according to Devon's recommendation and everything went well without any error messages. However, in the results, I do not understand why the interaction effect of Location and Status gives nearly the exact opposite log fold change. Is this something normal or it indicates something wrong with the analysis? Could it indicate a mislabeling of the samples? I am suspecting this already from the PCA plot of VST values. https://www.dropbox.com/s/g9jaey6zu1skeie/PCA.png Below is the relevant part of the code I used, starting from a SummarizedExperiment. Code:
# Remove sample 4 Code:
head(coef(dds)[, 12:19]) Thanks in advance, Youssef 
1 Attachment(s)
Below is the failed attachment from previous post
Attachment 3649 
Quote:
Long answer: Well, this will be long... This is normal, but a bit unfortunate. In DESeq2, we use a nonstandard way of setting up model matrices, which we call "extended model matrix". As explained in more detail in our paper, this is necessary to get proper shrinkage of logfoldchange estimates via ridge penality. In essence, using your example, and first explaining it without interaction: With standard model matrices, the intercept coefficient would be the expected log expression for a sample from patient level 1, location level 1 (ileum), status level 1 (NI), and then there are coefficients for all other patients, giving the differences in expression if the sample is any of the other patients, and one coefficient for the other location level (caecum) and one the other status level (I). The location coefficient, for example, is the difference between log expression for the other location level (caecum) and the first location level (ileum). For each factor, there is one less coefficient than there are levels. For extended model matrices, the intercept is the average over all levels rather then the value for all factors being the first level, and then, there is one coefficient for each level of each factor. These coefficients give the difference to the grand average (not to the first level as before). So, LocationIleum is the difference between ileum samples and the average over all samples, and LocationCaecum is the difference between caecum samples and the average. The difference caecumvsilleum is hence the difference between the LocationCaecum and the LocationIlleum coefficient  and this is what 'results' will calculate for you if you ask for this contrast. If you had a balanced design, these two coefficient would be exactly the same value with opposite signs (because the grand average would be in the middle of both), and half of what the coefficient with standard design matrices would be. As your design is not quite balanced, the values are only about the same, but this is not an issue. Now, for the interaction: An interaction is a difference of difference: In your case, it is the difference between (a) the difference between I and NI for caecum samples and (b) the difference between I and NI for illeum samples: interaction = (caecumI  caecumNI)  (illeumI  illeumNI) You could also ask for other interactions, e.g. interaction' = (caecumI  illeumI)  (caecumNI  illeumNI) and besides these two, there are two more ways to arrange the terms in this way. But all four of these are the same and differ only by sign, as you will notice quickly if you write the equations without parantheses. It is hence fully expected that all four interaction coefficients are the same, up to sign. There is only one interaction value  it's only due to the need for extended design matrices that the same information appears four times. And the lack of balances causes them to differ, but only slightly. Two final remarks:  The interaction in my double difference above is twice of what the coefficient says, again because these are the differences to the grand average, not to each other. (Mike: If you read this please doublecheck whether I'm right.)  Strictly speaking, we need extended design matrices only if we have more than two levels. Hence, normally, DESeq2 switches to standard design matrices (with socalled sum contrasts) if all factors have only two levels  to save users from exactly the confusion we have here. Unfortunately, your patient factor has more than two levels, and our code is not smart enough to realize that we need extended coding only this factor. 
Quote:
You need to figure out what the difference between the sample groups to the left and to the right is. Maybe there is some mislabelling, but it would have to be a rather catastrophic one, not just a single mislabelled sample. Or there is some dramaticv batch effect, such that something in the processing was done very differently in the left sample group than in the right one. 
Thank you for the detailed answer!

All times are GMT 8. The time now is 02:46 AM. 
Powered by vBulletin® Version 3.8.9
Copyright ©2000  2020, vBulletin Solutions, Inc.