lvclark / polyrad

Genotype Calling with Uncertainty from Sequencing Data in Polyploids 🍌🍓🥔🍠🥝

polyRAD: Genotype Calling with Uncertainty from Sequencing Data in Polyploids and Diploids

The ultimate goal of polyRAD is to enable research that strives to end world hunger and conserve biodiversity. On that note, I actively welcome input and contributions from diverse and underrepresented sources to help improve the software, and have adopted the Contributor Covenant Code of Conduct. If there is something I can do or add to make polyRAD more useful and approachable to you, post an issue, start a discussion or send me an email! If you would like to create a new function yourself or propose edits to the code or documentation, feel free to make a pull request! All levels of experience are welcome and you can contact me if you aren't sure how to get started with GitHub.

See the list of future features and issues for some of the things that I'm already thinking of adding to polyRAD.

polyRAD is available on CRAN or can be installed from this repository.

Table of Contents

• Purpose
• Why polyRAD
• Formats supported
• Installation
• Tutorial
• Citation
• Need help with R?
• Funding

Purpose

Genotypes derived from genotyping-by-sequencing (GBS) and restriction site-associated DNA sequencing (RAD-seq) have inherent uncertainty associated with them due to sampling error, i.e. some alleles might not get sequenced at all, or might not be sequenced in exact proportion to their copy number in the genome. This package imports read depth in a variety of formats output by various bioinformatics pipelines and estimates the probability of each possible genotype for each taxon and locus. Unlike similar pipelines, polyRAD can account for population structure and variable inheritance modes (autopolyploid, allopolyploid, intermediate). Genotypes and/or probability distributions can then be exported for downstream analysis such as genome-wide association, genomic selection, QTL mapping, or population structure analysis.

Starting with version 1.2, polyRAD also includes its own variant calling pipeline. It is designed for highly duplicated (typically allopolyploid) reference genomes, and can help correct and filter markers that appear to consist of multiple paralogous loci. Even in non-duplicated genomes, the same statistic can be used for filtering non-Mendelian loci, as well as loci with considerable amounts of amplification bias and/or overdispersion in the sequencing data. This filtering can be performed before genotype calling.

Why polyRAD?

If you're like me, you don't want to waste a lot of money sequencing your DNA samples at a higher depth than is necessary. You would rather spend that money adding more samples to the project, or using a different restriction enzyme to get more markers! You may have also noticed that some loci get sequenced at a much higher depth than others, which means that if you sequence the same library a second time, you aren't likely to get a lot of reads for the loci that need it most. So how can we get the maximum amount of information out of sequencing data where many loci are low depth? And, for example, if we only have five reads, how can we estimate allele dosage in a heterozygous octoploid?

The answer that polyRAD provides is a Bayesian genotype caller with many options for specifying genotype prior probabilties. When read depth is low, accurate priors make a big difference in the accuracy of genotype calls. And because some genotype calls are going to be uncertain no matter how sophisticated our algorithm is, polyRAD can export genotypes as continuous numeric variables reflecting the probabilities of all possible allele copy numbers. This includes genotypes with zero reads, where the priors themselves are used for imputation.

For an explanation of Bayesian genotype calling and posterior mean genotypes, see a seminar of mine.

Genotype priors in diversity panels and natural populations:

• Either assume no population structure (HWE), or let polyRAD infer population structure and model allele frequency gradients.
• The user specifies a rate of self-fertilization ranging anywhere from zero to one.
• If loci have known positions in a reference genome, polyRAD can search for loci in linkage disequilibrium and use those loci to update priors.

Genotype priors in biparental mapping populations:

• The user specifies the number of generations of backcrossing, intermating, and/or selfing. (For an F1 population, all three would be zero.)
• Based on likely parental genotypes and allele frequencies in the progeny, polyRAD determines the segregation pattern of each marker.
• If the loci have known positions in a reference genome or on a map, linked markers can be used for updating priors.

