Error running Hind/He script in R about polyrad

[lvclark]

Hi Paulo,

Just run library(VariantAnnotation) before running your script.

[PauloBaleeiro]

Unfortunately it is still not working. It might be the GRanges and IRanges. I'm using the zipfile cataloge.vcf file generated by Stacks, but I cannot work out where the problem is. I tried running without "which=...", but it didn't work either.

[lvclark]

What error are you getting?

[PauloBaleeiro]

Hi Lindsay,

This is the beginning of script followed by error:
myRAD <- VCF2RADdata("catalog.alleles.tsv.gz",

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

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

•                  yieldSize = NA)
  

Error in h(simpleError(msg, call)) :
error in evaluating the argument 'object' in selecting a method for function 'samples': [internal] _hts_rewind() failed

[lvclark]

Oh, you can't read a Stacks file with VCF2RADdata. There is a different function called readStacks.

…

On Sat, Apr 15, 2023, 1:14 AM Paulo Baleeiro ***@***.***> wrote: Hi Lindsay, This is the beginning of script followed by error: myRAD <- VCF2RADdata("catalog.alleles.tsv.gz", - svparam = ScanVcfParam(fixed = "ALT", info = NA, geno = "AD", - which = GRanges("06", IRanges(1, 937644))), - yieldSize = NA) Error in h(simpleError(msg, call)) : error in evaluating the argument 'object' in selecting a method for function 'samples': [internal] _hts_rewind() failed — Reply to this email directly, view it on GitHub <#37 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADG72KBEXKNKWG37HTUKZS3XBIVCNANCNFSM6AAAAAAW4NNFDA> . You are receiving this because you commented.Message ID: ***@***.***>

[PauloBaleeiro]

Hi Clark,

I'm a bit stuck in the "object" to be used. In one of your scripts you show the following path:
VCF2RADdata(file, phaseSNPs = TRUE, tagsize = 80, refgenome = NULL,
tol = 0.01, al.depth.field = "AD", min.ind.with.reads = 200,
min.ind.with.minor.allele = 10, possiblePloidies = list(2),
taxaPloidy = 2L, contamRate = 0.001,
samples = VariantAnnotation::samples(VariantAnnotation::scanVcfHeader(file)),
svparam = VariantAnnotation::ScanVcfParam(fixed = "ALT", info = NA,
geno = al.depth.field,
samples = samples),
yieldSize = 5000, expectedAlleles = 5e+05, expectedLoci = 1e+05,
maxLoci = NA)

I'm not sure which VCF files to use. Would it be "populations.snps.vcf? or "populations.haps.vcf"?
And what about the catalog.alleles.tsv?
I have tried many ways, but none worked.

Cheers

[lvclark]

Hi Paulo,

That looks like the usage section of the VCF2RADdata documentation. The argument defaults are shown, and you don't have to type the arguments if you are okay with the defaults.

I haven't used Stacks in a long time. If your VCF looks like this:
https://github.com/galaxyproject/tools-iuc/blob/main/tools/stacks2/test-data/populations/populations.haps.vcf
it won't work because the AD field is not included. However one like this should work: https://github.com/galaxyproject/tools-iuc/blob/main/tools/stacks2/test-data/populations/populations.snps.vcf

Your code (in its simplest version) would look like:

mydata <- VCF2RADdata("populations.snps.vcf")

You don't need catalog.alleles.tsv if you are using VCF2RADdata. That file is only needed for readStacks.

Best,

Lindsay

[PauloBaleeiro]

This is mine. I'll try once again. Thanks

[lvclark]

Should be okay. I was concerned about the missing AD field for missing genotypes (./.) but did a quick test and polyRAD handles those appropriately. If you continue to have issues, please provide your code and a small dataset so that I can reproduce the problem on my own machine.

[lvclark]

That being said, you will have to lower the values for the min.ind.with.reads and min.ind.with.minor.allele arguments to VCF2RADdata since you have so few samples. polyRAD's genotyping algorithms may perform poorly on a sample size this small but you can use the naïve protocol shown in https://lvclark.r-universe.dev/articles/polyRADtutorials/population_genetics.html (just search the page for mydataNaive and follow that protocol for genotyping).