Merge pull request #12 from brentp/dev

Dev

Author: brentp

.gitignore

@@ -6,4 +6,11 @@ src/slivarpkg/groups
 src/slivarpkg/pracode
 src/slivarpkg/sl_gnotate
 __tmp.groups
+sli.var
+*.zip
+*.bin
+src/slivarpkg/filter
+src/slivarpkg/make_gnotate
+src/slivarpkg/siset
+
 

.travis.yml

@@ -0,0 +1,26 @@
+language: c
+services: docker
+before_install:
+- docker pull brentp/musl-hts-nim
+script:
+- docker run -w /test -v `pwd`:/test brentp/musl-hts-nim scripts/ci-tests.sh
+branches:
+  except:
+  - gh-pages
+before_deploy:
+- git config --local user.name "brentp"
+- git config --local user.email "[email protected]"
+- export TRAVIS_TAG=${TRAVIS_TAG:-$(date +'%Y%m%d%H%M%S')-$(git log --format=%h -1)}
+- git tag $TRAVIS_TAG || true
+deploy:
+  # gem install travis
+  # travis setup releases --force
+  provider: releases
+  api_key:
+    secure: xIC7UCrPPB6lQqx9ZLcUP3OOdwF7NEClgE+HN3u0TfCbDkJvCdatrhezzRp09RhP1snaus56cl4U2oqDjJr+16E1jX/sHMfhF5tinreIRb3tgxlidtGEpiq8ljwK6cLB89Iq/yaihamqlG4Pw/R8AIlMCMyi99yHObTpaCYCptiOygUPkVtpDw5KkWoya9Ol7vrCyTVQK5rENEeivzpfPGtGNKHvNwCxI4f7njLI0VNHP1XLoAieb59ttDEgMIfjWI/T2mwjunIjgQtQKzVbBGUsG9M9eYhc5Jw0xZJQfoRiGn0p+ocxfVNX3H5oVM31sd/EEtuEXPfEO5UKxt5ro42ijd24kNjcVC8rVJxyPdqasjb7IQIXRX8CZqd9SF3boCGAbqSG39obdR2g7BoisL5VmZ/IcTCoj/zE5t9loyVLLJf24YViJwS39l9L51M4R4Oq5jHOtJGdrd/OhTYB9xAtc+DAEKiv1JkDmAJV34gAaZd12SnTCps2BtlttqqiTvXjgxIo62XgrXyBKtNzZ0+toYls2RBKdATBv615+zGQcxnXJkPuvl7DvaNJ2i1xCMNisPTSAFfVSciVulq1/k/IHrMYxY4ceRzrmkfUxzb4STETaNHAzq9cmF8MEa5LSLgSDCB0EqBHyll8JcBDlANBvgHUsUiwUbV0mKYUP6o=
+  skip_cleanup: true
+  draft: true
+  file: "$TRAVIS_BUILD_DIR/slivar_static"
+  on:
+    repo: brentp/slivar
+    tags: true

README.md

