2-way ANOVA-like analysis with binomial RVs?


I have a biological dataset and I'm trying to see if certain treatments are correlated with significant changes.

In each sample, I'm measuring 20 binomial RVs: B(p,n). For those curious, these are reads of sequence data (DNA), which can either be an A or a G (hence binomial). Over these 20 binomial RVs, p is expected to vary anywhere from 0 to 1 (ie. not always near 0.5). In almost all cases, n is very large, like 10^4.

I have 16 "conditions" total (4 treatments X 4 replicates). I'd like to use something like a 2-way ANOVA to see if the treatments significantly impact the 20 binomial RVs, while also taking into account the variability between replicates. Restated: Do any of the treatments significantly change the binomial variable p given the variation seen in the replicates? I'd like to answer this question for all 20 binomial RVs. (And for each RV, we have ~10^4 independent samples in each of the 4x4=16 conditions.)

I don't think a 2-way ANOVA (with replicates) is the right tool because my RVs are not normally distributed. But the general idea there is correct (2 factors: 1) treatment; 2) which of the 20 sites).

Thank you in advance for any ideas!