Estimation of L in Stairway Plots using RAD-seq

[Clemens Maylandt]

Dear Ryan,

I am estimating demographic histories from RAD-seq data using Stairway Plot. One ambiguity concerns how to define the sequence length L when using one SNP per RAD locus (--write-single-snp).

I see two possible definitions for a defined population:

1. L = number of retained RAD loci × mean fragment length (i.e., sequence length represented by loci contributing to the SFS, effectively excluding monomorphic sites after SFS construction)

2. L = total callable sites reported by Stacks populations after filtering (--write-single-snp), i.e. the “all sites” value (including monomorphic sites; essentially equivalent to all retained loci for a certain population × mean fragment length)

In my dataset, option (2) gives a substantially larger L (~4M vs ~2M) because SNP thinning removes redundant SNPs but does not reduce the callable sequence length of the retained loci. 

In this situation, would you recommend using the callable-site estimate from Stacks (“all sites” / “kept loci” (variant and invariant) × mean fragment length) for L, or scaling L according to the fraction of SNPs retained after SNP thinning and SFS construction?

Best regards,
Clemens