Presented here are scripts used in data processing and analysis for "Conservation genomics for threatened New Zealand Gentianella calcis (Gentianaceae) and implications for vulnerable limestone ecosystems."
As per the suggestion of Paris et al. 2017 we performed intial parameter trials in Stacks v2.64 using a representative subset of 12 samples.
In our parameter trials, the number of loci continued to increase up to M 6 in our trials, necessitating an alternative approach to optimise parameters. We used the method of Mastretta-Yanes et al. 2014 and performed additional trials using available replicate pairs of the previously used subset.
We estimated SNP error (defined as the number of differences in SNPs between sample replicate pairs divided by the total number of SNPs) using modified code from Mastretta-Yanes et al. to compare alongside number of loci produced for each parameter combination in Stacks to see which produced the greatest number of -R 80 loci while also reducing SNP error.
As paralogs may represent an additional source of error in RAD datasets, we used HDplot to SNPs with the expected proportion of heterozygous individuals (H) greater than 5 and a “read-ratio deviation” (D) less than -4 or greater than 4 as putative paralogs. Any RAD locus containing one of these SNPs was exlcuded from further analysis.
We used various R packages including adegenet v2.1.10, ggplot2, StAMMP v1.6.3, ade4 v.1.7-20, dartR v2.9.7, and vcfR v1.13.0 to import, manipulate, analyise and visualise VCF files produced from Stacks.
We used pcadapt v4.4.3 to identify putatively neutral and outlier SNPs in a version of our dataset. These were separated into their own genlight objects using dartR and separate analysis of each undertaken to compare their patterns.