gvcf fill ref allele from fasta

doc/bcftools.1

@@ -2,12 +2,12 @@
 .\"     Title: bcftools
 .\"    Author: [see the "AUTHORS" section]
 .\" Generator: DocBook XSL Stylesheets v1.76.1 <http://docbook.sf.net/>
-.\"      Date: 2015-03-20 11:57 GMT
+.\"      Date: 04/22/2015
 .\"    Manual: \ \&
 .\"    Source: \ \&
 .\"  Language: English
 .\"
-.TH "BCFTOOLS" "1" "2015\-03\-20 11:57 GMT" "\ \&" "\ \&"
+.TH "BCFTOOLS" "1" "04/22/2015" "\ \&" "\ \&"
 .\" -----------------------------------------------------------------
 .\" * Define some portability stuff
 .\" -----------------------------------------------------------------
@@ -41,7 +41,6 @@ Most commands accept VCF, bgzipped VCF and BCF with filetype detected automatica
 BCFtools is designed to work on a stream\&. It regards an input file "\-" as the standard input (stdin) and outputs to the standard output (stdout)\&. Several commands can thus be combined with Unix pipes\&.
 .SS "VERSION"
 .sp
-This manual page was last updated \fB2015\-03\-20 11:57 GMT\fR and refers to bcftools git version \fB1\&.2\-6\-ga8d7fe9+\fR\&.
 .SS "BCF1"
 .sp
 The BCF1 format output by versions of samtools <= 0\&.1\&.19 is \fBnot\fR compatible with this version of bcftools\&. To read BCF1 files one can use the view command from old versions of bcftools packaged with samtools versions <= 0\&.1\&.19 to convert to VCF, which can then be read by this version of bcftools\&.
@@ -1065,6 +1064,11 @@ tag to take values for \&.gen file: GT,PL,GL,GP
 .RS 4
 convert gVCF to VCF, expanding REF blocks into sites\&. Only sites with FILTER set to "PASS" or "\&." will be expanded\&.
 .RE
+.PP
+\fB\-f, \-\-fasta\-ref\fR \fIfile\fR
+.RS 4
+reference sequence in fasta format\&. Must be indexed with samtools faidx
+.RE
 .RE
 .sp
 .it 1 an-trap
@@ -1175,7 +1179,7 @@ comma\-separated list of fields in the input file\&. In the current version, the
 .PP
 \fB\-f, \-\-fasta\-ref\fR \fIfile\fR
 .RS 4
-reference sequence in fasta format
+reference sequence in fasta format\&. Must be indexed with samtools faidx
 .RE
 .PP
 \fB\-s, \-\-samples\fR \fILIST\fR

doc/bcftools.txt

@@ -606,6 +606,9 @@ Create consensus sequence by applying VCF variants to a reference fasta file.
     convert gVCF to VCF, expanding REF blocks into sites. Only sites
     with FILTER set to "PASS" or "." will be expanded.
 
+*-f, --fasta-ref* 'file'::
+    reference sequence in fasta format. Must be indexed with samtools faidx
+
 ==== HAPS/SAMPLE conversion:
 *--hapsample2vcf* 'prefix' or 'haps-file','sample-file'::
     convert from haps/sample format to VCF. The columns of .haps file are
@@ -674,7 +677,7 @@ Create consensus sequence by applying VCF variants to a reference fasta file.
     are not supported yet, missing data can be indicated with "--".
 
 *-f, --fasta-ref* 'file'::
-    reference sequence in fasta format
+    reference sequence in fasta format. Must be indexed with samtools faidx
 
 *-s, --samples* 'LIST'::
     list of sample names. See *<<common_options,Common Options>>*

test/convert.gvcf.vcf

@@ -22,38 +22,36 @@
 ##FILTER=<ID=HighSNVSB,Description="SNV strand bias value (SNVSB) exceeds 10">
 ##FILTER=<ID=HighREFREP,Description="Locus contains an indel allele occurring in a homopolymer or dinucleotide track with a reference repeat greater than 8">
 ##FILTER=<ID=HighDepth,Description="Locus depth is greater than 3x the mean chromosome depth">
