Answers.Body: <p>According to the <code>bcftools</code> man page, it is able to produce statistics using the command <code>bcftools stats</code>. Running this myself, the statistics look like what you're asking for:</p>
  <pre><code># This file was produced by bcftools stats (1.2-187-g1a55e45+htslib-1.2.1-256-ga356746) and can be plotted using plot-vcfstats.
  # The command line was: bcftools stats  OVLNormalised_STARout_KCCG_called.vcf.gz
  #
  # Definition of sets:
  # ID    [2]id   [3]tab-separated file names
  ID  0   OVLNormalised_STARout_KCCG_called.vcf.gz
  # SN, Summary numbers:
  # SN    [2]id   [3]key  [4]value
  SN  0   number of samples:  108
  SN  0   number of records:  333
  SN  0   number of no-ALTs:  0
  SN  0   number of SNPs: 313
  SN  0   number of MNPs: 0
  SN  0   number of indels:   20
  SN  0   number of others:   0
  SN  0   number of multiallelic sites:   0
  SN  0   number of multiallelic SNP sites:   0
  # TSTV, transitions/transversions:
  # TSTV  [2]id   [3]ts   [4]tv   [5]ts/tv    [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
  TSTV    0   302 11  27.45   302 11  27.45
  # SiS, Singleton stats:
  ...
  # IDD, InDel distribution:
  # IDD   [2]id   [3]length (deletions negative)  [4]count
  IDD 0   -9  1
  IDD 0   -2  4
  IDD 0   -1  6
  IDD 0   1   4
  IDD 0   2   1
  IDD 0   3   3
  IDD 0   4   1
  # ST, Substitution types:
  # ST    [2]id   [3]type [4]count
  ST  0   A&gt;C 2
  ST  0   A&gt;G 78
  ST  0   A&gt;T 2
  ST  0   C&gt;A 5
  ST  0   C&gt;G 0
  ST  0   C&gt;T 66
  ST  0   G&gt;A 67
  ST  0   G&gt;C 0
  ST  0   G&gt;T 1
  ST  0   T&gt;A 1
  ST  0   T&gt;C 91
  ST  0   T&gt;G 0
  # DP, Depth distribution
  # DP    [2]id   [3]bin  [4]number of genotypes  [5]fraction of genotypes (%)    [6]number of sites  [7]fraction of sites (%)
  DP  0   &gt;500    0   0.000000    333 100.000000
  </code></pre>
  <p>This is annotating what is explicitly in the VCF file, and that's SNVs (and INDELs). If you want a structural variant analysis (i.e. on a larger scale than single nucleotides), then you'll need to use something that does more than a summary of the VCF file.</p>
  <p>While Inversions, large-scale deletions, and breakpoints are part of the 4.1 specifications, they are unfortunately <a href="https://github.com/samtools/bcftools/issues/623" rel="nofollow noreferrer">not currently supported by bcftools</a>.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-05-29T05:16:36.277Z
Answers.LastActivityDate: 2022-11-04T03:21:17.493Z
Answers.LastEditDate: 2022-11-04T03:21:17.493Z
Answers.UserDisplayName: gringer
Answers.Body: <p>For structural variants  I think you can use <a href="https://github.com/fritzsedlazeck/SURVIVOR" rel="nofollow noreferrer">SURVIVOR</a> like <code>SURVIVOR stats</code> which was design specifically for such purpose (statistics on SV file).</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-16T20:19:02.353Z
Answers.LastActivityDate: 2019-09-16T20:19:02.353Z
Answers.UserDisplayName: Medhat Helmy
Answers.Body: <p>You might want to try <a href="https://github.com/pjotrp/bioruby-vcf" rel="nofollow noreferrer">Bio-VCF</a>. From the authors description </p>
  
  <blockquote>
    <p>Bio-vcf is a new generation VCF parser, filter and converter. Bio-vcf
    is not only very fast for genome-wide (WGS) data, it also comes with a
    really nice filtering, evaluation and rewrite language and it can
    output any type of textual data, including VCF header and contents in
    RDF and JSON.</p>
  </blockquote>
  
  <p>In addition it is a very fast parser. The rewrite DSL might help you customise your filtering and needs. </p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-03T10:55:52.707Z
