I want to use SNPs produced by Ipyrad, which is a python script for RADseq, using maize genome data via RAxML to examine the monophyly of a highly-variable focal species and its phylogenetic relationships with about 40 congeners.
Essentially phylogenomics using maximum likelihood trees.
An initial analysis suggested that samples from two outgroup species might share few SNPs with other samples, so the threshold for shared loci was set to two samples. As a result, the .u.snps.phy alignment has over 600,000 SNPs.
For subsequent use in RAxML, I would like to implement two filters to construct alignments with subsets of SNPs:
- SNPs that occur in at least one outgroup sample and in some (potentially maximal) number of other samples,
- SNPs constituting a random sample from an arbitrary subset of samples, i.e. potentially excluding the outgroup samples or other samples,
- merge the alignments into one
How can I do this efficiently, maximizing use of existing functionality in R or Python? Either code or a (likely successful) strategy would be helpful. I have access to all output formats from ipyrad as potential starting points.