Linear mixed model (repeated measurements) on HUGE dataset I have a huge dataset with measurements of 21k genes across 11k samples. The measurements, called expression, represent (roughly) how "active" that gene is in the cells of the sample (in this case tumor tissue from patients). In a given patient, a given gene may be mutated or not (wildtype).
The statistical question is whether having a mutation in gene X affects the expression of any other gene Y.
My questions are:

*

*Do the models I discuss below make sense?

*Will R handle a model this large?

*Because of the comparable number of samples and measurements I will need to employ some sort of shrinkage, right?

*Any other suggestions?

Because we are simultaneously measuring the expression for multiple genes in each patient, I believe this constitutes repeated measurements, so a linear mixed model could be used.
I have two tables:

*

*A mutation table tells me whether the type M_ij of gene i in patient j is mutant or wildtype.

*An expression table tells me the expression value X_ij of gene i in patient j.

A straightforward MANOVA would be:
X_1 + X_2 + ... + X_n ~ (M_1 + M_2 + ... + M+_m)^2

I'm not sure what is the right way to specify the mixed-effects model. Maybe
X_1 + X_2 + ... + X_n ~ (M_1 + M_2 + ... + M+_m)^2 + (1|M_1) + (1|M_2) + ... + (1|M_m)

I believe in a case such as this, where I am really comparing group means, I only need a random intercept model, right?
Finally, I still need to apply multiple comparisons as post-hoc analysis, right?
NOTE: The raw expression values are counts, which suggests a count model such as negative binomial. However, with such a large dataset I may be able to get away with a normal approximation. I would need to first normalize the gene expression across samples due to something called "read depth" which affects the scale of expression values. Also, there are packages that help with this, such as DESeq2 and limma-voom. These are details I will need to work out after I nail down my modeling scheme.
 A: It's a good idea to apply your knowledge of the subject matter intelligently first, before you focus on statistical analysis issues.
First, consider whether you really want to examine all mutations in all genes as predictors. Many genes are mutated only rarely. They are "passenger" mutations that don't really drive cancer but just happen to get mutated along the way as a tumor evolves. Even some genes that are mutated frequently, like TTN (coding for the muscle protein Titin), are frequently mutated simply because they are large or are replicated late in the cell cycle. You might be better off restricting yourself to documented cancer driver genes. See, for example, the Comprehensive Characterization of Cancer Driver Genes and Mutations published a few years ago.
Second, many differences in gene-expression values represent differences in cancer types rather than differences driven by gene mutations. Each cancer type originates from some normal tissue, so much of what gets expressed in a cancer type represents that normal tissue of origin. A breast cancer will have different gene expression than a liver cancer. And, as different types of cancer can have different sets of driver-gene mutations, what you find with a pan-cancer comparison of gene expression against gene mutations could just represent the corresponding cancer-type differences. A per-cancer-type analysis might be preferable.
Third, not all mutations have the same influence on gene expression. For example, in an analysis of TCGA data on 12 cancer types "mutation type (truncation versus amino acid-altering mutations) was the most important determinant of expression levels." Furthermore, many tumor mutations are changes in gene copy numbers, which can substantially affect gene expression, rather than point mutations. Such considerations really have to enter your analysis.
Fourth, you need to remove low-information noise before you start the analysis. Genes whose expression varies little among tumors of a certain type might well be removed prior to your analysis. Mutations that occur infrequently in a cancer type, even if they are in "driver" genes identified in other cancer types, might not be worth assessing.
Taking such steps will substantially reduce the scope of your analyses while improving interpretability. Only then should you worry about the statistical modeling mechanics.
Your concerns about negative-binomial modeling and relations between expression level and variance are well handled by established tools like the DESeq2 Bioconductor package that you mention and other helper packages like tximport. The way DESeq handles count data for differential expression analysis--including the individual samples explicitly in the negative binomial modeling--seems to make further per-tumor random effects unnecessary. For other types of analysis like clustering--which might be preferable for addressing your underlying biological questions--variance-stabilizing transformations of the counts are available.
Finally, even if you still end up with large numbers of gene-expression values and gene-mutation predictors, there are ways to proceed efficiently. Canonical correlation analysis is an established way to assess cross-correlations between multivariate data sets. Section 8.3 of Statistical Learning with Sparsity discusses sparse canonical correlation analysis, allowing that approach to be used on massive multivariate data sets. The bibliography cites a paper by Waaijenborg et al, which used that approach to analyze associations between gene expression and copy-number changes in cancer. You might find other ideas about how to proceed from the other sections of Chapter 8.
