search

brentp / slivar

genetic variant expressions, annotation, and filtering for great good.
MIT License
247 stars 23 forks source link
genomics rare-disease rare-variant-analysis variant-analysis variant-interpretation

readme

slivar: filter/annotate variants in VCF/BCF format with simple expressions Build Status

downloads

If you use slivar, please cite the paper

slivar is a set of command-line tools that enables rapid querying and filtering of VCF files. It facilitates operations on trios and groups and allows arbitrary expressions using simple javascript.

use-cases for slivar

slivar logo

slivar has sub-commands:

Table of Contents

Installation

get the latest binary from: https://github.com/brentp/slivar/releases/latest

slivar_static does not depend on any libraries and should work on any 64 bit linux system.

slivar_shared will require libhts.so (from htslib) to be in the usual places or in a directory indicated in LD_LIBRARY_PATH.

or use via docker from: brentp/slivar:latest

QuickStart

To get started quickly, grab a static binary for the latest release and then follow this example

So for hg38:

vcf=/path/to/your/vcf.vcf.gz
ped=/path/to/your/pedigree.ped
wget https://github.com/brentp/slivar/releases/download/v0.2.8/slivar
chmod +x ./slivar
wget https://raw.githubusercontent.com/brentp/slivar/master/js/slivar-functions.js
wget https://slivar.s3.amazonaws.com/gnomad.hg38.genomes.v3.fix.zip

# example command
./slivar expr --js slivar-functions.js -g gnomad.hg38.genomes.v3.fix.zip \
    --vcf $vcf --ped $ped \
    --info "INFO.gnomad_popmax_af < 0.01 && variant.FILTER == 'PASS'" \
    --trio "example_denovo:denovo(kid, dad, mom)" \
    --family-expr "denovo:fam.every(segregating_denovo)" \
    --trio "custom:kid.het && mom.het && dad.het && kid.GQ > 20 && mom.GQ > 20 && dad.GQ > 20" \
    --pass-only

The pedigree format is explained here

Commands

expr

expr allows filtering on (abstracted) trios and groups. For example, given a VCF (and ped/fam file) with 100 trios, slivar will apply an expression with kid, mom, dad identifiers to each trio that it automatically extracts.

expr can also be used, for example to annotate with population allele frequencies from a gnotate file without any sample filtering. See the wiki for more detail and the gnotate section for gnotation files that we distribute for slivar.

expr commands are quite fast, but can be parallelized using pslivar.

trio

when --trio is used, slivar finds all trios in a VCF, PED pair and let's the user specify an expression with indentifiers of kid, mom, dad that is applied to each possible trio. For example, a simple expression to call de novo variants:

variant.FILTER == 'PASS' && \                         # 
variant.call_rate > 0.95 && \                         # genotype must be known for most of cohort.
INFO.gnomad_af < 0.001 && \                           # rare in gnomad (must be in INFO [but see below])
kid.het && mom.hom_ref && dad.hom_ref && \            # also unknown
kid.DP > 7 && mom.DP > 7 && dad.DP > 7 && \           # sufficient depth in all
(mom.AD[1] + dad.AD[1]) == 0                          # no evidence for alternate in the parents

This requires passing variants that are rare in gnomad that have the expected genotypes and do not have any alternate evidence in the parents. If there are 200 trios in the ped::vcf given, then this expression will be tested on each of those 200 trios.

When trios are not sufficient, use Family Expressions which allow more heterogeneous family structures.

The expressions are javascript so the user can make these as complex as needed.

slivar expr \
   --pass-only \ # output only variants that pass one of the filters (default is to output all variants)
   --vcf $vcf \
   --ped $ped \
   # compressed zip that allows fast annotation so that `gnomad_af` is available in the expressions below.
   --gnotate $gnomad_af.zip \ 
   # any valid javascript is allowed in a file here. provide functions to be used below.
   --js js/slivar-functions.js \ 
   --out-vcf annotated.bcf \
   # this filter is applied before the trio filters and can speed evaluation if it is stringent.
   --info "variant.call_rate > 0.9" \ 
   --trio "denovo:kid.het && mom.hom_ref && dad.hom_ref \
                   && kid.AB > 0.25 && kid.AB < 0.75 \
                   && (mom.AD[1] + dad.AD[1]) == 0 \
                   && kid.GQ >= 20 && mom.GQ >= 20 && dad.GQ >= 20 \
                   && kid.DP >= 12 && mom.DP >= 12 && dad.DP >= 12" \
   --trio "informative:kid.GQ > 20 && dad.GQ > 20 && mom.GQ > 20 && kid.alts == 1 && \
           ((mom.alts == 1 && dad.alts == 0) || (mom.alts == 0 && dad.alts == 1))" \
   --trio "recessive:trio_autosomal_recessive(kid, mom, dad)"

Note that slivar does not give direct access to the genotypes, instead exposing hom_ref, het, hom_alt and unknown or via alts where 0 is homozygous reference, 1 is heterozygous, 2 is homozygous alternate and -1 when the genotype is unknown. It is recommended to decompose a VCF before sending to slivar

Here it is assumed that trio_autosomal_recessive is defined in slivar-functions.js; an example implementation of that and other useful functions is provided here. Note that it's often better to use --family-expr instead as it's more flexible than trio expressions.

Family Expressions

