refactor, add mosdepth, and fixes for bwa

14 files changed, 211 insertions, 106 deletions

diff --git a/default.nix b/default.nix
new file mode 100644
index 0000000..a9fb72b
--- /dev/null
+++ b/default.nix
@@ -0,0 +1,14 @@
+{nixpkgs ? import <nixpkgs> {}}:
+
+let
+  bionix = nixpkgs.lib.makeExtensible (self:
+    let callBionix = file: import file { bionix = self; nixpkgs = nixpkgs; };
+    in with self; {
+      bwa = callBionix ./tools/bwa.nix;
+      mosdepth = callBionix ./tools/mosdepth.nix;
+      platypus = callBionix ./tools/platypus.nix;
+      ref = callBionix ./references.nix;
+      samtools = callBionix ./tools/samtools.nix;
+      strelka = callBionix ./tools/strelka.nix;
+  });
+in bionix
diff --git a/references.nix b/references.nix
new file mode 100644
index 0000000..e7f1c7a
--- /dev/null
+++ b/references.nix
@@ -0,0 +1,18 @@
+{ bionix, nixpkgs }:
+
+with nixpkgs;
+
+rec {
+  grch38 = grch38-p12;
+  grch38-p12 = {
+    seq = stdenv.mkDerivation rec {
+        name = "seq-grch38.${version}";
+        version = "p12";
+        src = fetchurl {
+          url = "ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/GRCh38.p12.genome.fa.gz";
+          sha256 = "0ji2ggpmgnbpwbhq8mirj6h3lyy02nl2rnz7n892iq5cqpsblh4z";
+        };
+        buildCommand = "gunzip < $src > $out";
+      };
+  };
+}
diff --git a/test-tnpair.nix b/test-tnpair.nix
index b3a41cb..31570ee 100644
--- a/test-tnpair.nix
+++ b/test-tnpair.nix
@@ -1,23 +1,21 @@
-{pkgs ? import <nixpkgs> {}}:
-
-with pkgs;
+with (import <nixpkgs> {});
+with (import <bionix> {});
 with lib;
 
 let
   ref = ./example/ref.fa;
-  alignWithRG = rg: callPackage ./tools/bwa.nix { inherit ref; flags = "-R'@RG\\tID:${rg}\\tSM:${rg}'";};
-  sort = callPackage ./tools/samtools-sort.nix { };
-  callVariants = callPackage ./tools/strelka.nix { inherit ref; };
+  alignWithRG = rg: bwa.align { inherit ref; flags = "-R'@RG\\tID:${rg}\\tSM:${rg}'";};
+  sort = samtools.sort { };
+  callVariants = strelka.call { inherit ref; };
 
   tnpair = { tumour = {name = "mysample1"; files = {input1 = ./example/sample1-1.fq; input2 = ./example/sample1-2.fq;};};
-               normal = {name = "mysample2"; files = {input1 = ./example/sample2-1.fq; input2 = ./example/sample2-1.fq;};};};
+             normal = {name = "mysample2"; files = {input1 = ./example/sample2-1.fq; input2 = ./example/sample2-1.fq;};};};
 
   processPair = { tumour, normal }: rec {
     alignments = mapAttrs (_: x: sort (alignWithRG x.name x.files)) { inherit normal tumour; };
     variants = callVariants alignments;
   };
 
