---
author: Sam White
toc-title: Contents
toc-depth: 5
toc-location: left
date: 2018-10-15 13:55:44+00:00
layout: post
slug: vcf-splitting-with-bcftools
title: VCF Splitting with bcftools
categories:
  - 2018
  - Miscellaneous
tags:
  - .vcf
  - bcftools
  - htslib
  - tabix
---

[Steven asked for some help trying to split a VCF in to individual VCF files](https://github.com/RobertsLab/resources/issues/421).

VCF file (15GB): [SNP.TRSdp5g95FnDNAmaf05.vcf.gz](https://gannet.fish.washington.edu/seashell/eog/files/VCF_files/SNP.TRSdp5g95FnDNAmaf05.vcf.gz)

Skip to the Results section if you don't want to read through the tials and tribulations of getting this to work.

Here's an overview of how I managed to get this to work and what didn't work.





  * Installed [bcftools, htslib, and set up the bcftools plugins](https://samtools.github.io/bcftools/)



Figured out the VCF file needed to be sorted, bgzipped (part of htslib), and indexed with [tabix](https://www.htslib.org/doc/tabix.html), due to the following error when initially trying to process with VCF file using bcftools:

`[W::vcf_parse] contig '' is not defined in the header. (Quick workaround: index the file with tabix.)
Undefined tags in the header, cannot proceed in the sample subset mode.`

So, I did that:





  * Sort and bgzip:




    
    <code>cat SNP.TRSdp5g95FnDNAmaf05.vcf | \
    awk '$1 ~ /^#/ {print $0;next} {print $0 | "sort -k1,1 -k2,2n"}' | \
    bgzip --threads 20 > SNP.TRSdp5g95FnDNAmaf05.sorted.vcf.gz</code>







  * Index with tabix:




    
    <code>tabix --preset vcf SNP.TRSdp5g95FnDNAmaf05.sorted.vcf.gz</code>



This produced a separate file:





  * SNP.TRSdp5g95FnDNAmaf05.sorted.vcf.gz.tbi.



It seems as though this file must exist in the same directory as the source VCF for it to be utilized, although no commands work directly with this index file.

Then, tried [biostars solution](https://www.biostars.org/p/130456/#243638), but produces an error


    
    <code>
    #!/bin/bash
    for file in *.vcf.gz; do
      for sample in `bcftools query -l $file`; do
        bcftools view -c1 -Oz -s $sample -o ${file/.vcf*/.$sample.vcf.gz} $file
      done
    done
    </code>



Resulting error:


    
    <code>[E::bcf_calc_ac] Incorrect AN/AC counts at NC_035780.1:26174</code>



And empty split VCF files...

Tried tabix on unsorted bgzipped file yields this error:


    
    <code>[E::hts_idx_push] chromosome blocks not continuous</code>



Tried modified sort:


    
    <code>
    cat SNP.TRSdp5g95FnDNAmaf05.vcf | \
    awk '$1 ~ /^#/ {print $0;next} {print $0 | "sort -k1,1V -k2,2n"}' | \
    bgzip --threads 20 > SNP.TRSdp5g95FnDNAmaf05.sorted.vcf.gz
    </code>



Produces this error:


    
    <code>[E::bcf_calc_ac] Incorrect AN/AC counts at NC_035780.1:26174</code>



And empty split VCF files...

Changed to new version of "view" - trying "call" instead (it seems that bcftools view is deprecated?):


    
    <code>
    #!/bin/bash
    for file in *.vcf.gz; do
      for sample in `bcftools query -l $file`; do
        bcftools call \
        --consensus-caller \
        --output-type z \
        --threads 18 \
        --samples $sample 
        --output-file ${file/.vcf.gz/.$sample.vcf.gz} \
        $file
      done
    done
    </code>



Still results in empty output files.

Based off of the repeated error about AN/AC counts, tried to fill AN/AC values...


    
    <code>
    bcftools plugin fill-AN-AC SNP.TRSdp5g95FnDNAmaf05.sorted.vcf.gz \
    --output-type z \
    --threads 18 \
    --output SNP.TRSdp5g95FnDNAmaf05.sorted.ANACfill.vcf.gz
    </code>



And, ran this code:


    
    <code>
    #!/bin/bash
    for file in SNP.TRSdp5g95FnDNAmaf05.sorted.ANACfill.vcf.gz; do
      for sample in `bcftools query -l $file`; do
        bcftools call \
        --consensus-caller \
        --output-type z \
        --threads 18 \
        --samples $sample 
        --output-file ${file/.vcf.gz/.$sample.vcf.gz} \
        $file
      done
    done
    </code>



Still results in empty files...

Try original code again (expanded shortened arguments to improve readability):


    
    <code>
    #!/bin/bash
    for file in SNP.TRSdp5g95FnDNAmaf05.sorted.ANACfill.vcf.gz; do
      for sample in `bcftools query -l $file`; do
        bcftools view \
        --min-ac 1 \
        --output-type z \
        --samples $sample \
        --output-file ${file/.vcf*/.$sample.vcf.gz} \
        --threads 18 \
        $file
      done
    done
    </code>



P.S. I realize the outermost for loop is not necessary, but it was faster/easier to just quickly modify the code from that Biostars solution.



* * *





# RESULTS



That worked! Or, at least it is producing non-empty, split VCF files! I'll let Steven know and let him decide what impact (if any) the `fill-AN-AC` plugin had on the file(s)!

BTW, running with 18 threads on my computer, this took ~30mins to split each sample into its own VCF file. However, checking performance via `htop`, it certainly does not appear to be using 18 threads...

Output folder: [20181015_vcf_split](https://gannet.fish.washington.edu/Atumefaciens/20181015_vcf_split/)
