<p>I want to use SNPs produced by Ipyrad, which is a python script for RADseq, using maize genome data via <a href="https://sco.h-its.org/exelixis/web/software/raxml/index.html#documentation" rel="nofollow noreferrer" title="A maximum likelihood approach to tree building">RAxML</a> to examine the monophyly of a highly-variable focal species and its phylogenetic relationships with about 40 congeners.</p>
<p>Essentially phylogenomics using maximum likelihood trees.</p>
<p>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.</p>
<p>For subsequent use in RAxML, I would like to implement two filters to construct alignments with subsets of SNPs:</p>
<ol>
<li>SNPs that occur in at least one outgroup sample and in some (potentially maximal) number of other samples,</li>
<li>SNPs constituting a random sample from an arbitrary subset of samples, i.e. potentially excluding the outgroup samples or other samples,</li>
<li>merge the alignments into one</li>
</ol>
<p>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.</p>

CC BY-SA 4.0
2020-07-15T20:12:18.047Z
Its easy, just Biopython's AlignIO, all manipulations can be done therein. RADseq is restriction digest genome sequencing, bits of random contigs. I think the question is an example of maxium obstrufication, its quite simple really..
M__
CC BY-SA 3.0
2018-02-21T13:46:44.657Z
Sorrey but I was lost at the first sentence, could you add a link to what ipyrad, GBS and RAxML is? Also, could you explain what have you tried ? (It is better to show the effort done before asking random people in internet) You seem to have already a plan, so you might have already coded something...
llrs
CC BY-SA 3.0
2018-02-21T16:52:28.233Z
My bad.  Links are now available. The initial analysis wasn't appropriate and I didn't code it. I have no code yet, or I would have posted it. I do have some idea about what I want the code to accomplish. Since I haven't post-processed an alignment before, I aimed my question towards a non-random someone with knowledge of ways to post-process alignments (and specifically from ipyrad).  I naively tried a tibble (in R), with SNPs in columns, but it hung--must not be an appropriate approach with this much data. If I could work on rows and cols, it would be easy, but another approach must exist.
Peter Pearman
CC BY-SA 3.0
2018-02-21T21:58:28.193Z
I don't know the answer but have a look at Bioconductor packages (R) and Biopython (python).
llrs
<r><python><vcf><snp><radseq>
Post ipyrad filtering of SNP loci from GBS, based on SNP quality
# Post ipyrad filtering of SNP loci from GBS, based on SNP quality

(Asked by: username_0 on 2018-02-20)

I want to use SNPs produced by Ipyrad, which is a python script for RADseq, using maize genome data via [RAxML](https://sco.h-its.org/exelixis/web/software/raxml/index.html#documentation "A maximum likelihood approach to tree building") 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:

1. SNPs that occur in at least one outgroup sample and in some (potentially maximal) number of other samples,
2. SNPs constituting a random sample from an arbitrary subset of samples, i.e. potentially excluding the outgroup samples or other samples,
3. 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.

**Tags:** r, python, vcf, snp, radseq
> **Comments:**
> **username_2**: Sorrey but I was lost at the first sentence, could you add a link to what ipyrad, GBS and RAxML is? Also, could you explain what have you tried ? (It is better to show the effort done before asking random people in internet) You seem to have already a plan, so you might have already coded something...  
> **username_0**: My bad.  Links are now available. The initial analysis wasn't appropriate and I didn't code it. I have no code yet, or I would have posted it. I do have some idea about what I want the code to accomplish. Since I haven't post-processed an alignment before, I aimed my question towards a non-random someone with knowledge of ways to post-process alignments (and specifically from ipyrad).  I naively tried a tibble (in R), with SNPs in columns, but it hung--must not be an appropriate approach with this much data. If I could work on rows and cols, it would be easy, but another approach must exist.  
> **username_2**: I don't know the answer but have a look at Bioconductor packages (R) and Biopython (python).  
> **username_1**: Its easy, just Biopython's AlignIO, all manipulations can be done therein. RADseq is restriction digest genome sequencing, bits of random contigs. I think the question is an example of maxium obstrufication, its quite simple really..  

---

## 0 answers