In particular, by using population structure and linkage to inform genotype priors on a per-individual basis, high depth markers are used by polyRAD to improve the accuracy of genotyping at low depth markers. All pipelines allow autopolyploidy, allopolyploidy, or some mixture of the two. And because non-model organisms need some love, reference genomes are optional.

Formats supported

To hopefully answer the question, "Can I use polyRAD?":

polyRAD requires as input the sequence read depth at each allele for each sample. Alleles must also be grouped into loci. The bioinformatics pipeline that you used for SNP discovery did not have to assume polyploidy, as long as it faithfully reported allelic read depth. Genomic alignment information is optional. Right now there are data import functions for the following formats:

• Variant Call Format (VCF). The allele depth (AD) genotype field must be present. I have tested the import function on files produced by the TASSEL GBSv2 pipeline. It should also work for GATK and Tassel4-Poly.
• TagDigger. This is another piece of software I created, which reads FASTQ files and expected tag sequences and outputs a CSV file with read counts at samples x tags.
• UNEAK. The UNEAK pipeline outputs read depth in a file called HapMap.hmc.txt, which can be read by polyRAD. (Beware that read depth is capped at 127 by UNEAK; TagDigger can help you if you expect high depth to be common in your dataset.)
• Stacks. If you have catalog files (catalog.alleles.tsv etc.) and matches files (matches.tsv) generated by Stacks, they can be imported by polyRAD.
• TASSEL-GBSv2. Rather than running the entire pipeline, you can run GBSSeqToTagDBPlugin, TagExportToFastqPlugin, Bowtie2 or BWA as described in "Run Alignment Program(s)", then GetTagTaxaDistFromDBPlugin to get the TagTaxaDist and SAM files needed for import by polyRAD.
• DArTag. This uses a custom format developed in conjunction with Breeding Insight.

Currently there are export functions for the following software. Genotypes are exported as continuous variables (posterior mean genotypes) for these four formats. There are also functions to generate matrices of continuous or discrete genotypes, which can be used in custom export functions.

• GAPIT and FarmCPU
• rrBLUP
• TASSEL
• GWASpoly

There are export functions for discrete genotypes for the following software:

• polymapR
• Structure
• adegenet

Discrete genotypes can also be exported to VCF for import to a variety of software.

Genotype probabilities can be exported to:

• MAPpoly
• polymapR

Installation

polyRAD depends on some Bioconductor packages. Before attempting to install polyRAD, run

if (!requireNamespace("BiocManager", quietly = TRUE)){
  install.packages("BiocManager")
}

BiocManager::install("pcaMethods")

If you plan to import from or export to VCF, also run

BiocManager::install("VariantAnnotation")

polyRAD can then be installed from CRAN with

install.packages("polyRAD")

Alternatively, if there are new features not yet on the CRAN version that you want to use, you can install the development version here on GitHub at your own risk using the command

install.packages("polyRAD", repos = "https://lvclark.r-universe.dev")

Tutorial

The tutorial document for the package is available on R-universe.

A separate tutorial is available for the variant calling pipeline.

A video tutorial from January 2021 is available on YouTube.

A tutorial on using polyRAD for population genetics is also available. If you would like to have the dataset to follow along with that tutorial, run install.packages("polyRADtutorials", repos = "https://lvclark.r-universe.dev"). The tutorial covers genotype calling in a wild collection of diploid, triploid, and tetraploid Miscanthus sacchariflorus (a perennial grass) and some downstream applications in R.

The population genetics tutorial is also available as a video from April 2022. English closed captions are available for the video.

Note that the tutorials hosted on R-universe are derived from the development version of polyRAD. If you installed from CRAN and something from an R-universe tutorial isn't working, check the tutorial that came with your installation using browseVignettes("polyRAD").

Citation

If you use polyRAD, please cite this manuscript:

Clark LV, Lipka AE, and Sacks EJ (2019) polyRAD: Genotype calling with uncertainty from sequencing data in polyploids and diploids. G3 9(3) 663--673. doi:10.1534/g3.118.200913