-##reference=file:///illumina/scripts/clia/Genomes/Homo_sapiens/UCSC/hg19_rCRS/Sequence/WholeGenomeFasta/genome.fa
-##contig=<ID=chr22,length=51304566>
+##contig=<ID=22,length=450>
 ##SnvTheta=0.001
 ##IndelTheta=0.0001
-##MaxDepth_chr1=114.18
-##MaxDepth_chr10=131.73
-##MaxDepth_chr11=117.27
-##MaxDepth_chr12=116.97
-##MaxDepth_chr13=102.24
-##MaxDepth_chr14=101.55
-##MaxDepth_chr15=95.22
-##MaxDepth_chr16=111.33
-##MaxDepth_chr17=112.59
-##MaxDepth_chr18=121.86
-##MaxDepth_chr19=111.12
-##MaxDepth_chr2=121.83
-##MaxDepth_chr20=111.24
-##MaxDepth_chr21=98.43
-##MaxDepth_chr22=76.23
-##MaxDepth_chr3=120.09
-##MaxDepth_chr4=124.50
-##MaxDepth_chr5=119.82
-##MaxDepth_chr6=122.22
-##MaxDepth_chr7=120.27
-##MaxDepth_chr8=120.45
-##MaxDepth_chr9=102.48
-##MaxDepth_chrM=7005.66
-##MaxDepth_chrX=61.05
-##MaxDepth_chrY=37.17
+##MaxDepth_1=114.18
+##MaxDepth_10=131.73
+##MaxDepth_11=117.27
+##MaxDepth_12=116.97
+##MaxDepth_13=102.24
+##MaxDepth_14=101.55
+##MaxDepth_15=95.22
+##MaxDepth_16=111.33
+##MaxDepth_17=112.59
+##MaxDepth_18=121.86
+##MaxDepth_19=111.12
+##MaxDepth_2=121.83
+##MaxDepth_20=111.24
+##MaxDepth_21=98.43
+##MaxDepth_22=76.23
+##MaxDepth_3=120.09
+##MaxDepth_4=124.50
+##MaxDepth_5=119.82
+##MaxDepth_6=122.22
+##MaxDepth_7=120.27
+##MaxDepth_8=120.45
+##MaxDepth_9=102.48
+##MaxDepth_M=7005.66
+##MaxDepth_X=61.05
+##MaxDepth_Y=37.17
 ##FILTER=<ID=IndelSizeFilter,Description="Indel is outside reportable size range. Insertion range: [1,3], Deletion range: [1,11]">
 ##gvcftools_version="0.16"
-##gvcftools_cmdline="/illumina/scripts/clia/workflows/IsisWorkflow/IsisWorkflow_v2.0.13/bin/set_haploid_region --ref /illumina/scripts/clia/Genomes/Homo_sapiens/UCSC/hg19_rCRS/Sequence/WholeGenomeFasta/genome.fa --region-file /illumina/scripts/clia/workflows/IsisWorkflow/IsisWorkflow_v2.0.13/bin/../data/het_mask/ncbi37/male.bed"
 ##FILTER=<ID=HAPLOID_CONFLICT,Description="Locus has heterozygous genotype in a haploid region.">
 ##FORMAT=<ID=OPL,Number=.,Type=Integer,Description="Original PL value before ploidy correction">
 ##INFO=<ID=phastCons,Number=0,Type=Flag,Description="overlaps a phastCons element">
@@ -74,34 +72,34 @@
 ##INFO=<ID=CLNDSDBID,Number=.,Type=String,Description="Colon-delimited list of variant disease database identifier(s). Multiple values from a single database are pipe-delimited">
 ##INFO=<ID=CSQ,Number=A,Type=String,Description="Consequence type as predicted by VEP. Format: Allele|Gene|Feature|Feature_type|Consequence|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|EXON|INTRON|HGNC|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE|DISTANCE|CANONICAL|SIFT|PolyPhen|GMAF|ENSP|DOMAINS|CCDS|HGVSc|HGVSp|CELL_TYPE">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SAMPLE99
