Method principle

In brief, this method takes as input the SFSs of WS mutations (yWS) and SW mutations (ySW). From these SFSs, it computes the inversed SFSs:

• zWS = rev(yWS)
• zSW = rev(ySW)

which are used to create two variable:

• u = zWS / ySW
• v = zSW / yWS

for more precision, additional variables can be computed u^2, u^3, u^4 … v^2, v^3, v^4

Then one computes the multiple regression:

log(yWS/ySW) ~ x + u + v + u^2 + v^2 + u^3 +v^3

where x is the derived allele frequency (i/n).

The coefficient for x corresponds to B (like in the standard approach) and terms u and v account for polarization errors.

Method evaluation

To evaluate this method, I simulated datasets subject to different levels of gBGC (B ranging from 0 to 2), with levels of polarisation errors similar to those reported in humans. I then estimated B using multiple regression (MethodSG) with different numbers of factors (nf):

• MethodSG nf=1: log(yWS/ySW) ~ x + u + v
• MethodSG nf=2: log(yWS/ySW) ~ x + u + v + u^2 + v^2
• MethodSG nf=3: log(yWS/ySW) ~ x + u + v + u^2 + v^2 + u^3 +v^3

For a comparison, I also estimated B by simple regression:

• Method0: log(yWS/ySW) ~ x (i.e. without accounting for polarization errors)
• Method1: log(yWS/ySW) ~ x (computed on SNPs with a DAF < 0.9, to exclude the subset of SNPs that are the most strongly affected by polarization errors)

I also applied these methods on human popgen data (from the 1000 genomes project) to compare with the results obtained with the ML approach (Glémin et al. 2015).

Simulation parameters

I first ran simulations with polarization error rates corresponding to CpG sites (i.e. high polarization error rate for SW mutations), and a GC-content corresponding to synonymous codon position (i.e. GC-rich sequences).

# Polarization error rates for WS and SW mutations:
eSW=0.04 # matches observed error rate at CpG sites in humans
eWS=0.01 # matches observed error rate at other sites
# Ratio of mutation rates: SW/WS  and WS/(SS+WW)
rSW_WS=2
rWS_SSWW=2
# GC-content (match synonymous sites)
GC_content=0.65
# Range of B values
lB=c(0.001, 0.25, 0.5, 1, 2)
# Number of simulation replicates per set of parameters
nbReplicate=10

NbClassDAF=20

NbClassDAF=100

If we increase the number of bins of DAF, this does not improve the estimates of B. On the contrary when the number of SNPs is limited, then MethodSG becomes less accurate and more biased than Method1.

NbClassDAF=10

If we consider only 10 bins of DAF, MethodSG with nf=1 provides estimates of B that are more noisy than with 20 bins of DAF, but that remain unbiased even with limited SNP sample size.

Summary of simulations: step 1: CpG sites, GC_content=0.65

I provide below a summary of the accuracy of the different methods for this set of parameters (mimicking CpG sites at synonymous codon positions). I evaluated the accuracy according to two criteria:

• Correlation between the estimated and the true B => measure to what extent the variation in the estimated B reflects true variation
• Mean difference between the estimated and the true B => measure to what extent the estimates of B are biased

As expected, methodSG is sensitive to the number of bins of DAFs (NbClassDAF) and also to the median number of SNPs per DAF bin. The lower the number of SNPs, the higher the level of noise in methodSG estimates. Increasing the number of factors (from nf=1 to nf=3) leads to increase the level of noise.

For a given number of SNPs per DAF bin, increasing the number of DAF bins leads to stronger correlations between the estimated and true B. However, increasing the number of DAF bins implies that the number of SNPs per DAF bin is decreased. There is therefore a trade-off between increasing the number of DAF bins (to better capture the shape of the SFS) and increasing the number of SNPs per DAF bin (to reduce sampling noise).

Both Method0 and Method1 provide B estimates that are strongly correlated to the true B, even for small SNP datasets. Thus, eventhough these methods provide biased estimates (see below), they may be useful to capture variation of B along the genome (but see below for more details regarding biases linked to variation in GC-content).

As expected, Method0 overestimates B, and this overestimate increases as the number of bins of DAFs decreases.

Method1 is less biased than Method0, but still leads to overestimates B.

MethodSG is able to provide unbiased estimates of B, but only when the number of SNPs per DAF bin is sufficiently large. For a median number SNPs per DAF bin below 500, MethodSG tends to overestimate B (even more than Method1 for small SNP datasets). Overall MethodSG with nf=1 tends to perform better than with nf=2 or 3.

For this set of parameters, increasing the number of factors (from nf=1 to nf=3) does not reduce the bias in the estimate of B (on the contrary!).

Summary of simulations: step 2: non-CpG sites, GC_content=0.65

I then ran simulations with polarization error rates corresponding to non-CpG sites (i.e. with eSW=eWS=0.01), and a GC-content corresponding to synonymous codon position (i.e. GC-rich sequences).

# Polarization error rates for WS and SW mutations:
eSW=0.01 # matches observed error rate at non-CpG sites in humans
eWS=0.01 # matches observed error rate at non-CpG sites in humans

# Ratio of mutation rates: SW/WS  and WS/(SS+WW)
rSW_WS=2
rWS_SSWW=2
# GC-content (match synonymous sites)
GC_content=0.65

In terms of measuring variation in B, the results are similar to previous simulations: MethodSG is more noisy than Method0 or Method1 when the number of SNPs is limited.

For non-CpG sites, polarization errors are weaker, and hence B estimates are less biased. In particular, Method1 performs relatively well. MethodSG outperforms Method1 only when the number of SNPs is high

Summary of simulations: step 3: CpG sites (GC_content=0.4)

I then ran simulations with polarization error rates corresponding to CpG sites (i.e. with eSW=0.04 and eWS=0.01), and a GC-content corresponding to the genome-wide average (40% G+C).

# Polarization error rates for WS and SW mutations:
eSW=0.04 # matches observed error rate at CpG sites in humans
eWS=0.01 # matches observed error rate at other sites
# Ratio of mutation rates: SW/WS  and WS/(SS+WW)
rSW_WS=2
rWS_SSWW=2
# GC-content (match synonymous sites)
GC_content=0.4

Again, in terms of measuring variation in B, the results are similar to previous simulations: MethodSG is more noisy than Method0 or Method1 when the number of SNPs is limited.

Again, MethodSG is less biased than Method0 or Method1 when the number of SNPs is sufficiently large

Summary of simulations: step 4: non-CpG sites (GC_content=0.4)

Finally, I ran simulations with polarization error rates corresponding to non-CpG sites (i.e. with eSW=eWS=0.01), and a GC-content corresponding to the genome-wide average (40% G+C).

# Polarization error rates for WS and SW mutations:
eSW=0.01 # matches observed error rate at non-CpG sites in humans
eWS=0.01 # matches observed error rate at non-CpG sites in humans

# Ratio of mutation rates: SW/WS  and WS/(SS+WW)
rSW_WS=2
rWS_SSWW=2
# GC-content (match synonymous sites)
GC_content=0.4

Performances are unchanges in terms of accuracy.

But for this set of parameters, all methods are unbiased. Given that MethodSG is more noisy, Method0 and Method1 are better than MethodSG.

Impact of GC-content

The impact of polarization errors on the SFS (and hence on B estimates) depends on the difference between the number of WS vs SW mispolarized SNPs. It therefore depends on eWS, eSW, rSW_WS and on the GC-content of sequences. Given that the GC-content varies widely along the genome, this might induce variation in the estimate of B.

To assess how the different methods are affected by GC-content, I performed simulation for sequences of different GC-content (from 0.3 to 0.8, to mimick variation in GC3 at synonymous codons positions), with 3M, 300k or 100k SNPs. I set NbClassDAF=20.

