Hi again,

I already posted my question, but was asked to open a new one because it was suggested to use limma instead of DESeq2:
C: DESeq2 design formula: how to account for twin data (MZ and DZ)?

A summary of the relevant parts of my question:

I have to analyze some RNA-seq data, and I have got 400 samples of human whole blood (proportions of cell types available from FACS blood counts), all females, 40-86 years old. My phenotype of interest is disase as compared to healthy controls. There is also an intermediate stage, i.e. there are three disease states, but this intermediate stage is less interesting at the moment.

The problem: There are both monozygotic (MZ, n=220) and dizygotic (DZ, n=120) twins in the data. Additionally, there are singletons (n=60 samples with no twin pair). The different groups are not evenly distributed across my phenotype of interest, and the design is quite unbalanced in general:

• 10 disease samples (2 DZ twins, rest singletons)
• 20 intermediate samples (singletons)
• 370 healthy controls (MZ twins, DZ twins, singletons)

I would like to look at both differential expression, but also differential variability of expression (e.g. DiffVar) between the disease and healthy controls (I will probably have to use subsets of the data to have similar numbers of samples in each group, as I have only 10 disease samples, but 370 controls).

Additionally, we would like to have the possibility to look into twin effects as well, which is one of the reasons why I do not simply want to remove one of the two twins for each of them to have singletons only.

Thus I was wondering about how to correctly specify the design formula for these data.  Obviously - apart from covariates like cell subpopulations and perhaps age - I would like to take the twin pairs into account, as I expect a strong genetic effect. So, building blocks using the duplicateCorrelation() function could be useful, as also pointed out by Michael in my previous post. However, I would expect a much stronger effect for MZ twins than for DZ twins, so I do not think it is a good idea to treat all twin pairs equally (?) Adding the information of MZ and DZ as a covariate to the design formula will also not help, I suppose, as the effect occurs between the twin pairs, not across all individual samples within the group of MZs and DZs. It might work by adding an interaction term?

Any advice on how to build the blocks and/or define the design formula correctly would be highly appreciated.

More details about my doubts and concerns can be found in my original post (see link above).

Thank you very much for your help.

Simone

There's no reason why you need to use the same (sub)set of samples for differential expression and differential variability and testing for twin effects. So, I'm just going to address the issue of how to best perform a DE analysis with limma (assuming the counts have been voomed). Here's a few thoughts:

1. For starters, you should check whether the DZ and MZ twins actually have different correlations. This can be done by subsetting the normal samples into DZ and MZ subsets and running duplicateCorrelation on each subset separately. It may be that the correlations are very similar between the two twin types, in which case you don't have to worry.
2. If they are not similar, then you will just have to drop one twin, as duplicateCorrelation does not provide any facility for handling different random effects. I would probably drop one from each MZ pair, given that you would like to preserve as much DZ information in the few disease samples that you have.
3. Of course, duplicateCorrelation doesn't come for free. It requires more time to run and is slightly anticonservative. Given you have so many samples, I wonder whether you could obtain satisfactory results by simply ignoring one twin in any pair (DZ or MZ) and performing the DE analysis with singletons only. You can also try Ryan's suggestion, though I have no idea whether the correlation would be close to zero. I guess it depends on the contribution of genotype to expression.
4. There is no need to subset groups to equalize the sample sizes for a valid analysis with limma. Of greater concern - particularly for patient data - are hidden factors of variation in your groups. I can imagine that there is more variability within your normal samples due to genotype/lifestyle/etc., especially if disease samples are more prevalent in a subset of the population. This may require methods such as RUVseq or sva to identify and model such factors.