The Hind/He statistic is described in the following manuscript:

Clark LV, Mays W, Lipka AE, and Sacks EJ (2022) A population-level statistic for assessing Mendelian behavior of genotyping-by-sequencing data from highly duplicated genomes. BMC Bioinformatics 23: 101, doi:10.1186/s12859-022-04635-9.

The variant calling pipeline is described in version one of the preprint:

Clark LV, Mays W, Lipka AE, and Sacks EJ (2020) A population-level statistic for assessing Mendelian behavior of genotyping-by-sequencing data from highly duplicated genomes. bioRxiv doi:10.1101/2020.01.11.902890

Additional citations

Citable Zenodo DOI for the software:

polyRAD has also been presented in these posters:

Clark LV, Lipka AE, and Sacks EJ (2019) Improvements to Genotype Calling in Polyploids Using the polyRAD R Package. Plant and Animal Genome Conference XXVII, January 12-16, San Diego, California, USA. doi:10.13140/RG.2.2.18358.75847

Clark LV, Lipka AE, and Sacks EJ (2018) polyRAD: Genotype Calling with Uncertainty from Sequencing Data in Polyploids and Diploids. Plant and Animal Genome Conference XXVI, January 13-17, San Diego, California, USA. doi:10.13140/RG.2.2.27134.08001

Need help with R?

Video tutorials on R, Python, and linear algebra were produced as part of the NSF grant that funded polyRAD. They can be found here: https://www.youtube.com/channel/UCcyeu245nwXnUDoxvcBPqVg

Funding

This material is based upon work supported by the National Science Foundation under Grant No. 1661490.
Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

polyrad's Issues

Support pedigreed populations

I'd like to make a new pipeline that uses pedigree information. Genotype estimates of parents and offspring will iteratively influence genotype priors of parents and offspring. Even for biparental populations, this could perform much better than the existing pipeline, which doesn't handle segregation distortion well.

I probably won't tackle this until I have added support for multiploid populations (Issue #17), in order to avoid having to rewrite a lot of code after the fact.

If you have a good test dataset for this sort of population, please contact me!

"ExpectedHindHe" not found while all other functions of polyRAD work

Hello
I have followed the tutorial and succeeded in running the command lines but I have got a error message for this step:

"ExpectedHindHe(mydata, inbreeding = 0.94, ploidy = 2, reps = 10)
Error in ExpectedHindHe(mydata, inbreeding = 0.94, ploidy = 2, reps = 10) : impossible to find the function"

I have tried several times on several Rstudio. Could you please help me ?
Thanks!
my best
Amandine

Support multiploid populations

For certain crops like banana and yam there seems to be a need for genotyping within populations of mixed ploidy. I have started some work on this, but it will involve some restructuring of RADdata objects, enough to probably justify incrementing the major version number of the software (going from polyRAD 1.3 to polyRAD 2.0, since objects will not be backwards-compatible).

I plan to support cases where ploidy is known ahead of time e.g. by flow cytometry, but not ploidy estimation from the data. The Hind/He statistic can help identify accessions where ploidy was misidentified, but keep in mind that it is also highly influenced by inbreeding and hybridization. See https://github.com/delomast/tripsAndDipR for a package that estimates ploidy from sequencing data.

The possiblePloidies slot will continue to represent different inheritance modes across the genome, so different markers can have different inheritance models. I plan to add a taxaPloidy slot, which will be an integer vector of individual ploidies, used as a multiplier on top of possiblePloidies.

There will need to be some substantial changes to the biparental mapping population pipeline, in particular how ploidies of the parents and progeny are specified. The priorProbPloidies slot can possibly be eliminated.

Export to additional population genetics formats

Overview

There has been some interest in using genotypes called by polyRAD for population genetics analysis. Two pieces of software that have long supported polyploid genotypes are Structure and SPAGeDi. Although they were originally designed for relatively small sets of PCR markers, they both scale up to large sequence-based datasets. They each use a custom format, and it would be nice if polyRAD exported to those formats.