I considered a genome composed of 5 compartments (each corresponding to 20% of of the genome) with B=0.001, B=0.25, B=0.5, B=1, B=2. I computed the difference between the true B and the inferred B in each compartment, and then measured the mean of these differences.

I first considered a situation where eSW=eWS=0.01 and where rSW_WS=1 (i.e. no mutational bias, same error rate for WS and SW mutations):

All methods tends to be biased when the GC-content of sequences is far from the mutational equilibrium (50%) => this induces an artefactual correlation between B and GC-content.

MethodSG is able to provide unbiased estimates, but only if the number of SNPs is very large. With smaller SNP sample sizes, Method1 is more robust than MethodSG.

I then considered a more realistic case, with rWS_SSWW=2, and with error rates corresponding to that at CpG sites (eWS=0.01, eSW=0.04).

Both Method0 and Method1 overestimate B, and this bias increases with sequence GC-content => this induces an artefactual correlation between B and GC-content.

MethodSG with nf=1 provides (relatively) unbiased estimates of B for a wide range of GC-content, but only when the median number of SNPs per DAF bin is sufficiently large. When the number of SNPs is limited (~300 SNPs per DAF bin), MethodSG with nf=1 tends to overestimate B.

Finally, I considered a case with rWS_SSWW=2, and with error rates corresponding to that at non-CpG sites (eWS=0.01, eSW=0.01).

Again, with a large number of SNPs, MethodSG with nf=1 provides unbiased estimates, whatever the GC-content. But for smaller SNP dataset, Method1 performs as well as MethodSG (except for very GC-rich sequences, where MethodSG is still better than Method1).

Genome-wide B, NbClassDAF=20

I first estimated the genome-wide B, using either the entire dataset of non-coding SNPs, or CpG and non-CpG SNPs separately. I set the number of DAF bins (NbClassDAF) to 20.

I show below the observed SFSs and, and the corresponding plots [log(WS/SW) vs DAF] from which B is inferred:

NbClassDAF=20

Genome-wide B, NbClassDAF=100

I repeat the analysis with NbClassDAF=100.

NbClassDAF=100

Comparison with Glémin et al. 2015 (Table 1)

The table below provides estimates of B obtained by maximum likelihood (M1e = M1* = model M1 with polarization errors) or by Method1 or MethodSG (with nf=1, 2 or 3), for NbClassDAF=20:

Same table for NbClassDAF=100:

Overall, genome-wide estimates of B given by method SG (with nf=1 or nf=2) are pretty close to ML estimates, notably at CpG sites (which are the most affected by polarization errors). Increasing NbClassDAF from 20 to 100 does not improve much the match between ML and method SG estimates.

Estimate B per genomic window

rec="hapmap"
#rec="bherer"

NbClassDAF=20
nbinsNC=500

nbinsNC=50

Bins of contiguous SNPs: CpG vs non-CpG

In their Fig. 5, Glémin et al (2015) compared the values of B at CpG or non-CpG sites, estimated with model M1 (i.e. without accounting for polarization errors, Fig5B) or M1* (i.e. accounting for polarization errors, Fig5D). With M1*, they observed that for a same recombination rate, CpG and non-CpG sites have similar values of B (contrarily to model M1, where B is over-estimated for CpG sites):

Fig. 5BD, Glémin et al (2015)

I checked whether we could retrieve this result with method SG. For this, I classified non-coding SNPs into 50 or 500 bins of contiguous SNPs. I computed B on each bin, and I plotted B vs recombination rates.

nbinsNC=50

nbinsNC=50

=> Method1 clearly overestimate B at CpG sites (relative to the ML estimate). MethodSG gives results that are closer to ML estimates (notably with nf=2).

nbinsNC=500

nbinsNC=500

=> with nbinsNC=500, MethodSG with nf=1 seems to be as biased as Method1. With nf=2 or 3, estimates appear to be very noisy.