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:

  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.

  • $\begingroup$ 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... $\endgroup$ – llrs Feb 21 '18 at 13:46
  • 1
    $\begingroup$ 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. $\endgroup$ – Peter Pearman Feb 21 '18 at 16:52
  • $\begingroup$ I don't know the answer but have a look at Bioconductor packages (R) and Biopython (python). $\endgroup$ – llrs Feb 21 '18 at 21:58
  • $\begingroup$ 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.. $\endgroup$ – Michael Jul 15 at 20:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.