-chr22	1	.	N	.	0	LowGQX	END=16050039;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	.:.:0:0
-chr22	16050040	.	C	.	0	LowGQX	END=16050050;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:5:2:0
-chr22	16050051	.	C	.	0	LowGQX	END=16050056;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:10:4:0
-chr22	16050057	.	C	.	0	LowGQX	END=16050072;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:15:6:0
-chr22	16050073	.	G	.	0	LowGQX	END=16050080;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:21:8:0
-chr22	16050081	.	C	.	0	PASS	END=16050085;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:30:11:0
-chr22	16050086	.	G	.	0	PASS	END=16050102;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:42:15:0
-chr22	16050103	.	T	.	0	PASS	END=16050115;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:54:19:0
-chr22	16050116	.	G	C	23	LowGQX	SNVSB=0;SNVHPOL=2	GT:GQ:GQX:DP:DPF:AD	0/1:56:23:22:0:16,6
-chr22	16050117	.	T	.	0	PASS	END=16050131;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:69:24:0
-chr22	16050132	.	A	.	0	PASS	END=16050170;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:84:29:0
-chr22	16050171	.	G	.	0	PASS	END=16050215;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:90:31:0
-chr22	16050216	.	T	.	0	PASS	END=16050217;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:93:32:2
-chr22	16050218	.	T	.	0	PASS	END=16050234;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:102:35:0
-chr22	16050235	.	T	.	0	PASS	.	GT:GQX:DP:DPF	0/0:72:36:0
-chr22	16050236	.	T	.	0	PASS	END=16050261;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:87:30:0
-chr22	16050262	.	A	.	0	PASS	END=16050279;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:66:23:0
-chr22	16050280	.	G	.	0	PASS	END=16050287;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:87:30:0
-chr22	16050288	.	A	.	0	PASS	END=16050299;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:114:39:0
-chr22	16050300	.	G	.	0	PASS	END=16050309;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:150:51:0
-chr22	16050310	.	A	.	0	PASS	END=16050344;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:166:63:0
-chr22	16050345	.	C	.	0	PASS	END=16050346;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:220:74:0
-chr22	16050347	.	T	.	0	HighDepth	END=16050372;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:229:77:0
-chr22	16050373	.	T	.	0	HighDepth	END=16050407;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:241:81:0
-chr22	16050408	.	T	.	0	HighDepth	END=16050414;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:200:75:0
-chr22	16050415	.	G	.	0	PASS	END=16050418;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:217:73:0
-chr22	16050419	.	C	.	0	HighDepth	END=16050420;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:223:75:0
-chr22	16050421	.	T	.	0	PASS	.	GT:GQX:DP:DPF	0/0:223:75:1
-chr22	16050422	.	T	.	0	HighDepth	END=16050427;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:223:75:0
-chr22	16050428	.	T	.	0	PASS	END=16050450;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:178:60:0
-chr22	16050451	.	C	.	0	PASS	END=16050466;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:142:54:0
+22	1	.	N	.	0	LowGQX	END=9;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	.:.:0:0
+22	10	.	C	.	0	LowGQX	END=20;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:5:2:0
+22	21	.	C	.	0	LowGQX	END=26;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:10:4:0
+22	27	.	C	.	0	LowGQX	END=42;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:15:6:0
+22	43	.	G	.	0	LowGQX	END=50;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:21:8:0
+22	51	.	C	.	0	PASS	END=55;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:30:11:0
+22	56	.	G	.	0	PASS	END=72;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:42:15:0
+22	73	.	T	.	0	PASS	END=85;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:54:19:0
+22	86	.	G	C	23	LowGQX	SNVSB=0;SNVHPOL=2	GT:GQ:GQX:DP:DPF:AD	0/1:56:23:22:0:16,6
+22	87	.	T	.	0	PASS	END=101;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:69:24:0
+22	102	.	A	.	0	PASS	END=140;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:84:29:0
+22	141	.	G	.	0	PASS	END=185;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:90:31:0
+22	186	.	T	.	0	PASS	END=187;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:93:32:2
+22	188	.	T	.	0	PASS	END=204;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:102:35:0
+22	205	.	T	.	0	PASS	.	GT:GQX:DP:DPF	0/0:72:36:0
+22	206	.	T	.	0	PASS	END=231;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:87:30:0
+22	232	.	A	.	0	PASS	END=249;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:66:23:0
+22	250	.	G	.	0	PASS	END=257;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:87:30:0
+22	258	.	A	.	0	PASS	END=269;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:114:39:0
+22	270	.	G	.	0	PASS	END=279;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:150:51:0
+22	280	.	A	.	0	PASS	END=314;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:166:63:0
+22	315	.	C	.	0	PASS	END=316;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:220:74:0
+22	317	.	T	.	0	HighDepth	END=342;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:229:77:0
+22	343	.	T	.	0	HighDepth	END=377;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:241:81:0
+22	378	.	T	.	0	HighDepth	END=384;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:200:75:0
+22	385	.	G	.	0	PASS	END=388;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:217:73:0
+22	389	.	C	.	0	HighDepth	END=390;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:223:75:0
+22	391	.	T	.	0	PASS	.	GT:GQX:DP:DPF	0/0:223:75:1
+22	392	.	T	.	0	HighDepth	END=397;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:223:75:0
+22	398	.	T	.	0	PASS	END=420;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:178:60:0
+22	421	.	C	.	0	PASS	END=450;BLOCKAVG_min30p3a	GT:GQX:DP:DPF	0/0:142:54:0