I may eventually get to writing these export functions. If you would like to see them (or others), feel free to comment here! If you have some R experience and would like to write the functions yourself, see below. First-time contributors are welcome!

Contributing

You will need to make your own fork of polyRAD, clone that fork to your computer, commit your changes, push those commits to GitHub, and open a pull request. See the Github documentation if you need help understanding those tasks.

R code

Name the function Export_Structure or Export_SPAGeDi, as appropriate. Add it to the file R/data_export.R.

The first argument for the function should be a RADdata object named object. There should also be an argument called file indicating the path to write to. Another possible argument might indicate a particular ploidy to use for export; see other export functions, and check the documentation for Structure and SPAGeDi to see if variable ploidy is allowed.

The function will need to call GetProbableGenotypes with omit1allelePerLocus = FALSE and multiallelic = "correct". This will get you a list, the first item of which is a matrix with individuals and rows and alleles in columns, showing the copy number of that allele in that individual. This will then need to be converted to a different format, where each allele is assigned an integer, and the individual has a number of alleles up to its ploidy. For example:

	Allele_AG	Allele_AA	Allele_CA
Ind 1	0	3	1
Ind 2	2	2	0

Gets changed to

	Allele 1	Allele 2	Allele 3	Allele 4
Ind 1	2	2	2	3
Ind 2	1	1	2	2

Your function will make use of object$alleles2loc in order to group columns of the matrix output by GetProbableGenotypes into loci. A loop in R to process loci might look like:

geno <- GetProbableGenotypes(object, omit1allelePerLocus = FALSE)[[1]]

