Permutation strategy for testing coefficients of Generalized Linear Models

What is a good strategy for performing permutation test for coefficients of Generalized Linear Models (GLM) such as logistic regression? Specially when the response is binomial.

Following the questions here and here, and @Geln_b's informative answer (the former link), I was wondering if one could apply the same concept to logistic regression? If not, is there any other way to perform permutation test?

To provide more information about the problem I'm dealing with:

I have case-control cohorts of healthy/diseased individuals with respect to different diseases. The response variable $Y$ is binomial represting whether individual is healthy, and the predictors are genetic mutations.

For example, here is a model one would try to fit that includes pairwise interactions:

$Y = \beta_0 + \beta_1 \; mutation_1 + \beta_2 \; mutation_2 \; + \beta_3 \; mutation_1 * mutation_2 + \beta_4 \; age \; + \beta_5 \; gender \; + \beta_6 \; PC_1 + ... + \beta_9 \; PC_4 \;$

$PC_n$ are first four/five principal components to correct for population sub-structure e.g. ancestry (there are better ways to handle that, but doing this for now).

I'm interested in the interaction coefficient, and find it hard to set up a permutation strategy as I will be testing thousands of pairs of mutations separately (thousands of separate models - the matrix is fat $M << N$). And, I often hear that there is no solid/valid permutation strategy for interactions.

I was hoping to do something like what @GregSnow suggested here for a set of pairwise interactions (permutation distribution of multiple interactions in separate models as a multivariate gaussian distribution):

One possibility is to compute the mean vector from your permutation distribution (should be close to a vector of 0's) and the variance/covariance matrix. Then use that information to compute the Mahalanobis distance for all the points, the p-value would then be the proportion of points with a Mahalanobis distance greater than the original coefficients.

Since I'm dealing with diseases, I could use disease prevalence in population, and form the response as a probability instead of a binary variable $Y$. For example, disease A occurs in 0.1% of population, can we somehow use this as the response (like a probability of getting disease) instead of binomial variable, and would it make sense?