Answers.Text: Has Bio-vcf got any functionality for structural variation as the OP requested? A quick scan of the GitHub page reveals nothing.  Whilst many of these tools are useful for filtering SNVs and short indels most fail to do anything of worth for SVs.
Answers.UserDisplayName: Matt Bashton
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-29T05:23:32.083Z
Answers.LastActivityDate: 2017-05-29T05:23:32.083Z
Answers.UserDisplayName: eastafri
Answers.Body: <p>Let's say you have an INFO field called <code>SVTYPE</code> to indicate the structure variant type.</p>
  
  <p>Here is how you get the stats easily:</p>
  
  <ol>
  <li>Install <code>vcfstats</code>: <code>pip install vcfstats</code></li>
  <li>Define a macro to extract the <code>SVTYPE</code> information, say in <code>svtype.py</code>:</li>
  </ol>
  
  <pre class="lang-py prettyprint-override"><code>from vcfstats.macros import cat
  @cat
  def SVTYPE(variant):
      return variant.INFO['SVTYPE']
  </code></pre>
  
  <ol start="3">
  <li>Generate the statistics:</li>
  </ol>
  
  <pre class="lang-sh prettyprint-override"><code>vcfstats --vcf &lt;your vcf file&gt; \
     --outdir &lt;the output directory&gt; \
     --macro svtype.py \
     --formula 'COUNT(1) ~ SVTYPE[DEL,INS,INV...]' \ # desired svtypes
     --title 'Number of SV types'
  </code></pre>
  
  <p>Then in the output directory you will have the stats file named <code>Number_of_SV_Types.txt</code> and a plot for it as well.</p>
  
  <p>Check the details at <a href="https://github.com/pwwang/vcfstats" rel="nofollow noreferrer">https://github.com/pwwang/vcfstats</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-10-24T08:59:00.560Z
