Break_blocks conversion of gvcf -> vcf
1
0
Entering edit mode
21 months ago
dec986 ▴ 300

I have a list of VCFs that I'm converting from gVCF to VCF thus:

zcat HG00096.final_combined.g.vcf.gz | gvcftools-0.17.0/bin/break_blocks --region-file gsa_positions.chromEnd_plus1.bed --ref g1k_v37/human_g1k_v37.fasta --exclude-off-target > vcf/HG00096.bb.vcf 2> break_blocks_err/HG00096.bb.err

the bed file looks like

1   58813   58814   GSA-rs114420996 .   G   A   PASS    .   GT:GQ   ./.:0.0
1   565507  565508  GSA-rs9283150   .   G   A   PASS    .   GT:GQ   ./.:0.0
1   567091  567092  GSA-rs9326622   .   T   C   PASS    .   GT:GQ   ./.:0.0
1   726911  726912  GSA-1:726912    .   A   G   PASS    .   GT:GQ   0/0:0.27129138
1   727840  727841  GSA-rs116587930 .   G   A   PASS    .   GT:GQ   0/0:0.6141555

but the problem is that the output files are missing alternate alleles, e.g.

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HG00126
1   58814   .   G   <NON_REF>   .   .   .   GT:DP:GQ:MIN_DP:PL  0/0:31:84
:29:0,84,1192
1   565508  .   G   <NON_REF>   .   .   .   GT:DP:GQ:MIN_DP:PL  0/0:27:50
:27:0,50,857
1   567092  .   T   <NON_REF>   .   .   .   GT:DP:GQ:MIN_DP:PL  0/0:16:45
:16:0,45,675

I've also tried piping to extract_variants from Converting Gvcf Files Into Vcf but this isn't working either.

how can I get the output to not have the NON_REF so that it looks like a correct VCF?

VCF • 407 views
ADD COMMENT
0
Entering edit mode
21 months ago
ATpoint 55k

If the bed and the vcf have identical sort order simply take the ALT column from the BED to replace that column:

paste <(cut -f1-5 that.vcf) <(cut -f 7 that.bed | cat <(echo "ALT") /dev/stdin) <(cut -f7-) > fixed.vcf
ADD COMMENT
0
Entering edit mode

they don't have identical sort order :(

ADD REPLY
0
Entering edit mode

then sort both files...

ADD REPLY
1
Entering edit mode

the files are of different lengths. I just ended up using a Perl hash, which was much easier.

ADD REPLY

Login before adding your answer.

Traffic: 1919 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