test/gvcf.fa

@@ -0,0 +1,9 @@
+>22
+NNNNNNNNCCTTGGCCAAGTCACTTCCTCCTTCAGGAACATTGCAGTGGGCCTAAGTGCC
+TCCTCTCGGGACTGGTATGGGGACGGTCATGCAATCTGGACAACATTCACCTTTAAAAGT
+TTATTGATCTTTTGTGACATGCACGTGGGTTCCCAGTAGCAAGAAACTAAAGGGTCGCAG
+GCCGGTTTCTGCTAATTTCTTTAATTCCAAGACAGTCTCAAATATTTTCTTATTAACTTC
+CTGGAGGGAGGCTTATCATTCTCTCTTTTGGATGATTCTAAGTACCAGCTAAAATACAGC
+TATCATTCATTTTCCTTGATTTGGGAGCCTAATTTCTTTAATTTAGTATGCAAGAAAACC
+AATTTGGAAATATCAACTGTTTTGGAAACCTTAGACCTAGGTCATCCTTAGTAAGATCTT
+CCCATTTATATAAATACTTGCAAGTAGTAGTGCCATAATT

test/gvcf.fa.fai

@@ -0,0 +1 @@
+22	460	4	60	61

test/test.pl

@@ -174,7 +174,7 @@
 test_vcf_convert($opts,in=>'convert',out=>'convert.hls.samples',args=>'-h .,.,-');
 test_vcf_convert($opts,in=>'convert',out=>'convert.hs.hap',args=>'--hapsample -,.');
 test_vcf_convert($opts,in=>'convert',out=>'convert.hs.sample',args=>'--hapsample .,-');
-test_vcf_convert_gvcf($opts,in=>'convert.gvcf',out=>'convert.gvcf.out',args=>'--gvcf2vcf');
+test_vcf_convert_gvcf($opts,in=>'convert.gvcf',out=>'convert.gvcf.out',fa=>'gvcf.fa',args=>'--gvcf2vcf');
 test_vcf_convert_tsv2vcf($opts,in=>'convert.23andme',out=>'convert.23andme.vcf',args=>'-c ID,CHROM,POS,AA -s SAMPLE1',fai=>'23andme');
 test_vcf_consensus($opts,in=>'consensus',out=>'consensus.1.out',fa=>'consensus.fa',mask=>'consensus.tab',args=>'');
 test_vcf_consensus_chain($opts,in=>'consensus',out=>'consensus.1.chain',chain=>'consensus.1.chain',fa=>'consensus.fa',mask=>'consensus.tab',args=>'');
