A lot a SNPs or wrong analysis?

[Yesenia Vega]

Hello everyone. 

Im new with stacks and Im trying to get a de novo assamble.

I have ddRAD 2x110 pb. I tried to set the parameters according to Paris et al.  2017 and Rochette et al 2017. The best parameters for r80 were m4 M6 and n7. 

The final data set is for three species with a total of 91 samples in 26 populations (2-5 samples per population) and with those parameters my results were:

Coverage was 22.2X (8.6/50.6). Genotyped loci: 121 114.

After populations (-r 0.5 -p 7 --write-single-snp --min-maf 0.05 --max-obs-het 0.6) I get:

Assembled loci: 26 309, variants remained: 25 815.

To get the number of polymorphic loci I used: cat populations.hapstats.tsv | grep -v "^#" | wc -l I got: 334 580

To get the number of SNPs I used: cat populations.sumstats.tsv | grep -v "^#" | wc -l I get:  330 543

So my questions are, is it right the way am I getting the number of SNPs? Why the number is so different to the "variants remained"?

Thank you a lot in advance,

[seanc...@gmail.com]

Aloha Yesenia,

I just came across this issue tonight, and the reason is that the file populations.sumstats.tsv lists variants per population. That means that they're potentially listed multiple times, and your command is counting the same SNP/site several times over.

That command (cat populations.sumstats.tsv | grep -v "^#" | wc -l) was originally meant to be used when you group all individuals into a single population, as you would when testing r80 parameters.

Hope this helps!

- Sean

[Aron Katz]

Hi Yesenia, 

How were you able to retrieve coverage information for the r80 loci retained after running populations?

Best,

Aron