Problem using funnel plot to detect outliers

Question

I am trying to detect outliers from a list of genes. Each gene has x number of SNPs, plotted on the X axis, and a score derived from the mean score of each SNP in a test, plotted on the y axis. My goal is to identify genes with the most extreme mean scores, given the number of SNPs they contain.

I don't want to set a hard cut-off of mean score > 0.5 because I will miss the longer genes that are clearly outliers for their size but which are not above this hard cut-off.

I have been using a funnel plot to try and identify outliers.

I tried to get a null distribution for detecting outliers. All SNPs in the dataset have a value between 0.1-1, in steps of 0.1 (eg 0.1,0.2,0.3 etc). I took all SNPs in genes in my dataset and calculated the proportion of SNPs in each of these bins. Which looks like this:

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.5399651   0.24848 0.1235064   0.0553278   0.0226341   0.0075721   0.0020886   0.0003944   2.94e-05    2e-06

As you can see, 0.1 is by far the most common score.

I then simulated the mean score for each gene by randomly drawing from this distribution to try and get a null distribution:

#read in the probability distribution shown above
ps <-read.table("probability_distribution.txt", header = TRUE, check.names = FALSE)

n <- 1e4
set.seed(42)

sims <- sapply(1:8954, 
           function(k) 
             rowSums(
               replicate(k, sample(x=(1:10)/10, size=n, replace=TRUE, prob=ps))) / k)

Then I calculate the quantiles to try and identify outliers from this simulated distribution:

quants <- apply(sims, 2, quantile, probs = c(0.025, 0.975))

plot((mean_score) ~ total_number_snps, data = DF)
matlines(1:80, t(quants), col = "red", lty = 2)

The problem is that I when I plot the quantile ranges for the funnel plot there are still way too many genes above the cut-offs I set.

normal scale:

X axis log scale:

This is particularly the case for the middle of the distribution. I think this is because in the real data the scores of the SNPs in a gene are not independent like in the simulation. Therefore the simulated data appears not to be a good way to generate the confidence bands.

Unfortunately my statistics background is not strong enough to know if there is another way I could generate confidence bands that more closely match the observed data.

Maybe there is a simple way to correct the mean score by the total number of SNPs per gene, but as it doesn't seem to scale linearly I am not sure about this.

The only solution I can think of so far is to take bins of sites along the x-axis and calculate quantiles for each bin using the empirical distribution, then take outliers from each bin. Does this solution sound reasonable?

Any help is greatly appreciated.

Thanks in advance!

For additional context this issue emerged from following the answer to a previous post:

https://stats.stackexchange.com/questions/175834/detecting-outliers-along-the-distribution-in-a-scatter-plot