Answers.Text: The advantage of this solution is that it using a standard parser. The disadvantage is that it requires several steps (unlike `zcat var.vcf.gz | perl -ne 'print "$1\n" if /[;\t]SVTYPE=([^;\t]+)/' | sort | uniq -c
  ` that will make the summary in one go).
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-10-21T23:01:05.413Z
Answers.LastActivityDate: 2019-10-21T23:01:05.413Z
Answers.UserDisplayName: Panwen Wang
Answers.Body: <p>The &quot;my own parser&quot; solutions. The information I was searching for in part of column <code>INFO</code>, namely in variables <code>SVLEN</code> and <code>SVTYPE</code>.</p>
  <p>very fast SV types + counts (by @user172818 in comment) :</p>
  <pre><code>zcat var.vcf.gz | perl -ne 'print &quot;$1\n&quot; if /[;\t]SVTYPE=([^;\t]+)/' | sort | uniq -c
  </code></pre>
  <p>quite slow SV types + counts + sizes :</p>
  <pre><code>SV_colnames &lt;- c('CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'SAMPLE1')
  
  ssplit &lt;- function(s, split = '='){
      unlist(strsplit(s, split = split))
  }
  
  # note, capital letters just respect original naming conventions of the VCF file
  getSVTYPE &lt;- function(info){
      ssplit(grep(&quot;SVTYPE&quot;, info, value = T))[2]
  }
  
  getSVLEN &lt;- function(info){
      SVLEN &lt;- ssplit(grep(&quot;SVLEN&quot;, info, value = T))
      ifelse(length(SVLEN) == 0, NA, as.numeric(SVLEN[2]))
  }
  
  load_sv &lt;- function(file){
      vcf_sv_table &lt;- read.table(file, stringsAsFactors = F)
      colnames(vcf_sv_table) &lt;- SV_colnames
      # possible filtering
      # vcf_sv_table &lt;- vcf_sv_table[vcf_sv_table<span class="math-container">$FILTER == 'PASS',]
      vcf_sv_table_info &lt;- strsplit(vcf_sv_table$</span>INFO, ';')
      vcf_sv_table<span class="math-container">$SVTYPE &lt;- unlist(lapply(vcf_sv_table_info, getSVTYPE))
      vcf_sv_table$</span>SVLEN &lt;- unlist(lapply(vcf_sv_table_info, getSVLEN))
      return(vcf_sv_table)
  }
  
  sv_df &lt;- load_sv('my_sv_calls.vcf')
  table(sv_df$SVTYPE)
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-05-29T06:22:59.417Z
Answers.LastActivityDate: 2022-11-04T03:23:43.010Z
Answers.LastEditDate: 2022-11-04T03:23:43.010Z
Answers.UserDisplayName: Kamil S Jaron
Body: <p>I got a bunch of vcf files (v4.1) with structural variations of bunch of non-model organisms (i.e. there are no known variants). I found there are quite a some tools to manipulate vcf files like <a href="https://vcftools.github.io/" rel="noreferrer">VCFtools</a>, R package <a href="https://cran.r-project.org/web/packages/vcfR/index.html" rel="noreferrer">vcfR</a> or python library <a href="https://pypi.python.org/pypi/PyVCF" rel="noreferrer">PyVCF</a>. However none of them seems to provide a quick summary, something like (preferably categorised by size as well):</p>
  
  <pre><code>type    count
  DEL     x
  INS     y
  INV     z
  ....
  </code></pre>
  
  <p>Is there any tool or a function I overlooked that produces summaries of this style?</p>
  
  <p>I know that vcf file is just a plain text file and if I will dissect <code>REF</code> and <code>ALT</code> columns I should be able to write a script that will do the job, but I hoped that I could avoid to write my own parser.</p>
  
  <p>--- edit ---</p>
  
  <p>So far it seems that only tool that aims to do summaries (@gringer answer) is not working on vcf v4.1. Other tools would provide just partial solution by filtering certain variant type. Therefore I accept my own parser perl/R solutions, till there will be a <strong>working</strong> tool for <strong>stats</strong> of vcf with <strong>structural variants</strong>.</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-28T21:49:35.910Z
Comments.Text: What sort of summary are you looking for? Just raw counts (of number of insertions/deletions)? I'm just curious as I think personally I'd be interested to know how that variation is actually spread along a genome/contig/sequence spatially, would that be of interest to you?
Comments.UserDisplayName: Sam Nicholls
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-28T21:53:21.953Z
Comments.Text: I am comparing couple of relative species. However, I have not identified homologous regions yet (the annotation is so far under construction), therefore positioning of SVs is right now not very relevant. However, what can I do already is take a look if there is at least a difference in counts / sizes of SVs between species.
Comments.UserDisplayName: Kamil S Jaron
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-28T21:57:06.967Z
Comments.Text: Maybe even better than summary would be to get type of SV and its size for every SV call. Then I could a summary in R would be trivial.
Comments.UserDisplayName: Kamil S Jaron
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-28T22:05:03.170Z
Comments.Text: You've probably seen this already, but `vcftools` has a sibling called `bcftools`, which has [a query function, that allows users to query a VCF/BCF to pull out fields and information and output their own format](https://samtools.github.io/bcftools/bcftools-man.html#query). It might not do exactly what you want, but might get you very close (enough to maybe just need a little post-processing in R?).
Comments.UserDisplayName: Sam Nicholls
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-28T22:13:08.293Z
Comments.Text: Great, that sounds promising. I will take a look.
Comments.UserDisplayName: Kamil S Jaron
Tags: <vcf><structural-variation>
Title: Is there an easy way to create a summary of a VCF file (v4.1) with structural variations?
ThreadText: # Is there an easy way to create a summary of a VCF file (v4.1) with structural variations?
  
  (Asked by: username_0 on 2017-05-28)
  
  I got a bunch of vcf files (v4.1) with structural variations of bunch of non-model organisms (i.e. there are no known variants). I found there are quite a some tools to manipulate vcf files like [VCFtools](https://vcftools.github.io/), R package [vcfR](https://cran.r-project.org/web/packages/vcfR/index.html) or python library [PyVCF](https://pypi.python.org/pypi/PyVCF). However none of them seems to provide a quick summary, something like (preferably categorised by size as well):
  
  ```
  type    count
  DEL     x
  INS     y
  INV     z
  ....
  
  ```
  
  Is there any tool or a function I overlooked that produces summaries of this style?
  
  I know that vcf file is just a plain text file and if I will dissect `REF` and `ALT` columns I should be able to write a script that will do the job, but I hoped that I could avoid to write my own parser.
  
  --- edit ---
  
  So far it seems that only tool that aims to do summaries (@username_2 answer) is not working on vcf v4.1. Other tools would provide just partial solution by filtering certain variant type. Therefore I accept my own parser perl/R solutions, till there will be a **working** tool for **stats** of vcf with **structural variants**.
  
  **Tags:** vcf, structural-variation
  
  > **Comments:**
  > **username_1**: What sort of summary are you looking for? Just raw counts (of number of insertions/deletions)? I'm just curious as I think personally I'd be interested to know how that variation is actually spread along a genome/contig/sequence spatially, would that be of interest to you?  
  > **username_0**: I am comparing couple of relative species. However, I have not identified homologous regions yet (the annotation is so far under construction), therefore positioning of SVs is right now not very relevant. However, what can I do already is take a look if there is at least a difference in counts / sizes of SVs between species.  
  > **username_0**: Maybe even better than summary would be to get type of SV and its size for every SV call. Then I could a summary in R would be trivial.  
  > **username_1**: You've probably seen this already, but `vcftools` has a sibling called `bcftools`, which has [a query function, that allows users to query a VCF/BCF to pull out fields and information and output their own format](https://samtools.github.io/bcftools/bcftools-man.html#query). It might not do exactly what you want, but might get you very close (enough to maybe just need a little post-processing in R?).  
  > **username_0**: Great, that sounds promising. I will take a look.  
  
  ---
  
  ## 5 answers
  
  ### Answer by username_0 on 2017-05-29. Score: 6 | Accepted answer
  
  The "my own parser" solutions. The information I was searching for in part of column `INFO`, namely in variables `SVLEN` and `SVTYPE`.
  
  very fast SV types + counts (by @user172818 in comment) :
  
  ```
  zcat var.vcf.gz | perl -ne 'print "$1\n" if /[;\t]SVTYPE=([^;\t]+)/' | sort | uniq -c
  
  ```
  
  quite slow SV types + counts + sizes :
  
  ```
  SV_colnames <- c('CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'SAMPLE1')
  
  ssplit <- function(s, split = '='){
      unlist(strsplit(s, split = split))
  }
  
  # note, capital letters just respect original naming conventions of the VCF file
  getSVTYPE <- function(info){
      ssplit(grep("SVTYPE", info, value = T))[2]
  }
  
  getSVLEN <- function(info){
      SVLEN <- ssplit(grep("SVLEN", info, value = T))
      ifelse(length(SVLEN) == 0, NA, as.numeric(SVLEN[2]))
  }
  
  load_sv <- function(file){
      vcf_sv_table <- read.table(file, stringsAsFactors = F)
      colnames(vcf_sv_table) <- SV_colnames
      # possible filtering
      # vcf_sv_table <- vcf_sv_table[vcf_sv_table$FILTER == 'PASS',]
      vcf_sv_table_info <- strsplit(vcf_sv_table$INFO, ';')
      vcf_sv_table$SVTYPE <- unlist(lapply(vcf_sv_table_info, getSVTYPE))
      vcf_sv_table$SVLEN <- unlist(lapply(vcf_sv_table_info, getSVLEN))
      return(vcf_sv_table)
  }
  
  sv_df <- load_sv('my_sv_calls.vcf')
  table(sv_df$SVTYPE)
  
  ```
  
  
  ---
  
  ### Answer by username_2 on 2017-05-29. Score: 9
  
  According to the `bcftools` man page, it is able to produce statistics using the command `bcftools stats`. Running this myself, the statistics look like what you're asking for:
  
  ```
  # This file was produced by bcftools stats (1.2-187-g1a55e45+htslib-1.2.1-256-ga356746) and can be plotted using plot-vcfstats.
  # The command line was: bcftools stats  OVLNormalised_STARout_KCCG_called.vcf.gz
  #
  # Definition of sets:
  # ID    [2]id   [3]tab-separated file names
  ID  0   OVLNormalised_STARout_KCCG_called.vcf.gz
  # SN, Summary numbers:
  # SN    [2]id   [3]key  [4]value
  SN  0   number of samples:  108
  SN  0   number of records:  333
  SN  0   number of no-ALTs:  0
  SN  0   number of SNPs: 313
  SN  0   number of MNPs: 0
  SN  0   number of indels:   20
  SN  0   number of others:   0
  SN  0   number of multiallelic sites:   0
  SN  0   number of multiallelic SNP sites:   0
  # TSTV, transitions/transversions:
  # TSTV  [2]id   [3]ts   [4]tv   [5]ts/tv    [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
  TSTV    0   302 11  27.45   302 11  27.45
  # SiS, Singleton stats:
  ...
  # IDD, InDel distribution:
  # IDD   [2]id   [3]length (deletions negative)  [4]count
  IDD 0   -9  1
  IDD 0   -2  4
  IDD 0   -1  6
  IDD 0   1   4
  IDD 0   2   1
  IDD 0   3   3
  IDD 0   4   1
  # ST, Substitution types:
  # ST    [2]id   [3]type [4]count
  ST  0   A>C 2
  ST  0   A>G 78
  ST  0   A>T 2
  ST  0   C>A 5
  ST  0   C>G 0
  ST  0   C>T 66
  ST  0   G>A 67
  ST  0   G>C 0
  ST  0   G>T 1
  ST  0   T>A 1
  ST  0   T>C 91
  ST  0   T>G 0
  # DP, Depth distribution
  # DP    [2]id   [3]bin  [4]number of genotypes  [5]fraction of genotypes (%)    [6]number of sites  [7]fraction of sites (%)
  DP  0   >500    0   0.000000    333 100.000000
  
  ```
  
  This is annotating what is explicitly in the VCF file, and that's SNVs (and INDELs). If you want a structural variant analysis (i.e. on a larger scale than single nucleotides), then you'll need to use something that does more than a summary of the VCF file.
  
  While Inversions, large-scale deletions, and breakpoints are part of the 4.1 specifications, they are unfortunately [not currently supported by bcftools](https://github.com/samtools/bcftools/issues/623).
  
  
  ---
  
  ### Answer by username_3 on 2019-09-16. Score: 2
  
  For structural variants I think you can use [SURVIVOR](https://github.com/fritzsedlazeck/SURVIVOR) like `SURVIVOR stats` which was design specifically for such purpose (statistics on SV file).
  
  
  ---
  
  ### Answer by username_4 on 2017-05-29. Score: 2
  
  You might want to try [Bio-VCF](https://github.com/pjotrp/bioruby-vcf). From the authors description
  
  > Bio-vcf is a new generation VCF parser, filter and converter. Bio-vcf
  > is not only very fast for genome-wide (WGS) data, it also comes with a
  > really nice filtering, evaluation and rewrite language and it can
  > output any type of textual data, including VCF header and contents in
  > RDF and JSON.
  
  In addition it is a very fast parser. The rewrite DSL might help you customise your filtering and needs.
  
  > **Comments:**
  > **username_5**: Has Bio-vcf got any functionality for structural variation as the OP requested? A quick scan of the GitHub page reveals nothing.  Whilst many of these tools are useful for filtering SNVs and short indels most fail to do anything of worth for SVs.  
  
  ---
  
  ### Answer by username_6 on 2019-10-21. Score: 1
  
  Let's say you have an INFO field called `SVTYPE` to indicate the structure variant type.
  
  Here is how you get the stats easily:
  
  1. Install `vcfstats`: `pip install vcfstats`
  2. Define a macro to extract the `SVTYPE` information, say in `svtype.py`:
  
  ```py
  from vcfstats.macros import cat
  @cat
  def SVTYPE(variant):
      return variant.INFO['SVTYPE']
  
  ```
  
  3. Generate the statistics:
  
  ```sh
  vcfstats --vcf <your vcf file> \
     --outdir <the output directory> \
     --macro svtype.py \
     --formula 'COUNT(1) ~ SVTYPE[DEL,INS,INV...]' \ # desired svtypes
     --title 'Number of SV types'
  
  ```
  
  Then in the output directory you will have the stats file named `Number_of_SV_Types.txt` and a plot for it as well.
  
  Check the details at <https://github.com/pwwang/vcfstats>
  
  > **Comments:**
  > **username_0**: The advantage of this solution is that it using a standard parser. The disadvantage is that it requires several steps (unlike `zcat var.vcf.gz | perl -ne 'print "$1\n" if /[;\t]SVTYPE=([^;\t]+)/' | sort | uniq -c ` that will make the summary in one go).  
  
  ---
  
  