-  #results = map processPair tnpairs;
   tnpairResult = processPair tnpair;
 
   testNaming = stdenv.mkDerivation {
diff --git a/tools/bwa-index.nix b/tools/bwa-index.nix
index 8ca0eec..b5e3779 100644
--- a/tools/bwa-index.nix
+++ b/tools/bwa-index.nix
@@ -15,5 +15,6 @@ stdenv.mkDerivation {
     bwa index ref.fa
     mkdir $out
     mv ref.fa.* $out
+    grep '^>[^ \t]*_alt$' ref.fa | tr -d '^>' > $out/idxbase.alt || true
   '';
 }
diff --git a/tools/bwa-mem.nix b/tools/bwa-mem.nix
new file mode 100644
index 0000000..5a2772e
--- /dev/null
+++ b/tools/bwa-mem.nix
@@ -0,0 +1,40 @@
+{ stdenv
+, callPackage
+, lib
+, bc
+, bwa
+, samtools ? null
+, ref
+, bamOutput ? true
+, flags ? null
+}:
+
+{ input1
+, input2 ? null
+}:
+
+assert bamOutput -> samtools != null;
+
+with lib;
+
+let index = callPackage ./bwa-index.nix { inherit bwa stdenv lib; } ref;
+
+in stdenv.mkDerivation {
+  name = "bwa-mem";
+  buildInputs = [ bwa bc ] ++ optional bamOutput samtools;
+  buildCommand = ''
+    ln -s ${ref} ref.fa
+    for f in ${index}/* ; do
+      ln -s $f
+    done
+    cores=$(echo $NIX_BUILD_CORES ${optionalString bamOutput "- 1"} | bc)
+    if [[ $cores -lt 1 ]] ; then
+      >&2 echo "not enough build cores"
+      exit 1
+    fi
+    bwa mem ${optionalString (flags != null) flags} -t $cores ref.fa ${input1} \
+      ${optionalString (input2 != null) input2} \
+      ${optionalString bamOutput "| samtools view -b"} \
+      > $out
+  '';
+}
diff --git a/tools/bwa.nix b/tools/bwa.nix
index 20a308f..df61cba 100644
--- a/tools/bwa.nix
+++ b/tools/bwa.nix
@@ -1,34 +1,8 @@
-{ stdenv
-, callPackage
-, lib
-, bwa
-, samtools ? null
-, ref
-, bamOutput ? true
-, flags ? null
-}:
+{ bionix, nixpkgs }:
 
-{ input1
-, input2 ? null
-}:
+with nixpkgs;
 
-assert bamOutput -> samtools != null;
-
-with lib;
-
-let index = callPackage ./bwa-index.nix { inherit bwa stdenv lib; } ref;
-
-in stdenv.mkDerivation {
-  name = "bwa-mem";
-  buildInputs = [ bwa ] ++ optional bamOutput samtools;
-  buildCommand = ''
-    ln -s ${ref} ref.fa
-    for f in ${index}/* ; do
-      ln -s $f
-    done
-    bwa mem ${optionalString (flags != null) flags} -t $NIX_BUILD_CORES ref.fa ${input1} \
-      ${optionalString (input2 != null) input2} \
-      ${optionalString bamOutput "| samtools view -b"} \
-      > $out
-  '';
+{
+  align = attrs: callPackage ./bwa-mem.nix attrs;
+  index = attrs: callPackage ./bwa-index.nix attrs;
 }
diff --git a/tools/mosdepth-depth.nix b/tools/mosdepth-depth.nix
new file mode 100644
index 0000000..3340872
--- /dev/null
+++ b/tools/mosdepth-depth.nix
@@ -0,0 +1,14 @@
+{ stdenv
+, mosdepth
+, flags ? null}:
+
+input:
+
+stdenv.mkDerivation {
+  name = "mosdepth-depth";
+  buildInputs = [ mosdepth ];
+  buildCommand = ''
+    mkdir $out
+    mosdepth -t $NIX_BUILD_CORES ${optionalString (flags != null) flags} $out/out ${input}
+  '';
+}
diff --git a/tools/mosdepth-plot.nix b/tools/mosdepth-plot.nix
new file mode 100644
index 0000000..464e495
--- /dev/null
+++ b/tools/mosdepth-plot.nix
@@ -0,0 +1,15 @@
+{ stdenv
+, mosdepth
+, python
+, flags ? null}:
+
+input:
+
+stdenv.mkDerivation {
+  name = "mosdepth-depth";
+  buildInputs = [ python ];
+  buildCommand = ''
+    python ${mosdepth.src}/scripts/plot-dist.py ${input}/*global.dist.txt
+    mv dist.html $out
+  '';
+}
diff --git a/tools/mosdepth.nix b/tools/mosdepth.nix
new file mode 100644
index 0000000..447e530
--- /dev/null
+++ b/tools/mosdepth.nix
@@ -0,0 +1,8 @@
+{ bionix, nixpkgs }:
+
+with nixpkgs;
+
+{
+  depth = attrs: callPackage ./mosdepth-depth.nix attrs;
+  plot = attrs: callPackage ./mosdepth-plot.nix attrs;
+}
diff --git a/tools/platypus-callVariants.nix b/tools/platypus-callVariants.nix
new file mode 100644
index 0000000..3e150d2
--- /dev/null
+++ b/tools/platypus-callVariants.nix
@@ -0,0 +1,33 @@
+{ stdenv
+, callPackage
+, lib
+, platypus
+, ref
+, index ? callPackage ./samtools-faidx.nix {}
+, bamIndex ? callPackage ./samtools-index.nix {}
+, flags ? null
+}:
+
+inputs:
+
+with lib;
+
+let filename = path: last (splitString "/" path);
+
+in stdenv.mkDerivation {
+  name = "platypus";
+  buildInputs = [ platypus ];
+  buildCommand = ''
+    ln -s ${ref} ref.fa
+    ln -s ${index ref} ref.fa.fai
+    ${concatMapStringsSep "\n" (p: "ln -s ${p} ${filename p}.bam") inputs}
+    ${concatMapStringsSep "\n" (p: "ln -s ${bamIndex p} ${filename p}.bai") inputs}
+    ls -l
+    platypus callVariants \
+      --nCPU=$NIX_BUILD_CORES \
+      --refFile=ref.fa \
+      ${optionalString (flags != null) flags} \
+      -o $out \
+      --bamFiles=${concatMapStringsSep "," (p: "${filename p}.bam") inputs}
+  '';
+}
diff --git a/tools/platypus.nix b/tools/platypus.nix
index 3e150d2..88d88d7 100644
--- a/tools/platypus.nix
+++ b/tools/platypus.nix
@@ -1,33 +1,7 @@
-{ stdenv
-, callPackage
-, lib
-, platypus
-, ref
-, index ? callPackage ./samtools-faidx.nix {}
-, bamIndex ? callPackage ./samtools-index.nix {}
-, flags ? null
-}:
+{ bionix, nixpkgs }:
 
-inputs:
+with nixpkgs;
 
-with lib;
-
-let filename = path: last (splitString "/" path);
-
-in stdenv.mkDerivation {
-  name = "platypus";
-  buildInputs = [ platypus ];
-  buildCommand = ''
-    ln -s ${ref} ref.fa
-    ln -s ${index ref} ref.fa.fai
-    ${concatMapStringsSep "\n" (p: "ln -s ${p} ${filename p}.bam") inputs}
-    ${concatMapStringsSep "\n" (p: "ln -s ${bamIndex p} ${filename p}.bai") inputs}
-    ls -l
-    platypus callVariants \
-      --nCPU=$NIX_BUILD_CORES \
-      --refFile=ref.fa \
-      ${optionalString (flags != null) flags} \
-      -o $out \
-      --bamFiles=${concatMapStringsSep "," (p: "${filename p}.bam") inputs}
-  '';
+{
+  call = attrs: callPackage ./platypus-callVariants.nix attrs;
 }
diff --git a/tools/samtools.nix b/tools/samtools.nix
new file mode 100644
index 0000000..6fdbef5
--- /dev/null
+++ b/tools/samtools.nix
@@ -0,0 +1,9 @@
+{ bionix, nixpkgs }:
+
+with nixpkgs;
+
+{
+  index = attrs: callPackage ./samtools-index.nix attrs;
+  sort = attrs: callPackage ./samtools-sort.nix attrs;
+  faidx = attrs: callPackage ./samtools-faidx.nix attrs;
+}
diff --git a/tools/strelka-call.nix b/tools/strelka-call.nix
new file mode 100644
index 0000000..ec7a764
--- /dev/null
+++ b/tools/strelka-call.nix
@@ -0,0 +1,40 @@
+{ stdenv
+, callPackage
+, lib
+, strelka
+, ref
+, index ? callPackage ./samtools-faidx.nix {}
+, bamIndex ? callPackage ./samtools-index.nix {}
+, flags ? null
+}:
+
+{normal, tumour}:
+
+with lib;
+
+let
+  filename = path: last (splitString "/" path);
+  inputs = [ normal tumour ];
+
+in stdenv.mkDerivation {
+  name = "strelka";
+  buildInputs = [ strelka ];
+  buildCommand = ''
+    ln -s ${ref} ref.fa
+    ln -s ${index ref} ref.fa.fai
+    ${concatMapStringsSep "\n" (p: "ln -s ${p} ${filename p}.bam") inputs}
+    ${concatMapStringsSep "\n" (p: "ln -s ${bamIndex p} ${filename p}.bai") inputs}
+
+    configureStrelkaSomaticWorkflow.py \
+      --normalBam ${filename normal}.bam \
+      --tumourBam ${filename tumour}.bam \
+      --ref ref.fa \
+      --runDir $TMPDIR
+
+    ./runWorkflow.py \
+      -m local \
+      -j $NIX_BUILD_CORES
+
+    cp -r results $out
+   '';
+}
diff --git a/tools/strelka.nix b/tools/strelka.nix
index ec7a764..a115740 100644
--- a/tools/strelka.nix
+++ b/tools/strelka.nix
@@ -1,40 +1,7 @@
-{ stdenv
-, callPackage
-, lib
-, strelka
-, ref
-, index ? callPackage ./samtools-faidx.nix {}
-, bamIndex ? callPackage ./samtools-index.nix {}
-, flags ? null
-}:
+{ bionix, nixpkgs }:
 
-{normal, tumour}:
+with nixpkgs;
 
-with lib;
-
-let
-  filename = path: last (splitString "/" path);
-  inputs = [ normal tumour ];
-
-in stdenv.mkDerivation {
-  name = "strelka";
-  buildInputs = [ strelka ];
-  buildCommand = ''
-    ln -s ${ref} ref.fa
-    ln -s ${index ref} ref.fa.fai
-    ${concatMapStringsSep "\n" (p: "ln -s ${p} ${filename p}.bam") inputs}
-    ${concatMapStringsSep "\n" (p: "ln -s ${bamIndex p} ${filename p}.bai") inputs}
-
-    configureStrelkaSomaticWorkflow.py \
-      --normalBam ${filename normal}.bam \
-      --tumourBam ${filename tumour}.bam \
-      --ref ref.fa \
-      --runDir $TMPDIR
-
-    ./runWorkflow.py \
-      -m local \
-      -j $NIX_BUILD_CORES
-
-    cp -r results $out
-   '';
+{
+  call = attrs: callPackage ./strelka-call.nix attrs;
 }