@@ -472,8 +472,8 @@ sub test_vcf_convert_gvcf
 {
     my ($opts,%args) = @_;
     bgzip_tabix_vcf($opts,$args{in});
-    test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools convert $args{args} $$opts{tmp}/$args{in}.vcf.gz | grep -v ^##bcftools");
-    test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools view -Ob $$opts{tmp}/$args{in}.vcf.gz | $$opts{bin}/bcftools convert $args{args} | grep -v ^##bcftools");
+    test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools convert $args{args} -f $$opts{path}/$args{fa} $$opts{tmp}/$args{in}.vcf.gz | grep -v ^##bcftools");
+    test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools view -Ob $$opts{tmp}/$args{in}.vcf.gz | $$opts{bin}/bcftools convert $args{args} -f $$opts{path}/$args{fa} | grep -v ^##bcftools");
 }
 sub test_vcf_convert_tsv2vcf
 {

vcfconvert.c

@@ -71,6 +71,7 @@ struct _args_t
 
 static void destroy_data(args_t *args)
 {
+    if ( args->ref ) fai_destroy(args->ref);
     if ( args->convert) convert_destroy(args->convert);
     if ( args->filter ) filter_destroy(args->filter);
     free(args->samples);
@@ -1073,8 +1074,8 @@ static int tsv_setter_aa(tsv_t *tsv, bcf1_t *rec, void *usr)
 
 static void tsv_to_vcf(args_t *args)
 {
-    if ( !args->ref_fname ) error("Missing the --ref option\n");
-    if ( !args->sample_list ) error("Missing the --samples option\n");
+    if ( !args->ref_fname ) error("--tsv2vcf requires the --fasta-ref option\n");
+    if ( !args->sample_list ) error("--tsv2vcf requires the --samples option\n");
 
     args->ref = fai_load(args->ref_fname);
     if ( !args->ref ) error("Could not load the reference %s\n", args->ref_fname);
@@ -1125,7 +1126,6 @@ static void tsv_to_vcf(args_t *args)
     if ( hts_close(in_fh) ) error("Close failed: %s\n", args->infname);
     free(line.s);
 
-    fai_destroy(args->ref);
     bcf_hdr_destroy(args->header);
     hts_close(out_fh);
     tsv_destroy(tsv);
@@ -1167,6 +1167,11 @@ static void vcf_to_vcf(args_t *args)
 
 static void gvcf_to_vcf(args_t *args)
 {
+    if ( !args->ref_fname ) error("--gvcf2vcf requires the --fasta-ref option\n");
+
+    args->ref = fai_load(args->ref_fname);
+    if ( !args->ref ) error("Could not load the fai index for reference %s\n", args->ref_fname);
+
     open_vcf(args,NULL);
     htsFile *out_fh = hts_open(args->outfname,hts_bcf_wmode(args->output_type));
     if ( !out_fh ) error("Failed to open: %s\n", args->outfname);
@@ -1202,10 +1207,15 @@ static void gvcf_to_vcf(args_t *args)
             continue;
         }
         bcf_update_info_int32(hdr,line,"END",NULL,0);
-        int pos;
+        int pos, len;
         for (pos=line->pos; pos<itmp[0]; pos++)
         {
             line->pos = pos;
+            char *ref = faidx_fetch_seq(args->ref, (char*)bcf_hdr_id2name(hdr,line->rid), line->pos, line->pos, &len);
+            if ( !ref ) error("faidx_fetch_seq failed at %s:%d\n", bcf_hdr_id2name(hdr,line->rid), line->pos+1);
+            // we have already checked above that there is only one allele,
+            // so fine to just update alleles with the ref allele from the fasta
+            bcf_update_alleles_str(hdr, line, &ref[0]);
             bcf_write(out_fh,hdr,line);
         }
     }
@@ -1243,7 +1253,8 @@ static void usage(void)
     fprintf(stderr, "       --vcf-ids               output VCF IDs in second column instead of CHROM:POS_REF_ALT\n");
     fprintf(stderr, "\n");
     fprintf(stderr, "gVCF conversion:\n");
-    fprintf(stderr, "       --gvcf2vcf              \n");
+    fprintf(stderr, "       --gvcf2vcf              expand gVCF reference blocks\n");
+    fprintf(stderr, "   -f, --fasta-ref <file>      reference sequence in fasta format\n");
     fprintf(stderr, "\n");
     fprintf(stderr, "HAP/SAMPLE conversion (output from SHAPEIT):\n");
     fprintf(stderr, "       --hapsample2vcf <...>   <prefix>|<haps-file>,<sample-file>\n");