writeVcf save invalid vcf with metadata out of order
1
0
Entering edit mode
@alessandropastore-10879
Last seen 3.9 years ago
If I read a valid vcf with the last 3.8 Bioconductor version 1.28.1

param_T <- ScanVcfParam(sample = samples_names[1],  geno = c("AD", "MQ", "DP", "AF"), info = c("DP","AF", "CALLERS"))
  Tumor <- readVcf(file = vcf_file[i], genome = "GRCH37", param = param_T)
  

 

head of  vcf_file[i
##fileformat=VCFv4.2
##contig=<ID=1,length=249250621>
##contig=<ID=2,length=243199373>
##contig=<ID=3,length=198022430>
##contig=<ID=4,length=191154276>
##contig=<ID=5,length=180915260>
##contig=<ID=6,length=171115067>

 

and then I save it with 

writeVcf(Tumor, paste0("./vcf/",samples_names[1], ".vcf"), index = F)

the header is out of order ##fileformat=VCFv4.2 is not on the first line and there is a Strange SAMPLE field 

Thanks for you support!

##SAMPLE=Tumor
##commandline="/home/pastore/data/bcbio/anaconda/bin/freebayes -f /lila/data/abdelwao-lab/pastorea/bcbio/genomes/Hsapiens/GRCh37/seq/GRCh37.fa --genotype-qualities --strict-vcf --ploidy 2 --targets /lila/data/abdelwao-lab/pastorea/Project/Ben/histiocytosis/hempact/work/freebayes/1/B02_816_T1-1_0_249250621-regions-nolcr.bed --min-repeat-entropy 1 --no-partial-observations --min-alternate-fraction 0.05 --pooled-discrete --pooled-continuous --report-genotype-likelihood-max --allele-balance-priors-off /lila/data/abdelwao-lab/pastorea/Project/Ben/histiocytosis/hempact/work/align/C_K0Y3K8_T001_d/C_K0Y3K8_T001_d-sort.bam /lila/data/abdelwao-lab/pastorea/Project/Ben/histiocytosis/hempact/work/align/PON/PON-sort.bam"
##fileDate=20181115
##fileformat=VCFv4.2
##phasing=none
variantannotation • 463 views
ADD COMMENT
0
Entering edit mode
@valerie-obenchain-4275
Last seen 8 months ago
United States

Thanks for reporting the inconsistent ordering of the header lines in writeVcf(). I've checked in a fix to VariantAnnotation 1.28.2 and 1.29.1 to enforce 'fileformat' as the first line.

As for the SAMPLE line in the header, this must be a carry over from your original file, e.g.,

header(scanVcfHeader(file))

If this is not the case, please send me the vcf file (valerie.obenchain@roswellpark.org) so I can try to reproduce this.

Valerie

ADD COMMENT

Login before adding your answer.

Traffic: 451 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6