How many SNPs vs how many loci?

[Brenna Levine]

Hi Everyone,

I am new to Stacks, and have what I'm sure are some fairly basic questions, but which really have me frazzled!

I'm trying to compare the results of different runs with different parameters, and I want to determine the total numbers of loci and snps in my data under each run. Would a line count of "batch1.catalog.tags.tsv" give me my total number of ddRAD loci, and a line count of "batch1.catalog.snps.tsv" give me my total number of snps? I have also performed the "populations" analysis and am a bit confused on how the results in the sumstats_summary.tsv file jive with the output in these other files.

For instance, if I do line counts as above, I find that there are 297,112 lines in the tags.tsv file and 24,593 lines in the snps.tsv file. However, the number of variant sites listed in the sumstats_summary.tsv file is only 16551.

For my populations run, I do not specify -m or -r, and all of my samples are from the same population. I've been through the manual, journal articles, and this online chat room for answers, but I'm still confused. Could someone  please provide me some insight into these discrepancies?

Many thanks ahead of time, and please let me know if any of you need more information!

Brenna

[Mark Ravinet]

Hi Brenna,

Running:

zcat batch1.catalog.tags.tsv | wc -l

Will give you the number of ddRADs that Stacks has been able to identify. Similarly, doing the same for the SNPs file (-1 because the first line is a header) will give you the number of SNPs identified. 

The reason these figures do not match the populations summary output is that the populations module applies filters to the dataset, as you correctly guessed. Some loci will be excluded for multiple reasons. Firstly, many of the RAD loci identified using Stacks will be monomorphic and so typically the number of loci included in the populations analysis is much lower. Some loci may also only be found in a single individual - these too should be dropped. Some loci may also be incompatible - i.e. more than two alleles present. 

In later versions of Stacks, the pipeline outputs at populations.log file. If you check this and the output to the Std Error, you will see more information on how/why loci are dropped. In principle though, this is exactly what you want - only a subset of your final dataset will be worth analysing.

Hope that helps

Mark

[Brenna Levine]

Hi Mark,

Thank you so much for your thorough response! This cleared a lot of confusion up for me.

I have two additional questions: 

(1) Could you (or someone else) explain what reasons would exclude loci from the populations module? As a reminder, I currently have neither the -m nor -r parameters set for this module.

(2) Is there a good rule of thumb for setting the -r parameter (i.e., the minimum percentage of individuals required to process a locus) in the populations module? If it helps, these loci are for a parentage/kinship study.

Thanks again for your help!

Brenna

--
Stacks website: http://creskolab.uoregon.edu/stacks/
---
You received this message because you are subscribed to a topic in the Google Groups "Stacks" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stacks-users/rhZjMOfSs1A/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stacks-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Mark Ravinet]

Hi Brenna,

I do not know exactly why but I guess a large number of RAD loci are monomorphic - i.e. no polymorphism. They will be automatically dropped because they are meaningless for analysis conducted in the populations module. If you are performing denovo assembly, some loci will only be found in a single individual and these too will be dropped because they provide no information on polymorphism. Finally some loci are dropped because they are incompatible - i.e. more than two alleles per individual. 

As for the -r parameter, I would run populations multiple times with different values of -r to see how many loci you retain. I don't have a great deal of experience with kinship studies but I would think that you certainly wouldn't want anything less than -r 0.5. 

Hope that helps

Mark

[Brenna Levine]

Hi Mark,

Thanks for your answer. Much appreciated!

Brenna

[Ross]

Sorry to interrupt in this post..

May I know, if there is any output files that indicate number of loci obtained from single individual? 

[Julian Catchen]

The number of loci assembled will be printed on to the command line by
ustacks or pstacks (which processes one individual at a time). If you
are running one of the wrappers, the number of loci assembled will be
captured in the denovo_map.log or ref_map.log files. Finally, you can
easily count the number of loci from the *.tags.tsv file for each
individual using some UNIX, or, you can look at the web interface, which
will report the number of loci and the number of polymorphic loci for
each sample.

[Ross]

TQ julian