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

[Ryan Gutenkunst]

Hello Clemens,

I’ve not used Stacks before, but based on another user’s reports… I would use the all sites length, then multiply it by the ratio of SNPs after thinning to before thinning. I assume that “all sites” is the total number of base pairs for which calls could be made, and you’re then artificially reducing the SNP count by the ratio of after to before thinning.

Best,

Ryan

--
You received this message because you are subscribed to the Google Groups "dadi-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dadi-user+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dadi-user/83a1adca-d94c-4dfd-8054-b57c7a62f812n%40googlegroups.com.