This is a test version of Biostars. For the public version, visit https://www.biostars.org.
PopGenome software: error calculating neutrality_stats

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.

snp genome popgenome r
• 5,389 views • link
updated by Ram Ram • written by anita.b0321 • anita.b0321

2 answers

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

• 0 views • link
updated by Ram Ram • written by bastianxpfeifer • bastianxpfeifer

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.

• 0 views • link
updated by Ram Ram • written by anita.b0321 • anita.b0321

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

• 0 views • link
updated by Ram Ram • written by bastianxpfeifer • bastianxpfeifer

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? lipc20052006@163.com

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
• 0 views • link
updated by Ram Ram • written by lipc20052006 • lipc20052006

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

• 0 views • link
updated by Ram Ram • written by bastianxpfeifer • bastianxpfeifer

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

• 0 views • link
updated by Ram Ram • written by anita.b0321 • anita.b0321

Log in to answer this question.

More posts like this