@@ -1,9 +1,22 @@
-# slivar: filter/annotate variants in VCF/BCF format with simple expressions
+# slivar: filter/annotate variants in VCF/BCF format with simple expressions [![Build Status](https://travis-ci.org/brentp/slivar.svg?branch=master)](https://travis-ci.org/brentp/slivar)
+
+slivar is a set of command-line tools that enables rapid querying and filtering of VCF files. 
+It facilitates operations on trios and [groups](#groups) and allows arbitrary expressions using simple javascript.
+
+#### use-cases for `slivar`
+
++ annotate variants with combined exomes + whole genomes at > 30K variants/second using only a 1.5GB compressed annotation file
++ call *denovo* variants with a simple expression that uses *mom*, *dad*, *kid* labels that is applied to each trio in a cohort (as inferred from a pedigree file).
+  `kid.alts == 1 && mom.alts == 0 && dad.alts == 0 && kid.DP > 10 && mom.DP > 10 && dad.DP > 10`
++ define and filter on arbitrary groups with labels. For example, 7 sets of samples each with 1 normal and 3 tumor time-points:
+  `normal.AD[0] = 0 && tumor1.AB  < tumor2.AB && tumor2.AB < tumor3.AB`
++ filter variants with simple expressions:
+  `variant.call_rate > 0.9 && variant.FILTER == "PASS" && INFO.AC < 22 && variant.num_hom_alt == 0`
 
 slivar has sub-commands:
 + [expr](#expr): trio and group expressions and filtering
-+ [gnotate](#gnotate): rapidly annotate a VCF/BCF with gnomad
-+ filter: filter a VCF with an expression
++ [gnotate](#gnotate): filter and/or annotate a VCF/BCF files
++ [make-gnotate](#gnotate): make a compressed zip file of annotations for use by slivar
 
 # Table of Contents
 
@@ -47,28 +60,35 @@ will be tested on each of those 200 trios.
 The expressions are javascript so the user can make these as complex as needed.
 
 
-```
+```bash
 slivar expr \
    --pass-only \ # output only variants that pass one of the filters (default is to output all variants)
    --vcf $vcf \
    --ped $ped \
-   --gnomad $gnomad_zip \ # a compressed gnomad zip that allows fast annotation so that `gnomad_af` is available in the expressions below.
-   --load functions.js \ # any valid javascript is allowed here.
+   # 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/functions.js \ 
    --out-vcf annotated.bcf \
-   --info "variant.call_rate > 0.9" \ # this filter is applied before the trio filters and can speed evaluation if it is stringent.
+   # 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.alts == 1 && mom.alts == 0 && dad.alts == 0 \
                    && 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:recessive_func(kid, mom, dad)"
+   --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 `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 `functions.js`; an example implementation of that
+and other useful functions is provided [here](https://github.com/brentp/slivar/blob/master/js/functions.js
+
 #### 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
@@ -91,7 +111,7 @@ sample7	sample8	sample9	sample12
 ```
 
 where `sample10` will be available as "sibling" in the first family and an expression like:
-```
+```bash
 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.
@@ -106,7 +126,7 @@ 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:
-```
+```bash
 normal.alts == 0 && normal.DP > 10 \
   && tumor1.AB > 0 \
   && tumor1.AB < tumor2.AB \
@@ -117,18 +137,49 @@ normal.alts == 0 && normal.DP > 10 \
 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](https://github.com/brentp/slivar/wiki/groups-in-slivar)
+
 ### Gnotate
 
-This uses a compressed, reduced representation of gnomad allele frequencies **and FILTERs** to reduce from the 600+ GB of data for the
-**whole genome and exome** to a 1.5GB file distributed [here](https://s3.amazonaws.com/gemini-annotations/gnomad-2.1.zip).
-The zip file encodes the popmax_AF (whichever is higher between whole genome and exome) and the FILTER for every variant in gnomad.
-It can annotate at faster than 10K variants per second.
+The `gnotate` sub-command allows filtering and/or annotating.
+More extensive documentation and justification for annotating with `gnotate` are [here](https://github.com/brentp/slivar/wiki/gnotate)
+
+`gnotate` uses a compressed, reduced representation of a single value pulled from a (population VCF) along with a boolean that indicates a
+non-pass filter. This can, for example, reduce the 600+ GB of data for the **whole genome and exome** from gnomad to a 1.5GB file
+distributed [here](https://s3.amazonaws.com/gemini-annotations/gnomad-2.1.zip).
+The zip file encodes the popmax_AF (whichever is higher between whole genome and exome) and the presence of FILTER for every variant
+in gnomad.  
+
+It can annotate at faster than 30K variants per second (limited by speed of parsing the query VCF).
+
+```
+slivar gnotate --vcf $input_vcf -o $output_bcf --threads 3 --gnotate encoded.zip
+```
+It's also possible to use `gnotate` as a filtering command without specifying any `--gnotate` arguments.
+
+
+#### make-gnotate
+
+Users can make their own `gnotate` files like:
+
+```bash
+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](https://github.com/brentp/slivar/wiki/make-gnotate)
 
-slivar gnotate --vcf $input_vcf -o $output_bcf --threads 3 -g gnomad-2.1.zip
 
 ## Installation
 
 get the latest binary from: https://github.com/brentp/slivar/releases/latest
+This will require libhts.so (from [htslib](https://htslib.org)) to be in the usual places or in a directory indicated in `LD_LIBRARY_PATH`.
 
 or use via docker from: [brentp/slivar:latest](https://hub.docker.com/r/brentp/slivar)
 

docker/Dockerfile

@@ -30,7 +30,8 @@ RUN cd / && \
     echo 'PATH=/nim/bin:/root/.nimble/bin:$PATH' >> ~/.bashrc && \
     echo 'PATH=/nim/bin:/root/.nimble/bin:$PATH' >> ~/.bash_profile
 
-RUN cd / && export PATH=/nim/bin:/root/.nimble/bin:$PATH && \
+RUN cd /   \
+    && export PATH=/nim/bin:/root/.nimble/bin:$PATH && \
     source scl_source enable devtoolset-2  && \
     nimble install -y nimgen && \
     cd / && git clone -b dev https://github.com/brentp/duktape-nim && \
@@ -40,4 +41,4 @@ RUN cd / && export PATH=/nim/bin:/root/.nimble/bin:$PATH && \
     cd slivar && \
     nimble install --verbose -y && \
     nim c -d:release -o:/usr/bin/slivar --passC:-flto src/slivar && \
-    rm -rf /slivar /duktape-nim && slivar trio -h
+    rm -rf /slivar /duktape-nim && slivar expr -h

js/functions.js

@@ -0,0 +1,25 @@
+function trio_denovo(kid, dad, mom) {
+    // alts are 0:hom_ref, 1:het, 2:hom_alt, -1:unknown
+    if(!(kid.alts == 1 && mom.alts == 0 && dad.alts == 0)){ return false}
+    // sufficient depth in all
+    if(kid.DP < 7 || mom.DP < 7 || dad.DP < 7) { return false; }
+    // no evidence for alternate in the parents
+    if((mom.AD[1] + dad.AD[1]) > 0) { return false; }
+    // check the kid's allele balance.
+    if(kid.AB < 0.2 || kid.AB > 0.8) { return false; }
+    return true
+}
+
+function trio_autosomal_dominant(kid, dad, mom) {
+    // affected samples must be het.
+    if(!(kid.affected == (kid.alts == 1) && mom.affected == (mom.alts == 1) && dad.affected == (dad.alts == 1))) { return false; }
+    if(kid.DP < 7 || mom.DP < 7 || dad.DP < 7) { return false; }
+    return kid.affected
+}
+
+function trio_autosomal_recessive(kid, dad, mom) {
+    if(!(kid.affected == (kid.alts == 2) && mom.affected == (mom.alts == 2) && dad.affected == (dad.alts == 2))) { return false; }
+    if(kid.DP < 7 || mom.DP < 7 || dad.DP < 7) { return false; }
+    if(kid.GQ < 10 || mom.GQ < 10 || dad.GQ < 10) { return false; }
+    return kid.affected
+}

scripts/ci-tests.sh

@@ -0,0 +1,9 @@
+#!/bin/sh
+set -e
+apk add -q bash
+mkdir test-tmp && cd test-tmp
+git clone --depth 1 https://github.com/samtools/htslib
+git clone --depth 1 https://github.com/samtools/bcftools
+cd bcftools && make -j4 install && cd ../.. && rm -rf test-tmp
+nimble install -y && nimble test
+nim c -d:release -o:/test/slivar_static -d:nsb_static src/slivar

slivar.nimble

@@ -18,7 +18,7 @@ license       = "MIT"
 
 # Dependencies
 
-requires "hts >= 0.2.7",  "binaryheap", "zip", "https://github.com/brentp/duktape-nim#dev", "https://github.com/brentp/bpbio", "https://github.com/brentp/nim-minizip"
+requires "hts >= 0.2.7", "nimgen", "binaryheap", "zip", "https://github.com/brentp/duktape-nim#dev", "https://github.com/brentp/bpbio", "https://github.com/brentp/nim-minizip"
 srcDir = "src"
 installExt = @["nim"]
 

src/slivar.nim

@@ -5,7 +5,7 @@ import ./slivarpkg/evaluator
 import ./slivarpkg/groups
 
 import ./slivarpkg/gnotate
-import ./slivarpkg/sl_gnotate
+import ./slivarpkg/make_gnotate
 import ./slivarpkg/filter
 import strutils
 import hts/vcf
@@ -19,7 +19,7 @@ proc expr_main*(dropfirst:bool=false) =
   let doc = """
 slivar -- variant expression for great good
 
-Usage: slivar expr [options --pass-only --out-vcf <path> --vcf <path> --ped <path> --trio=<expression>... --group-expr=<expression>... --info=<expression>]
+Usage: slivar expr [options --pass-only --out-vcf <path> --vcf <path> --ped <path> --trio=<expression>... --group-expr=<expression>... --info=<expression> --gnotate=<path>...]
 
 About:
 
@@ -48,10 +48,10 @@ Options:
   --pass-only                only output variants that pass at least one of the filters [default: false]
   --trio <string>...         an expression applied to each trio where "mom", "dad", "kid" labels are available from trios inferred from
                              a ped file.
-  --group-expr <string>...   expressions applied to the groups defined in the alias option.
+  --group-expr <string>...   expressions applied to the groups defined in the alias option [see: https://github.com/brentp/slivar/wiki/groups-in-slivar].
   --info <string>            a filter expression using only variables from  the info field and variant attributes. If this filter
                              does not pass, the trio and alias expressions will not be applied.
-  -g --gnomad <path>          optional path compressed gnomad allele frequencies distributed at: https://github.com/brentp/slivar/releases
+  -g --gnotate <path>...     optional paths compressed gnote file (made with slivar make-gnotate)
   """
 
   var args: Table[string, docopt.Value]
@@ -74,7 +74,7 @@ Options:
     ivcf:VCF
     ovcf:VCF
     groups: seq[Group]
-    gno:Gnotater
+    gnos:seq[Gnotater]
     samples:seq[Sample]
 
   if not open(ivcf, $args["--vcf"], threads=1):
@@ -96,21 +96,25 @@ Options:
   if $args["--alias"] != "nil":
     groups = parse_groups($args["--alias"], samples)
 
-  if $args["--gnomad"] != "nil":
-    doAssert gno.open($args["--gnomad"])
-    gno.update_header(ivcf)
+  if $args["--gnotate"] != "nil":
+    for p in @(args["--gnotate"]):
+      var gno:Gnotater
+      doAssert gno.open(p)
+      gno.update_header(ivcf)
+      gnos.add(gno)
 
   ovcf.copy_header(ivcf.header)
   var
     trioTbl: TableRef[string,string]
     grpTbl: TableRef[string, string]
+    iTbl: TableRef[string, string]
 
   if $args["--trio"] != "nil":
     trioTbl = ovcf.getExpressionTable(@(args["--trio"]), $args["--vcf"])
   if $args["--group-expr"] != "nil":
     grpTbl = ovcf.getExpressionTable(@(args["--group-expr"]), $args["--vcf"])
   doAssert ovcf.write_header
-  var ev = newEvaluator(samples, groups, trioTbl, grpTbl, $args["--info"], gno, field_names=id2names(ivcf.header))
+  var ev = newEvaluator(samples, groups, iTbl, trioTbl, grpTbl, $args["--info"], gnos, field_names=id2names(ivcf.header))
 
   if $args["--js"] != "nil":
     var js = $readFile($args["--js"])
@@ -157,8 +161,8 @@ proc main*() =
 
   var dispatcher = {
     "expr": pair(f:expr_main, description:"trio and group expressions and filtering"),
-    "gnotate": pair(f:sl_gnotate.main, description:"rapidly annotate a VCF/BCF with gnomad"),
-    "filter": pair(f:filter.main, description:"filter a vcf with javascript expressions"),
+    "gnotate": pair(f:filter.main, description:"filter and/or annotate a VCF/BCF"),
+    "make-gnotate": pair(f:make_gnotate.main, description:"make a gnotate zip file for a given VCF"),
     }.toTable
 
   stderr.write_line "slivar version: " & slivarVersion

src/slivarpkg/duko.nim

@@ -81,6 +81,16 @@ proc check*(d:Dukexpr): bool {.inline.} =
     result = d.ctx.duk_get_boolean(-1)
   d.ctx.pop()
 
+proc asfloat*(d:Dukexpr): float32 {.inline.} =
+  ## evaluate a (previously compiled) expression and return a float32
+  discard d.ctx.duk_push_heapptr(d.vptr)
+  if d.ctx.duk_pcall(0) != 0:
+    var err = $d.ctx.duk_safe_to_string(-1)
+    raise newException(ValueError, "error from duktape: " & $err & " for expression:" & d.expr & "\n")
+  else:
+    result = d.ctx.duk_get_number(-1).float
+  d.ctx.pop()
+
 proc check*(ctx: DTContext, expression: string): bool {.inline.} =
     ## evaluate the expression in the current context
     if ctx.duk_peval_string(expression):