Genotype calling for an autotetraploid species without reference genome

[jorisbertrand]

Dear Lindsay,
first of all we would like to thank you for the promising initiative that PolyRAD offers for population geneticists wanting to deal with natural populations of non-model polyploid species.
We would like to use PolyRAD to call genotypes from an autotetrapoloid species (106 individuals genotyped through a nGBS protocol from 9 localities).
SNPs were called with STACKS (v.2.41) in de novo mode (as we do not have any reference genome available) for this species. As far as we know, the populations display a significant level of genetic structure and moderate to low degrees of genetic diversity.
PolyRAD seems very complete and what we want to do is rather quite simple: input data from STACKS, infer the most likely genotypes and output the results to feed other population genomics softwares.
The appropriate plan seems to follow what you showed in the video tutorial: i) quality check and filtering, ii) use the IteratePopStruct() function to infer genotypes and iii) export the results as a .vcf (or a structure file).
If we are correct, we may want to have a quick feedback on the values we obtained at different steps but before that, we would like to know whether or not it is possible to output RADdata objects (that have been called with the readStacks() function) as a vcf. We did some tests and got the following error message:
Error in RADdata2VCF(mydataPopStruct, "test.vcf") :
Complete haplotype information not provided; unable to determine SNP positions. Use refgenome argument in VCF2RADdata.
We also tried to Export_Structure(). It works but we are not sure to understand the values of the genotypes as displayed in the file.
All the best,

[lvclark]

Hi Joris,

I'm glad that you're finding polyRAD helpful! VCFs have columns for chromosome and position, so if that information is not available, polyRAD will not export to VCF. As I recall with the Stacks import option, the full tag sequences are not imported, just variable nucleotides, so that's something to keep in mind as well when interpreting allele names.

The Structure file should have four rows for each individual, and one column for each locus. Each allele is assigned a number, and the most probable genotype is shown. So if you have 1,1,2,3, that means two copies of one allele, one copy of a different allele, and then one copy of a third allele.

Best,

Lindsay

[jorisbertrand]

Hi Lindsay,
Yes! adegenet would definitely be an option to perform analyses such as PCA or DAPC. Also being able to import the data in Genodive could be great as this software is able to perform polyploid dosage corrections. I guess including such features in PolyRAD could allow to make it more broadly used in our field.
Just let me know if you need some more information about our data.
Thank you for your precious help,
Cheers
Joris

[lvclark]

FYI, polyRAD does PCA, and the results are available in the $PCA slot. One could also use GetWeightedMeanGenotypes to generate a matrix that could be passed to the adegenet::find.clusters and adegenet::dapc functions, bypassing the need to import as genind.

[lvclark]

And if importing from polyRAD, the procedure for correction of allele dosages in polyploids in GenoDive would be unneccessary, since polyRAD has already estimated allele dosages using a similar procedure. Supposedly GenoDive can import from Structure format, although I'm not sure if it would require one to add a column indicating population assignments. If it doesn't work, let me know. A GenoDive export function would be more involved for me to create than a genind export function would be, and I'm working on a grant this week so I don't want to make more work for myself than is necessary.

[lvclark]

Dear Joris,

I just added a Export_adegenet_genind function to the version of polyRAD on GitHub. If you get a chance to test it out this week, let me know how it goes! I need to update polyRAD on CRAN very soon anyway, so if there are any tweaks needed for the new function I might as well include them.

Best,

Lindsay

[jorisbertrand]

Thanks a lot! we are gonna try. Also we managed to open the structure file in GenoDive. By default, it was automatically trying to open it as a 'Migrate IAM' file, so we add to import it as a text file, change the format to 'Structure' and then use the menu Data > Transformation to set up the ploidy level. Keep in touch!