Is this logistf approach the right way to deal with complete separation in my logistic regression?

Question

I am performing lots of association tests between genotypes and binary diseases. The genotypes can be very rare, and my tests often have extreme case control imbalance (e.g. 500 controls for every case among 50K samples). When running logistic regression using glm() in R, my test often produce the error: glm.fit: fitted probabilities numerically 0 or 1 occurred.

Here is some dummy data that replicates the issue. There are more predictors than a single "target genotype" in my real data (multiple genotypes at a multi-allelic site), but I leave them out here for simplicity:

set.seed(001)
reference_genotypes_control <- sample(c(rep(0,5),rep(1,495), rep(2,39500)))
target_genotypes_controls <- sample(c(rep(0,39750),rep(1,249), rep(2,1)))
disease_state_controls <- rep(0,40000)

controls <- data.frame(reference_genotypes=reference_genotypes_control,
                       target_genotypes=target_genotypes_controls,
                       disease_state=disease_state_controls)

reference_genotypes_cases <- sample(c(rep(1,5),rep(2,495)))
target_genotypes_cases <- rep(0,500)
disease_state_cases <- rep(1,500)

cases <- data.frame(reference_genotypes=reference_genotypes_cases,
                       target_genotypes=target_genotypes_cases,
                       disease_state=disease_state_cases)

my_data <- rbind(controls,cases)

Genotypes can be 0,1, or 2. I am interested in the impact of the target genotype on the binary disease state. I leave reference genotype out of the model specification as reference.

Regular old logistic regression gives the glm.fit error:

glm(disease_state ~ target_genotypes, data=my_data, family = 'binomial')
summary(regular_logistic)

which gives:

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.1581  -0.1581  -0.1581  -0.1581   2.9625  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -4.376      0.045 -97.235   <2e-16 ***
target_genotypes  -13.180    249.390  -0.053    0.958    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5388.3  on 40499  degrees of freedom
Residual deviance: 5382.0  on 40498  degrees of freedom
AIC: 5386

Number of Fisher Scoring iterations: 16

Am I having issues because the target genotype is always 0 in the cases and this is considered complete separation issue? My coefficients and standard error are very large and p value is close to 1 for target_genotypes, which I have read is a side effect of complete separation. I also checked complete separation with brglm2.

library("brglm2")
glm(disease_state ~ target_genotypes,data=my_data,family=binomial("logit"),method="detect_separation")

which gives:

**Separation: TRUE** 
Existence of maximum likelihood estimates
     (Intercept) target_genotypes 
               0             -Inf 
0: finite value, Inf: infinity, -Inf: -infinity

I read that logistf can work with data that produces the separation problem, so I tried to run:

logistf(disease_state ~ target_genotypes, data = my_data, firth=TRUE,pl=TRUE)

which gives:

Model fitted by Penalized ML
Coefficients:
                      coef   se(coef) lower 0.95  upper 0.95   Chisq          p
(Intercept)      -4.374772 0.04498013  -4.464213 -4.28785813     Inf 0.00000000
target_genotypes -1.836535 1.41644734  -6.670910  0.07811345 3.41453 0.06462481
                 method
(Intercept)           2
target_genotypes      2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

Likelihood ratio test=3.41453 on 1 df, p=0.06462481, n=40500
Wald test = 1.681115 on 1 df, p = 0.1947764>

Is the logistf approach I show here an appropriate way to deal with my apparent separation problem?