# Problem using funnel plot to detect outliers

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.

• Not very obvious if I don't have the data to experiment, but I think if you have the low and high quantile values for every total value somewhere (red points/line), you can see if for that total value the mean value is lower, between, or higher that those quantiles. – AntoniosK Oct 13 '15 at 21:52