VariantAnnotation
VariantAnnotation copied to clipboard
Bioconductor
`readVcf` doesn't handle quoted `=` in `hapmap_exome_chr22.vcf.gz` correctly
The header line for GATKCommandLine causes problems:
##GATKCommandLine=<ID=ApplyRecalibration,Version=3.3-0-g37228af,Date="Thu May 07 02:55:12 EDT 2015",Epoch=1430981712637,CommandLineOptions="analysis_type=ApplyRecalibration input_file=[] showFullBamList=false read_buffer_size=null phone_home
=NO_ET gatk_key=/isilon/sequencing/CIDRSeqSuiteSoftware/gatk/GATK_2/lee.watkins_jhmi.edu.key tag=NA read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/isilo
n/sequencing/GATK_resource_bundle/1.5/b37/human_g1k_v37_decoy.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=
1000 baq=OFF baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quanti
ze_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=nu
ll disable_auto_index_creation_and_locking_when_reading_rods=true no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false nu
m_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=
false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 logging_level=INFO log_to_file=null help=false version=false input=[(RodBinding name=input source=/isilon/sequencing/Seq_Proj//Amos_LungCa_SeqWholeExome
Plus_040413_1/MULTI_SAMPLE/Amos_LungCa_SeqWholeExomePlus_040413_1.GATK-3.3-0.2050samples.raw.HC.INDEL.vcf)] recal_file=(RodBinding name=recal_file source=/isilon/sequencing/Seq_Proj//Amos_LungCa_SeqWholeExomePlus_040413_1/MULTI_SAMPLE/Amos_L
ungCa_SeqWholeExomePlus_040413_1.GATK-3.3-0.2050samples.HC.INDEL.recal) tranches_file=/isilon/sequencing/Seq_Proj/Amos_LungCa_SeqWholeExomePlus_040413_1/MULTI_SAMPLE/Amos_LungCa_SeqWholeExomePlus_040413_1.GATK-3.3-0.2050samples.HC.INDEL.tran
ches out=org.broadinstitute.gatk.engine.io.stubs.VariantContextWriterStub ts_filter_level=99.0 lodCutoff=null ignore_filter=null ignore_all_filters=false excludeFiltered=false mode=INDEL filter_reads_with_N_cigar=false filter_mismatching_bas
e_and_quals=false filter_bases_not_stored=false">
Calling readVcf gives me:
fl <- system.file("extdata", "hapmap_exome_chr22.vcf.gz", package="VariantAnnotation")
library(VariantAnnotation)
first <- readVcf(fl)
meta(header(first))$GATKCommandLine$CommandLineOptions
## [1] "\"analysis_type"
Oops, look like it just cut off after the = inside the quote. This causes no end of problems during a write round-trip:
writeVcf(first, out)
roundtrip <- readVcf(out, row.names=FALSE)
## [W::bcf_hdr_parse_line] Incomplete header line, trying to proceed anyway:
## [##GATKCommandLine=<ID=ApplyRecalibration,Version=3.3-0-g37228af,Date="Thu May 07 02:55:12 EDT 2015",Epoch=1430981712637,CommandLineOptions="analysis_type>
Session information
R Under development (unstable) (2023-11-29 r85646)
Platform: aarch64-apple-darwin22.5.0
Running under: macOS Ventura 13.6.1
Matrix products: default
BLAS: /Users/luna/Software/R/trunk/lib/libRblas.dylib
LAPACK: /Users/luna/Software/R/trunk/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] VariantAnnotation_1.49.2 Rsamtools_2.19.2
[3] Biostrings_2.71.1 XVector_0.43.0
[5] SummarizedExperiment_1.33.1 Biobase_2.63.0
[7] GenomicRanges_1.55.1 GenomeInfoDb_1.39.5
[9] IRanges_2.37.0 S4Vectors_0.41.3
[11] MatrixGenerics_1.15.0 matrixStats_1.2.0
[13] BiocGenerics_0.49.1
loaded via a namespace (and not attached):
[1] KEGGREST_1.43.0 rjson_0.2.21 lattice_0.22-5
[4] vctrs_0.6.5 tools_4.4.0 bitops_1.0-7
[7] generics_0.1.3 curl_5.2.0 parallel_4.4.0
[10] tibble_3.2.1 fansi_1.0.6 AnnotationDbi_1.65.2
[13] RSQLite_2.3.4 blob_1.2.4 pkgconfig_2.0.3
[16] Matrix_1.6-4 BSgenome_1.71.1 dbplyr_2.4.0
[19] lifecycle_1.0.4 GenomeInfoDbData_1.2.11 compiler_4.4.0
[22] stringr_1.5.1 progress_1.2.3 codetools_0.2-19
[25] yaml_2.3.8 RCurl_1.98-1.13 pillar_1.9.0
[28] crayon_1.5.2 BiocParallel_1.37.0 DelayedArray_0.29.0
[31] cachem_1.0.8 abind_1.4-5 tidyselect_1.2.0
[34] digest_0.6.33 stringi_1.8.3 restfulr_0.0.15
[37] dplyr_1.1.4 biomaRt_2.59.0 fastmap_1.1.1
[40] grid_4.4.0 cli_3.6.2 SparseArray_1.3.2
[43] magrittr_2.0.3 S4Arrays_1.3.1 GenomicFeatures_1.55.1
[46] XML_3.99-0.16 utf8_1.2.4 rappdirs_0.3.3
[49] filelock_1.0.3 prettyunits_1.2.0 bit64_4.0.5
[52] httr_1.4.7 bit_4.0.5 png_0.1-8
[55] hms_1.1.3 memoise_2.0.1 BiocIO_1.13.0
[58] BiocFileCache_2.11.1 rtracklayer_1.63.0 rlang_1.1.2
[61] glue_1.6.2 DBI_1.2.0 xml2_1.3.6
[64] R6_2.5.1 GenomicAlignments_1.39.0 zlibbioc_1.49.0
Not sure if this helps. If we use TabixFile(fl) for readVcf, the round trip will work. Is the embedded "=" possible in a valid VCF?
The problem seems to be in Rsamtools:::.bcfHeaderAsSimpleList. There is a regular expression check that excludes a header line containing a period: header line beginning GATKCommandLine.SelectVariants is excluded, IMHO erroneously. Pressing on...
Somewhere in here, in the !idx lines, there is mishandling of embedded '='
rex <- "^([[:alnum:]]+)=<(.*)>"
lines <- grep(rex, h, value = TRUE)
if (length(lines) > 0) {
if (any(duplicated(lines))) {
warning("duplicate INFO and FORMAT header lines will be ignored")
lines <- lines[duplicated(lines) == FALSE]
}
}
tags <- sub(rex, "\\1", lines)
keyval0 <- sub(rex, "\\2", lines)
keyval1 <- rep(NA_character_, length(keyval0))
keyval <- list()
idx <- tags %in% c("INFO", "FORMAT", "FILTER", "ALT")
keyval1[idx] <- strsplit(keyval0[idx], ",(?=(ID|Number|Type)=[[:alnum:]]*)|,(?=Description=\".*?\")",
perl = TRUE)
keyval[idx] <- lapply(which(idx), function(i, keyval1) strsplit(keyval1[[i]],
"(?<=[ID|Number|Type|Description])=", perl = TRUE), keyval1)
keyval1[!idx] <- strsplit(keyval0[!idx], ",(?=[[:alnum:]]+=)",
perl = TRUE)
keyval[!idx] <- lapply(which(!idx), function(i, keyval1) {
strsplit(keyval1[[i]], "(?<=[[:alnum:]])=", perl = TRUE)
}, keyval1)
The issue comes down to parsing this record (or any like it) into a nice DataFrame ... the code didn't anticipate a record of this type and it clearly fails to process it correctly, truncating the info after analysis_type. I do not have time right now to fix it.
##GATKCommandLine.SelectVariants=<ID=SelectVariants,Version=3.4-46-gbc02625,Date="Wed Sep 23 12:04:45 GMT-08:00 2015",Epoch=1443038685033,CommandLineOptions="analysis_type=SelectVariants input_file=[] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/projects/cidr/Amos/amos_cidr/Analysis_Pipeline_Files/human_g1k_v37_decoy.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=1000 baq=OFF baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=false no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 logging_level=INFO log_to_file=null help=false version=false variant=(RodBinding name=variant source=/projects/cidr/Amos/amos_cidr/MultiSampleVCF/withGenotypeRefinement_Redone/Amos_LungCa_SeqWholeExomePlus_040413_1.GATK-3.3-0.2047samples.BEDsuperset.VQSR.1KG.ExAC3.REFINED.06192015.vcf.gz) discordance=(RodBinding name= source=UNBOUND) concordance=(RodBinding name= source=UNBOUND) out=/projects/geneva/gcc-fs2/stephanie/Code/Bioconductor/amos_example_data/hapmap_exome_chr22.vcf.gz sample_name=[] sample_expressions=null sample_file=[amos_hapmaps.txt] exclude_sample_name=[] exclude_sample_file=[] exclude_sample_expressions=[] selectexpressions=[] invertselect=false excludeNonVariants=true excludeFiltered=false preserveAlleles=true removeUnusedAlternates=false restrictAllelesTo=ALL keepOriginalAC=false keepOriginalDP=false mendelianViolation=false invertMendelianViolation=false mendelianViolationQualThreshold=0.0 select_random_fraction=0.0 remove_fraction_genotypes=0.0 selectTypeToInclude=[] selectTypeToExclude=[] keepIDs=amos_chr22_subset.txt excludeIDs=null fullyDecode=false justRead=false maxIndelSize=2147483647 minIndelSize=0 maxFilteredGenotypes=2147483647 minFilteredGenotypes=0 maxFractionFilteredGenotypes=1.0 minFractionFilteredGenotypes=0.0 setFilteredGtToNocall=false ALLOW_NONOVERLAPPING_COMMAND_LINE_SAMPLES=false filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">
I consulted the standard, section 1.4, where it says
The values of optional fields must be written as quoted strings, even for numeric values.
I wrote three test cases that seem to be consistent with the standard
h = c(
'key1=<ID=foo1,clo="a=1,i=[],s=false",h=baz>',
'key2=<ID=foo2,clo="a=1 i=[] s=false",h=baz>',
'key3=<ID=foo3, date="Thu May 07", clo="a=1,i=[],s=false",h=baz>'
)
the tag / key-value pairs are extracted like in the code @vjcitn notes
rex <- "^([[:alnum:]]+)=<(.*)>"
tags <- sub(rex, "\\1", h)
keyval0 <- sub(rex, "\\2", h)
Within keyval0, it seems we want to find a 'tag' [[:alnum:]]+, an equal sign =, and EITHER an alphanumeric value [[:alnum:]]+ OR a quoted (with double quotes?) string "[^"]+". So extract tags and values
rex0 <- '[[:alnum:]]+=([[:alnum:]]+|"[^"]+")'
m <- gregexpr('[[:alnum:]]+=([[:alnum:]]+|"[^"]+")', keyval0)
keyval1 <- regmatches(keyval0, m)
This gives us the correct (?) key/value pairs for each test case
> keyval1
[[1]]
[1] "ID=foo1" "clo=\"a=1,i=[],s=false\""
[3] "h=baz"
[[2]]
[1] "ID=foo2" "clo=\"a=1 i=[] s=false\""
[3] "h=baz"
[[3]]
[1] "ID=foo3" "date=\"Thu May 07\""
[3] "clo=\"a=1,i=[],s=false\"" "h=baz"
and the individual tags / values can be extracted with relatively simple regular expressions
rex1 <- '^([[:alnum:]]+)=(.*)'
## nested tags
lapply(keyval1, sub, pattern = rex1, replacement = "\\1")
## nested keys
lapply(keyval1, sub, pattern = rex1, replacement = "\\2")
The example @LTLA provides seems to violate the standard, since Version=3.3-0-g37228af is an optional string, and it is not quoted. I added a test case.
h = c(
'key1=<ID=foo1,clo="a=1,i=[],s=false",h=baz>',
'key2=<ID=foo2,clo="a=1 i=[] s=false",h=baz>',
'key3=<ID=foo3, date="Thu May 07", clo="a=1,i=[],s=false",h=baz>',
'key4=<ID=foo3, version=3.14-15, date="Thu May 07", clo="a=1,i=[],s=false",h=baz>'
)
We could relax the pattern for finding nested values, e.g., by matching EITHER any string not containing a comma [^,]+ OR a quoted string "[^"]+"
rex0 <- '[[:alnum:]]+=([^,]+|"[^"]+")'
Thanks for looking into this. If you have time to introduce this enhancement, that would be great. Otherwise, eventually I will study and attempt to do it.
Another simple fix would be to just update hapmap_exome_chr22.vcf.gz so that it doesn't have this non-compliant header. I only stumbled across this because I was rummaging around VariantAnnotation's VCFs to find some test files.
IMO if the file isn't standards-compliant, then perhaps we shouldn't work too hard to accommodate it.
I am going to take the step of removing offending records from hapmap_exome_chr22.vcf.gz now.