Sample.idError with seqVCF2GDS

[alexisregelson]

Hello, I am trying to use seqVCF2GDS and am getting the following error:

library(SeqArray)
library(data.table)

seqVCF2GDS(high_mod_vcf, "r4_chr1_high_mod.gds", parallel=6L)
Mon Nov 6 16:09:06 2023
Variant Call Format (VCF) Import:
file(s):
r4_PASS_chr1_updated_varID_dups_drop_updated_IDs_nhw_hwe6_noNHWrelateds_high_mod_impact.vcf (198.8M)
file format: VCFv4.2
the number of sets of chromosomes (ploidy): 2
the number of samples: 14,306
genotype storage: bit2
compression method: LZMA_RA
# of samples: 14306
calculating the total number of variants ...
the total number of variants for import: 3,632
Writing to 6 files:
r4_chr1_high_mod_tmp01_ad336f56fc72 [1..606]
r4_chr1_high_mod_tmp02_ad3315e862b7 [607..1,212]
r4_chr1_high_mod_tmp03_ad33613818b1 [1,213..1,818]
r4_chr1_high_mod_tmp04_ad33473817c6 [1,819..2,424]
r4_chr1_high_mod_tmp05_ad334e0fea8c [2,425..3,030]
r4_chr1_high_mod_tmp06_ad33607634f8 [3,031..3,632]
Done (Mon Nov 6 16:09:10 2023).
Output:
r4_chr1_high_mod.gds
Merging:
opening 'r4_chr1_high_mod_tmp01_ad336f56fc72' ... [done]
opening 'r4_chr1_high_mod_tmp02_ad3315e862b7' ... [done]
opening 'r4_chr1_high_mod_tmp03_ad33613818b1' ... [done]
opening 'r4_chr1_high_mod_tmp04_ad33473817c6' ... [done]
opening 'r4_chr1_high_mod_tmp05_ad334e0fea8c' ... [done]
opening 'r4_chr1_high_mod_tmp06_ad33607634f8' ... [done]
Digests:
sample.idError: segfault from C stack overflow

Do the sampel IDs need to be in a particular format? I created my vcf with plink and used double-id option. IDs are in format: A-[Cohort]-[A#####]. A .gds file is outputed, but I don't know if it's is incorrect due to the segfault.

gds <- seqOpen(r4_chr1_high_mod.gds)
gds
Object of class "SeqVarGDSClass"
File: r4_chr1_high_mod.gds (294.4K)

• [ ] *
  |--+ description [ ] *
  |--+ sample.id { Str8 14306 LZMA_ra(2.94%), 12.6K }
  |--+ variant.id { Int32 3632 LZMA_ra(12.7%), 1.8K }
  |--+ position { Int32 3632 LZMA_ra(62.3%), 8.8K }
  |--+ chromosome { Str8 3632 LZMA_ra(1.62%), 125B }
  |--+ allele { Str8 3632 LZMA_ra(24.4%), 4.0K }
  |--+ genotype [ ] *
  | |--+ data { Bit2 2x14306x3632 LZMA_ra(0.95%), 242.2K }
  | |--+ extra.index { Int32 3x0 LZMA_ra, 18B } *
  | --+ extra { Int16 0 LZMA_ra, 18B }
  |--+ phase [ ]
  | |--+ data { Bit1 14306x3632 LZMA_ra(0.02%), 1.3K }
  | |--+ extra.index { Int32 3x0 LZMA_ra, 18B } *
  | --+ extra { Bit1 0 LZMA_ra, 18B }
  |--+ annotation [ ]
  | |--+ id { Str8 3632 LZMA_ra(28.1%), 16.0K }
  | |--+ qual { Float32 3632 LZMA_ra(0.92%), 141B }
  | |--+ filter { Int32 3632 LZMA_ra(0.92%), 141B }
  | |--+ info [ ]
  | | --+ PR { Bit1 3632 LZMA_ra(18.9%), 93B } *
  | --+ format [ ]
  --+ sample.annotation [ ]

sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS: /cvmfs/priv.accre.vanderbilt.edu/mirror/optimized/sandy_bridge/easybuild/software/MPI/intel/2019.1.144/impi/2018.4.274/R/3.6.0/lib64/R/lib/libR.so
LAPACK: /cvmfs/priv.accre.vanderbilt.edu/mirror/optimized/sandy_bridge/easybuild/software/MPI/intel/2019.1.144/impi/2018.4.274/R/3.6.0/lib64/R/modules/lapack.so

locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] data.table_1.14.8 SeqArray_1.26.2 gdsfmt_1.22.0

loaded via a namespace (and not attached):
[1] zlibbioc_1.32.0 compiler_3.6.0 IRanges_2.20.2
[4] XVector_0.26.0 parallel_3.6.0 GenomicRanges_1.38.0
[7] GenomeInfoDbData_1.2.2 RCurl_1.95-4.12 Biostrings_2.54.0
[10] S4Vectors_0.24.4 BiocGenerics_0.32.0 GenomeInfoDb_1.22.1
[13] bitops_1.0-6 stats4_3.6.0

Thank you,
Alexis

[zhengxwen]

See:
the total number of variants for import: 3,632
This number is too small, parallel=6L does not help at all.
I guess parallel=6L might trigger a bug when merging the data files when the number of variants is too small.

seqVCF2GDS(high_mod_vcf, "r4_chr1_high_mod.gds", parallel=1)

It might solve your problem.

[alexisregelson]

Hello,

I've now tried this with a vcf with a 200k+ varaints. I have successfully converted this vcf to a gds using SNPRelate. However, I am using another software that specifically needs the gds file in SeqArray format, not SNPRelate. But I am still getting the same error: sample.idError: segfault from C stack overflow.

Alexis

[zhengxwen]

Your R version and gdsfmt versions are old.
The recent update was made with a focus on R (>= v4.0).
I suggest using SeqArray GDS format instead of SNPRelate GDS.