for(L in seq_len(nLoci(object)){
  theseal <- which(object$alleles2loc == L)
  thesegeno <- geno[, theseal]
  # do some processing here
}

See the MakeGTstrings internal function in src/PrepVCFexport.cpp. For SPAGeDi in particular this function could simply have an argument added to it to start numbering from 1 rather than 0. For Structure, MakeGTstrings should probably not be used directly but can provide some inspiration for what needs to be done.

See the Structure or SPAGeDi documentation for details on file format needed. Your function should export a text file in that format.

Test the function on exampleRAD after running IterateHWE. Install Structure or SPAGeDi on your computer and see if it can import the file.

Documentation and integration

• Update man/ExportGAPIT.Rd to document your function.
• Add the function name to the export call in NAMESPACE
• Update the "Formats supported" section of README.md
• Update the "Summary of available functions" section of vignettes/polyRADtutorial.Rmd.
• Run vignettes/render_vignettes.R , with vignettes as the working directory, to recompile the vignette.
• Add yourself as a contributor to the DESCRIPTION file! Use the "ctb" role code in Authors@R and Author.
• Run R CMD build and R CMD check to make sure the package builds and passes checks.

HidHe is not found

I am following tutorial for polyRad and noticed that function HindHe is not found when I wanna check skewness of markers based on Mendelian ratio. which library should be loaded?

Best practice for Multiparent population

Hi,
THanks @lvclark for developing this package.

We have a multiparent heterozygous autotetraploid population with four families, one shared pollen donor and four mums, so offspring is half siblings. Our input in polyRAD is a filtered VCF with ok data.

We are estimating genotype probabilities one family at the time following the tutorial. We use SubsetByTaxon and run one family at the time (quality control and PipelineMapping2Parents). Our PCAs per family look very similar to the example. We are exporting genotype probabilities for Mappoly and PolymapR, so we finish with four files x2 formats, one per family.
Offspring per family is uneven and we could benefit by "joining" the dataset somehow.

Questions:
=Is there a better way to work with these families than one at the time?
=We have replicates for the parents, can we set donor/recurrent parent with several samples?

THanks

Error message on R when reading isoloci

Hello Lindsay,

Thanks for the polyRAD package, it is well documented and the tutorials are awesome!

Following the steps for "Variant and Genotype Calling in Highly Duplicated Genomes" from https://lvclark.r-universe.dev/articles/polyRAD/isolocus_sorting.html, to call variants on my own dataset on a Miscanthus F2 diploid population, I get to the steps before importing the sorted dataset into polyRAD without any problems, but when trying to import the sorted .csv file generated by process_isoloci.py in R to polyRAD with the function readProcessIsoloci, I get the following error:

myRAD <- readProcessIsoloci("20230315_out_sam_multi_1_sorted.csv", possiblePloidies = list(2), taxaPloidy = 2)
Reading file...
Filtering and sorting loci...
Building RADdata object...
Error in sum(sapply(seq_len(nsites), function(i) polyRADsubmat[splitnuc1[[a1]][i], :
invalid 'type' (list) of argument

I don't quite understand the error message. Could you help me figure out what am I doing wrong? I tried changing default arguments, but still gives me the same error message no matter what I do.

Thanks,

Export to and import from updog

It would be nice to have some convenience functions to convert between a RADdata object and the input and output of updog. This would allow users to take advantage of the file import and export options in polyRAD, while performing the genotype calling itself in updog (more accurate than polyRAD in some cases but much slower).

If you would like to add this feature and make a pull request, just comment here and I will give any help and guidance that I can. In particular see the multidog and format_multidog functions in updog. See also the checklist for pull requests.

How to use user input VCF file?

Hi there,

I had an error when I was trying to call VCF2RADdata, the command is:
VCF2RADdata(myVCF, possiblePloidies = list(c(2,2)), expectedLoci = 100, expectedAlleles = 500)
the error is:
"Error in VCF2RADdata(myVCF, possiblePloidies = list(c(2, 2)), expectedLoci = 100, :
Need to adjust min.ind.with.reads for number of samples in dataset."

Could you please tell me how to fix this?
Besides, does your method work for one single individual? Thanks.

Best,
Yudi

RADdata2VCF with asSNPs=TRUE gives error

Trying to export VCF file after running genotyping and filtering on my data. No reference was provided when using VCF2RADdata to import the original data. Here's my function call to export the data:

RADdata2VCF(mydata_75, file = "polyrad_out_75.vcf", asSNPs=TRUE, hindhe=TRUE)

Here's the error I get:

  Complete haplotype information not provided; unable to determine SNP positions.  Use refgenome argument in VCF2RADdata.

I can successfully run the function and output the VCF if I wrap the code giving the error [e.g. if(is.null(varok) || varok)] with a "if (asSNPs==FALSE){}" check.

Not sure if this is the desired behavior or if it is an error on my part?

Thanks,

Tyler

TestOverdispersion: Genotype quality too low for this protocol

Hello,

I am having issues with this function. What other protocol do you suggest if the genotypes are too low? Attached is the structure of my RAD object.

Error running Hind/He script in R

Hello Lindsay,

First of all thank you for your effort putting all these together.

I study a group of aquatics Eriocaulon, and just recently we found out we might have a mix of diploid and polyploid depending the population and one specific case where I might have one pop with different ploidies.

I have been reading your recent paper on Hind/He and I tried running the script. Unfortunately I have been unable to run it, as a constant error pops up. I've ran the following command and got the following error:

myRAD <- VCF2RADdata("populations.haps.vcf",

•                  svparam = ScanVcfParam(fixed = "ALT", info = NA, geno = "AD",
  

•                                         which = GRanges("06", IRanges(1, 937644))),
  

•                  yieldSize = NA)
  

Error in ScanVcfParam(fixed = "ALT", info = NA, geno = "AD", which = GRanges("06", :
could not find function "ScanVcfParam"

I've looked at this page for answers but I couldn't work it out. https://support.bioconductor.org/p/9137213/

I'm doing a De novo analysis on Stacks program, and I understood I needed to use the haplotype.vcf file from Stacks' folder, after a De novo analysis is finished. If there is anything you could tell me or point me in the right direction, I very much appreciate it.

Thanks very much for your time.

number of loci

Hi, I am using polyRAD to generate allele dosage in alfalfa. I used NGSEP (https://sourceforge.net/p/ngsep/wiki/Home/) to generate the VCF. However the number of loci generated by VCF2RADdata is lower compared with the original VCF. Do you know the way to obtain the same number of sites? Thank you.

myRAD7 <- VCF2RADdata("logan_8.vcf.bgz",
                      refgenome = "MtrunA17r5.0-20161119-ANR.fasta",
                      al.depth.field = "ACN",
                      possiblePloidies = list(4))
> myRAD7
## RADdata object ##
272 taxa and 3807 loci
3811212 total reads
Assumed sample cross-contamination rate of 0.001

Possible ploidies:
Autotetraploid (4)

logan8 <- readVcf("logan_8.vcf")
> logan8
class: CollapsedVCF 
dim: 6862 272 
rowRanges(vcf):
  GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
info(vcf):
  DataFrame with 12 columns: CNV, TA, TID, TGN, TCO, TACH, NS, MAF, AN, AFS, TYPE, FS
info(header(vcf)):
        Number Type    Description                                                                                                 
   CNV  1      Integer Number of samples with CNVs around this variant                                                             
   TA   1      String  Variant annotation based on a gene model                                                                    
   TID  1      String  Id of the transcript related to the variant annotation                                                      
   TGN  1      String  Name of the gene related to the variant annotation                                                          
   TCO  1      Float   One based codon position of the start of the variant. The decimal is the codon position                     
   TACH 1      String  Description of the aminoacid change produced by a non-synonymous mutation. String encoded as reference am...
   NS   1      Integer Number of samples genotyped                                                                                 
   MAF  1      Float   Minor allele frequency                                                                                      
   AN   1      Integer Number of alleles in called genotypes                                                                       
   AFS  R      Integer Allele counts over the population for all alleles, including the reference                                  
   TYPE 1      String  Type of variant                                                                                             
   FS   1      Float   Phred-scaled p-value using Fisher's exact test to detect strand bias                                        
geno(vcf):
  List of length 7: GT, PL, GQ, DP, ADP, BSDP, ACN
geno(header(vcf)):
        Number Type    Description                                                                                                 
   GT   1      String  Genotype                                                                                                    
   PL   G      Integer Phred-scaled genotype likelihoods rounded to the closest integer                                            
   GQ   1      Integer Genotype quality                                                                                            
   DP   1      Integer Read depth                                                                                                  
   ADP  R      Integer Counts for observed alleles, including the reference allele                                                 
   BSDP 4      Integer Number of base calls (depth) for the 4 nucleotides in called SNVs sorted as A,C,G,T                         
   ACN  R      Integer Predicted copy number of each allele taking into account the prediction of number of copies of the region...

Best,

Cesar

Error in RADdata2VCF

Hello Lindsay,

Thanks for the polyRAD package. I am trying to use it to estimate polyploid genotypes for my supposedly polyploid RADseq data (which are samples from populations in the wild), for which there is not a reference genome available. My goal is to use the estimated polyploid genotypes for further analysis starting from a VCF file.

I imported data from STACKS using readStacks with the 'catalog', 'matches files' and the 'sumstatsFile'. I was able to follow the tutorial "Estimating genotype probabilities in a diversity panel" but I have an issue when trying to export the RADdata object to VCF format using the RADdata2VCF function. I get this error message:

Error in RADdata2VCF(myRADdata, file = "file.vcf",  : 
     Complete haplotype information not provided; unable to determine SNP positions.  Use refgenome argument in VCF2RADdata.

It seems I need the 'refgenome argument', but to my understanding I am lacking this.
Do have any idea how can I come around this?

Thanks

Export_MAPpoly error

Hello

I am using polyRAD for calling genotypes in alfalfa crop. I have F1 progenies SNP data. I am having issues exporting the data out of polyRAD to run analysis in Mappoly. Here is the error I am getting:

Error in object$posteriorProb[[ploidyIndex, as.character(pld.p)]] : subscript out of bound

Thanks for your help!

Harpreet