What can I do if I want to apply Fisher's method (or some other method of combining p-values), but my p-values have unknown (positive) dependence? I have p-values from several tests (k~50) that I'd like to combine. I don't have a model of their correlation structure. They aren't independent. In the worst-case scenario, my effective sample size might be more more like 15.
I'm dealing with marginally significant results (About 10% remain significant after a two-stage Benjamini Hochberg FDR correction for α=0.05). So, I need to be very careful to both avoid a false positive and to avoid wasting statistical power.
Option 1: Bonferroni argument
The Bonferroni correction is considered conservative when considering multiple tests, where any single positive test is taken as a positve result.
You usually divide your p-value threshold by $k$ for this, but it is equivalent to multiple the p-value by $k$. This is really a linear approximation when $p\ll 1$. The real quantity you want to use is $\tilde p = 1-(1-p)^k$, which gives you the probability of seeing at least one test pass under the null distribution.
By this argument, if I think that my effective sample size is no less than $k/3$, I could correct my p-values to $\tilde p = 1-(1-p)^3$, then apply Fisher's method. On my data, this gives something like $p_{combined}$ = 5×10^-3.
Option 2: Fisher's method argument
My understand is that Fisher's method: (1) converts the p-values to the sum-of-squares $\sum_{i=1}^k\|z_i\|^2$ quantity of a χ² distribution (2) sums these up, then (3) tests whether the result looks significant for a χ² with DoF = $2k$. On my data, this gives something absurd like 1×10^-11, which is certainly spurious.
My intuition is this: Say that my data secretly consisted of individual tests replicated three times each. Each of this would give a (possibly noisy) estimate of the same p-value.
I could hypothetically average these together into $k/3$ tests, then apply Fishers method, testing against a χ² with DoF = $2k/3$.
Instead of testing whether $\sum_{i=1}^k\|z_i\|^2$ looks significant with a χ² with $2k$ DoF, I coulld test whether $(\sum_{i=1}^k\|z_i\|^2)/3$ looks significant with a χ² with $2k/3$ DoF. On my data, this gives $p_{combined}$ =  1×10^-5.
Option 3: Preprocessing with Benjamini Hochberg
Say I trust the corrected p-values froma two-stage Benjamini Hochberg FDR correction (α=0.05). I think this is fine, since this FDR correction depends on the distribution of p-values, and doesn't require independence? (correct me if this is wrong).
I can run FDR correction first, then apply Fisher's method to combine the corrected p-values. On my data, this gives $p_{combined}$ = 2×10^-3.
Question: Are any of these methods considered kosher and, if so, does anyone know a reference to back them up? (Other solutions very welcome, of course).

Combining p-values within statistical framework is a flourishing research area, mostly in genetics or in public health, when dealing with meta-analyses for example.
Basically, p-value combination can be reached using two kind of methods, asymptotic and resampling-based approaches, respectively. I assume you are comfortable enough with these families of methods to move forward and present 2 good options for your case (3 actually, but the latter is just a combination of the two first). I will focus my attention on the ACAT (Liu et Lin, 2019) and minP (Westfall & Young, 1993) procedures, while providing more a practical guide  focusing on general concepts. I strongly recommend you to read original papers if you want technical details.
Aggregated CAUCHY Association Test a.k.a ACAT
ACAT is a very recent method, initially proposed by Liu & Xie, 2019. Basically, the method is similar to Fisher's method since it is built upon a linear combination of weighted p-values. Following the original paper, correlation structure does not impact much the Type I error rate because of the heaviness of the Cauchy distribution tails. One other argument in favor of ACAT is to maximize the power by not adjusting any degrees of freedom of the test (unlike Fisher's method). Finally, by providing elegant asymptotic results from the Cauchy distribution (the sum of standard Cauchy random variables is a standard Cauchy), computation times are fast which is very advantageous in practice.
The method can be implemented in your favorite programming language, but if you are working in R, the method is available here: https://github.com/yaowuliu/ACAT
minP
minP is a re-sampling procedure. minP is a two-stage method, using minimal p-value as test statistic (Dudoit & Van der Laan, 2008). Roughly speaking you need to resample your data to obtain your test statistic then resample them to obtain your p-values. Two advantages I see here, the resampling nature of the procedure respects your correlation structure observed in your original data (even if it's very strong), while maximizing the power by choosing the minimal p-value. Sandrine Dudoit has written lot of good papers on standard permutation-based approaches. It is worth noting that in presence of a large number of p-values, the method can be practically unfeasible, even though adaptations can be used in order to reach faster results (kurtosis adapted approaches for example, Lee et al., 2014).
Hybrid approach , MinP-CCT-MinP
This is the more recent method. Since I did not work with it, I will just mention that Chen, 2022 showed that minP is more powerful than ACAT under certain correlation structures, such as autoregressive models for example.
The code is available here: https://github.com/zchen2020/Robust-P-value-combination-tests
Finally I provide a short guide of good practices when dealing with combination of p-values.
Good practices
Generally speaking, in many cases, these approaches will be enough to deal with a large set of problems.

*

*If you have a large number of p-values and you do not expect large correlations across your p-values, ACAT should be a good option. You will obtain fast results with a good control in your distribution tail (small p-values).

*However, if you do expect strong correlations, ACAT should be not enough. Using resampling-based approaches seems the way-to-go.

N.B: Recent methods use saddlepoint approximation to reach more control in distribution tails. From my knowledge, the method is not implemented for general purpose, but only for very specific applications.


*Finally, if you expect specific correlation structures, probably using hybrid technique will lead to powerful results with good Type I error control.

Anyway, I could illustrate the pros and cons of the methods simulating some data. Let me know if you are interested.
Hope this helps.