Trios are a nice abstraction for cohorts consisting of only trios, but for more general uses, there is --family-expr for example, given either a duo, or a quartet, we can find variants present only in affected samples with:

 --family-expr "aff_only:fam.every(function(s) { return s.het == s.affected && s.hom_ref == !s.affected && s.GQ > 5 })"

Note that this does not explicitly check for transmission or non-transmission between parents and off-spring so it is less transparent than the trio mode, but more flexible.

Groups

A trio is a special-case of a group that can be inferred from a pedigree. For more specialized use-cases, a group can be specified. For example we could, instead of using --trio, use a group file like:

#kid    mom dad
sample1 sample2 sample3
sample4 sample5 sample6
sample7 sample8 sample9

Where, here we have specified 3 trios below a header with their "labels". This can be accomplished using --trio, but we can for example specify quartets like this:

#kid    mom dad sibling
sample1 sample2 sample3 sample10
sample4 sample5 sample6 sample11
sample7 sample8 sample9 sample12

where sample10 will be available as "sibling" in the first family and an expression like:

kid.alts == 1 && mom.alts == 0 && dad.alts == 0 and sibling.alts == 0

could be specified and it would automatically be applied to each of the 3 families.

Another example could be looking at somatic variants with 3 samples, each with a normal and 4 time-points of a tumor:

#normal tumor1  tumor2  tumor3  tumor4
ss1 ss8 ss9 ss10    ss11
ss2 ss12    ss13    ss14    ss15    
ss3 ss16    ss17    ss18    ss19

where, again each row is a sample and the ID's (starting with "ss") will be injected for each sample to allow a single expression like:

normal.hom_ref && normal.DP > 10 \
  && tumor1.AB > 0 \
  && tumor1.AB < tumor2.AB \
  && tumor2.AB < tumor3.AB \
  && tumor3.AB < tumor4.AB

to find a somatic variant that has increasing frequency (AB is allele balance) along the tumor time-points. More detail on groups is provided here

Sample Expressions

Users can specify a boolean expression that is tested against each sample using e.g.:

--sample-expr "hi_quality:sample.DP && sample.GQ > 10"

Each sample that passes this expression will be have its sample id appended to the INFO field of hi_quality which is added to the output VCF.

make-gnotate

Users can make their own gnotate files like:

slivar make-gnotate --prefix gnomad \
    --field AF_popmax:gnomad_popmax_af \
    --field nhomalt:gnomad_num_homalt \
    gnomad.exomes.r2.1.sites.vcf.gz gnomad.genomes.r2.1.sites.vcf.gz

this will pull AF_popmax and nhomalt from the INFO field and put them into gnomad.zip as gnomad_popmax_af and gnomad_num_homalt respectively. The resulting zip file will contain the union of values seen in the exome and genomes files with the maximum value for any intersection. Note that the names (gnomad_popmax_af and gnomad_num_homalt in this case) should be chosen carefully as those will be the names added to the INFO of any file to be annotated with the resulting gnomad.zip

More information on make-gnotate is in the wiki

compound-het

This command is used to find compound heterozygous variants (with phasing-by-inheritance) in trios. It is used after filtering to rare(-ish) heterozygotes.

See a full description of use here

NOTE that by default, this command limits to a subset of impacts; this is adjustable with the --skip flag. See more on the wiki

tsv

This command is used to convert a filtered and annotated VCF to a TSV (tab-separated value file) for final examination. An example use is:

slivar tsv -p $ped \
    -s denovo -s x_recessive \
    -c CSQ \
    -i gnomad_popmax_af -i gnomad_nhomalt \
    -g gene_desc.txt -g clinvar_gene_desc.txt \
    $vcf > final.tsv

where denovo and x_recessive indicate the INFO fields that contain lists of samples (as added by slivar) that should be extracted. and gnomad_popmax_af and gnomad_nhomalt are pulled from the INFO field. The -c arugment (CSQ) tells slivar that it can get gene, transcript and impact information from the CSQ field in the INFO. And the -g arguments are tab-delimited files of gene -> description where the description is added to the text output for quick inspection. Run slivar tsv without any arguments for examples on how to create these for pLI and clinvar.

Also see the wiki

duo-del

slivar duo-del finds structural deletions in parent-child duos using non-transmission of alleles. this can work to find deletions in exome data using genotypes, thereby avoiding the problems associated with depth-based CNV calling in exomes.

see: https://github.com/brentp/slivar/wiki/finding-deletions-in-parent-child-duos

Data Driven Cutoffs

slivar ddc is a tool to discover data-driven cutoffs from a VCF and pedigree information. It generates an interative VCF so a user can see how mendelian violation and transmissions are effected by varying cutoffs for values in the INFO and FORMAT fields.

See the wiki for more details.

Attributes

How it works

slivar embeds the duktape javascript engine to allow the user to specify expressions. For each variant, each trio (and each sample), it fills the appropriate attributes. This can be intensive for VCFs with many samples, but this is done as efficiently as possible such that slivar can evaluate 10's of thousand of variants per second even with dozens of trios.

Summary Table

slivar outputs a summary table with rows of samples and columns of expression where each value indicates the number of variants that passed the expression in each sample. By default, this goes to STDOUT but if the environment variable SLIVAR_SUMMARY_FILE is set, slivar will write the summary to that file instead.

Gnotation Files

Users can create their own gnotation files with slivar make-gnotate, but we provide: