PopGenome software: error calculating neutrality_stats

Question

[anita.b0321]

Dear all,

I have a question regarding the PopGenome package for R which I have just begun to get comfortable with.

I have been loading my VCF files individually by chromosome and calculating stats on them. I've done it for every other chromosome except chromosome 8. All files have been processed the same exact way. However, when I use the neutrality_stats function I get the following error only with chromosome 8:

|Error in if (outgroup[1] == FALSE) { : argument is of length zero

This is weird as it only happens with this one chromosome. Below are the commands I use:

GENOME.class_gla<-readVCF("chr8.vcf.gz",1000,"8",1,146364022,gffpath = "chr8.gff3")
genes_gla<-splitting.data(GENOME.class_gla,subsites="gene")
genes_gla<-neutrality.stats(genes_gla,FAST=TRUE)

How can I solve this? If anyone has any suggestions, I would be so grateful.

Answer 1

[bastianxpfeifer]

Hi,

rm(list=ls())
gc()

quit R (q()) and save the empty workspace.

Note, you will delete all objects within your R session with these commands.

Try again those PopGenome calls. I guess the error will disappear.

Best,
Bastian

[anita.b0321]

Hi Bastian,

Thanks for your suggestion. It did not solve the problem however. I really don't know the cause, but I will start by simply recreating the VCF which was a subset of the original larger file.

edit Recreating the VCF did nothing to help also. Still getting the same error.

[bastianxpfeifer]

Thats weird, indeed. FAST=FALSE in the neutrality.stats() module does not fix the problem ?

Ok, probably the following is a fast and dirty fix:

GENOME.class_gla<-readVCF("chr8.vcf.gz",1000,"8",1,146364022,gffpath = "chr8.gff3")
GENOME.class_gla <- set.outgroup(GENOME.class_gla, FALSE)
genes_gla<-splitting.data(GENOME.class_gla,subsites="gene")
genes_gla<-neutrality.stats(genes_gla, FAST=TRUE)

All the best,
Bastian

[lipc20052006]

Hi Bastian

I'm trying to read a HapMap file with readData, but I get an error,

> GENOME.class<-readData("D:/R/a/b", format="HapMap")
|            :            |            :            | 100 %
|
Read 5764 items
Error in CCC$matrix : $ operator is invalid for atomic vectors

These are the first ten lines and 16 columns of the hapmap file

Can you help me? [e-mail]

rs#  alleles  chrom  pos    strand  assembly  center  protLSID  assayLSID  panelLSID  QCcode  L1  L2  L3  L4  L5
 S1   A/G      1      65513  +       NA        NA      NA        NA         NA         NA      AA  AA  AA  GG  AA
 S2   G/A      1      86025  +       NA        NA      NA        NA         NA         NA      AA  GG  GG  GG  GG
 S3   A/G      1      98468  +       NA        NA      NA        NA         NA         NA      AA  AA  AA  AA  AA
 S4   G/A      1      2E+05  +       NA        NA      NA        NA         NA         NA      GG  GG  GG  GG  GG
 S5   G/A      1      2E+05  +       NA        NA      NA        NA         NA         NA      GG  GG  GG  GG  GG
 S6   G/A      1      2E+05  +       NA        NA      NA        NA         NA         NA      GG  GG  GG  GG  GG
 S7   G/A      1      2E+05  +       NA        NA      NA        NA         NA         NA      GG  AA  AA  GG  AA
 S8   G/A      1      2E+05  +       NA        NA      NA        NA         NA         NA      GG  GG  GG  GG  GG
 S9   G/A      1      2E+05  +       NA        NA      NA        NA         NA         NA      GG  AA  AA  GG  AA
 S10  A/C      1      2E+05  +       NA        NA      NA        NA         NA         NA      AA  AA  AA  AA  AA

Answer 2

[bastianxpfeifer]

Thats weird, indeed. FAST=FALSE in the neutrality.stats() module does not fix the problem?

Ok, probably the following is a fast and dirty fix:

GENOME.class_gla<-readVCF("chr8.vcf.gz",1000,"8",1,146364022,gffpath = "chr8.gff3")
GENOME.class_gla <- set.outgroup(GENOME.class_gla, FALSE)
genes_gla<-splitting.data(GENOME.class_gla,subsites="gene")
genes_gla<-neutrality.stats(genes_gla, FAST=TRUE)

All the best,
Bastian

[anita.b0321]

Hi Bastian,

Setting FAST=FALSE did not fix the problem. And unfortunately neither did the setting the outgroup to FALSE. However, I did manage to get it working by adding the command include.unknown=TRUE when loading the VCF.

Thanks for your help. This really is an awesome tool.

-Anita