Dear all,
I'm trying to obtain annotations for SNPs I analysed in GWASTools. Some time ago it worked fine by using the line
vcfWrite(genoData,vcf.file=vcf.file,id.col="snpID",snp.exclude=snp.excl)
to write the vcf file from the data objects in GWASTools and reading the vcf file with
vcf <- readVcf(vcf.file, genome="hg19")
that is needed to use the VariantAnnotation package.
Now I did the same thing with another set of SNPs and I got the error message:
Fehler: scanVcf: _DNAencode(): invalid DNAString input character: '7' (byte value 55)
path: /home/forschung/data/omni5/GWASTools_corrX/allSNPs.vcf
When I read the "old" vcf file that I created with an older version some months ago, it still works. So I figured that the two vcf files must be different. When I compared them, I found out that the "REF" and "ALT" columns of the old file are empty (containing ".") but there are entries for these columns in the new file.
Could this be the reason for the error? Why is the output different from before? If this is the problem, can I prevent vcfWrite from writing these columns?
Thank you for your help!
> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.2 LTS
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] LC_CTYPE=de_DE.UTF-8 LC_NUMERIC=C LC_TIME=de_DE.UTF-8 LC_COLLATE=de_DE.UTF-8 LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=de_DE.UTF-8
[7] LC_PAPER=de_DE.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] org.Hs.eg.db_3.4.1 TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 GenomicFeatures_1.28.3 AnnotationDbi_1.38.1
[5] BiocInstaller_1.26.0 VariantAnnotation_1.22.1 Rsamtools_1.28.0 Biostrings_2.44.1
[9] XVector_0.16.0 SummarizedExperiment_1.6.3 DelayedArray_0.2.7 matrixStats_0.52.2
[13] GenomicRanges_1.28.3 GenomeInfoDb_1.12.2 IRanges_2.10.2 S4Vectors_0.14.3
[17] GWASTools_1.22.0 Biobase_2.36.2 BiocGenerics_0.22.0
loaded via a namespace (and not attached):
[1] quantsmooth_1.42.0 Rcpp_0.12.11 lattice_0.20-35 zoo_1.8-0 digest_0.6.9 lmtest_0.9-35
[7] plyr_1.8.3 MatrixModels_0.4-1 gdsfmt_1.12.0 RSQLite_1.1-2 ggplot2_2.2.1 zlibbioc_1.22.0
[13] rlang_0.1.1 lazyeval_0.2.0 SparseM_1.77 rpart_4.1-11 Matrix_1.2-10 splines_3.4.0
[19] BiocParallel_1.10.1 GWASExactHW_1.01 RCurl_1.95-4.8 biomaRt_2.32.1 munsell_0.4.2 compiler_3.4.0
[25] rtracklayer_1.36.3 mgcv_1.8-17 nnet_7.3-12 tibble_1.3.3 GenomeInfoDbData_0.99.0 DNAcopy_1.50.1
[31] XML_3.98-1.8 GenomicAlignments_1.12.1 MASS_7.3-47 bitops_1.0-6 grid_3.4.0 nlme_3.1-131
[37] gtable_0.1.2 DBI_0.7 scales_0.4.1 ncdf4_1.16 mice_2.30 sandwich_2.3-4
[43] tools_3.4.0 BSgenome_1.44.0 survival_2.41-3 colorspace_1.2-4 memoise_0.2.1 logistf_1.22
[49] quantreg_5.33
You could look at the REF column for non-DNA characters. As an aside, it's really unfortunate that two Bioconductor tools are communicating through files.
You make a good point about the file communication, but unfortunately VariantAnnotation only takes VCF files as input, so it's up to every other package author to write a coercion method to a VCF object. Based on my experience with doing this for the SeqArray package, that's a substantial amount of work, since the VCF objects are quite complex. Realistically it's not going to happen for GWASTools anytime soon.
A method in VariantAnnotation for creating a VCF object from simpler inputs (integer/character vectors for the fixed columns, a matrix for the genotypes) would make it easier for other packages to communicate with VariantAnnotation.
If it's easier to write (or at least maintain) your own VCF writer than to construct a VCF object, then indeed the constructor must be very difficult to use. I will have a go at a coercion from GWASTools objects to VCF with an eye towards the difficulties.