Repository: vierstralab/footprint-tools (language: python)

Issue:
Pull in release v1.2.1 changes



Discussion / hints:

Fix patch:
diff --git a/.bumpversion.cfg b/.bumpversion.cfg
new file mode 100644
index 0000000..cc625d3
--- /dev/null
+++ b/.bumpversion.cfg
@@ -0,0 +1,6 @@
+[bumpversion]
+current_version = 1.3.0
+commit = True
+tag = True
+
+[bumpversion:file:setup.py]
diff --git a/README.md b/README.md
index 9200c7d..ae8418e 100644
--- a/README.md
+++ b/README.md
@@ -7,16 +7,7 @@ are computed from a series emperically fit negative binomial distribution.
 
 ## Requirements
 
-``footprint-tools`` requires Python 3.6+ and depends on the following additional
-packages:
-
-* cython
-* numpy
-* scipy
-* statsmodels
-* pysam
-* genome_tools (http://www.github.com/jvierstra/genome_tools)
-* pwlf (https://github.com/cjekel/piecewise_linear_fit_py)
+``footprint-tools`` requires Python 3.6+
 
 We also recommend these non-Python analysis tools:
 
@@ -47,7 +38,8 @@ conda env create -f conda-env.yml
 conda activate footprint-tools
 
 # Run commands
-ftd-learn-dispersion... 
+ftd --version
+ftd {commands}
 ```
 
 ## Documentation & usage
diff --git a/conda-env.yml b/conda-env.yml
index d3c69a6..4fc1dfe 100644
--- a/conda-env.yml
+++ b/conda-env.yml
@@ -13,10 +13,14 @@ dependencies:
 - pandas
 - statsmodels
 - simplejson
-- pysam
 - tqdm
+- click
+- matplotlib
 
 - pip:
-  - pwlf
-  - genome_tools
-  - .
+- git+https://github.com/jvierstra/pysam.git
+- pwlf
+- genome_tools
+- yapsin
+- click-option-group
+- .
diff --git a/docs/source/api.rst b/docs/source/api.rst
index 876543b..d3958c9 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -44,7 +44,6 @@ The following code example retrieves reads overlapping a SNV resolved to each al
 
 	import pysam
 
-
 	vcf_file = pysam.VariantFile("variants.vcf.gz")
 	var = vcf_file.fetch(chrom, start, end).next()
 
diff --git a/examples/compute_deviation.slurm b/examples/compute_deviation.slurm
deleted file mode 100644
index dd025f9..0000000
--- a/examples/compute_deviation.slurm
+++ /dev/null
@@ -1,295 +0,0 @@
-#!/bin/bash
-
-# 2018 Jeff Vierstra
-# jvierstra@altius.org
-
-set -e -o pipefail
-
-######## CHANGE VARIABLES BELOW ######## 
-
-bias_model_file=/home/jvierstra/proj/dnase-perspective/cleavage_model/vierstra_et_al.txt
-
-prefix=$1 # Dataset name
-fasta_file=$2 # Genome FASTA file w/ index (.fai)
-bam_file=$3 # Alignment file w/ index (.bai)
-interval_file=$4 # File (BED-format) containing regions where footprint detection will occur (i.e., hotspots)
-output_dir=$5 # Path to directory where output will be stored (created by script if non-existent)
-
-# Change this if a standard UCSC chrom.sizes file doesn't follow this naming convention
-chrom_sizes_file="${fasta_file%.fa}.chrom.sizes"
-
-
-# Size of chunks when parallelzing (=# of cluster jobs)
-chunksize=10000
-
-# FPR thresholds to call footprints
-fps_thresh=(0.2 0.1 0.05 0.01)
-
-
-echo $prefix
-echo $fasta_file
-echo $bam_file
-echo $interval_file
-echo $output_dir
-
-######## DO NOT CHANGE BELOW THIS LINE ######## 
-
-mkdir -p ${output_dir}/logs
-
-# load modules
-module load bedops/2.4.3
-module load python/2.7.11
-module load gcc/5.3.0
-
-###
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.learn_dispersion_model
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --mem=4G
-#SBATCH --cpus-per-task=8
-#SBATCH --partition=queue1 
-
-set -e -o pipefail
-
-TMPDIR=/tmp/\$SLURM_JOB_ID
-mkdir -p \${TMPDIR}
-
-zcat -f ${interval_file} > \${TMPDIR}/intervals.bed
-
-ftd-learn-dispersion-model \
-	--bm ${bias_model_file} --half-win-width 5 \
-	--processors \${SLURM_CPUS_PER_TASK} \
-	${bam_file} ${fasta_file} \${TMPDIR}/intervals.bed \
-> ${output_dir}/dm.json
-
-rm -rf \${TMPDIR}
-__SCRIPT__
-
-###
-
-JOB0=$(sbatch --export=ALL \
-	--job-name=${prefix}.learn_dm \
-	${output_dir}/slurm.learn_dispersion_model)
-
-echo $JOB0
-
-###
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_deviation_chunk
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%A.%a.out
-#SBATCH --mem=8G
-#SBATCH --cpus-per-task=4
-#SBATCH --partition=queue1
-
-set -e -o pipefail
-
-INPUT_FILES=(\`cat ${output_dir}/slurm.compute_deviation_chunk.params | head -n \${SLURM_ARRAY_TASK_ID} | tail -n 1\`)
-
-ftd-compute-deviation \
-	--bm ${bias_model_file} --half-win-width 5 \
-	--smooth-half-win-width 50 --smooth-clip 0.01 \
-	--dm ${output_dir}/dm.json --fdr-shuffle-n 50 \
-	--processors \${SLURM_CPUS_PER_TASK} \
-	${bam_file} ${fasta_file} \${INPUT_FILES[0]} \
-| sort --buffer-size=8G -k1,1 -k2,2n > \${INPUT_FILES[0]}.out
-__SCRIPT__
-
-###
-
-zcat -f ${interval_file} | split -l ${chunksize} -a 4 -d - ${output_dir}/interval.
-ls ${output_dir}/interval.* > ${output_dir}/slurm.compute_deviation_chunk.params
-
-JOB1=$(sbatch --export=ALL \
-	--job-name=${prefix}.compute_deviation \
-	--array=1-$(wc -l < ${output_dir}/slurm.compute_deviation_chunk.params) \
-	--depend=afterok:${JOB0##* } \
-	${output_dir}/slurm.compute_deviation_chunk)
-
-echo $JOB1
-
-###
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_deviation_merge_sort
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --mem=32G
-#SBATCH --cpus-per-task=1
-#SBATCH --partition=queue0
-
-set -e -o pipefail
-
-INPUT_FILES=(\`cat ${output_dir}/slurm.compute_deviation_chunk.params\`)
-OUTPUT_FILES=("\${INPUT_FILES[@]/%/.out}")
-
-sort -k1,1 -k2,2n -S 32G -m \${OUTPUT_FILES[@]} > ${output_dir}/interval.all.bedgraph
-__SCRIPT__
-
-###
-
-JOB2=$(sbatch --export=ALL \
-	--job-name=${prefix}.merge_sort \
-	--depend=afterok:${JOB1##* } \
-	${output_dir}/slurm.compute_deviation_merge_sort)
-
-echo $JOB2
-
-###
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_deviation_working_tracks
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --mem=32G
-#SBATCH --cpus-per-task=1
-#SBATCH --partition=queue0
-
-set -e -o pipefail
-
-TMPDIR=/tmp/\$SLURM_JOB_ID
-mkdir -p \${TMPDIR}
-
-sort-bed --tmpdir \${TMPDIR} --max-mem 32G ${output_dir}/interval.all.bedgraph | starch - > ${output_dir}/interval.all.bedgraph.starch
-
-bgzip -c ${output_dir}/interval.all.bedgraph > ${output_dir}/interval.all.bedgraph.gz
-tabix -0 -p bed ${output_dir}/interval.all.bedgraph.gz
-
-rm -rf \${TMPDIR}
-__SCRIPT__
-
-# Make footprint calls
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_deviation_footprints
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --mem=8G
-#SBATCH --cpus-per-task=1
-#SBATCH --partition=queue1 
-
-set -e -o pipefail
-
-TMPDIR=/tmp/\$SLURM_JOB_ID
-mkdir -p \${TMPDIR}
-
-thresholds=(${fps_thresh[@]})
-thresh=\${thresholds[\${SLURM_ARRAY_TASK_ID}-1]}
-
-cat ${output_dir}/interval.all.bedgraph \
-	| awk -v OFS="\t" -v thresh="\${thresh}" '\$8 <= thresh { print \$1, \$2-3, \$3+3; }' \
-	| sort-bed --tmpdir \${TMPDIR} --max-mem 8G - \
-	| bedops -m - \
-	| awk -v OFS="\t" -v thresh="\${thresh}" '{ \$4="."; \$5=thresh; print; }' \
-> ${output_dir}/interval.all.fps.\${thresh}.bed
-
-bgzip -c ${output_dir}/interval.all.fps.\${thresh}.bed > ${output_dir}/interval.all.fps.\${thresh}.bed.gz
-tabix -0 -p bed ${output_dir}/interval.all.fps.\${thresh}.bed.gz
-
-rm -rf \${TMPDIR}
-__SCRIPT__
-
-# Make UCSC Browser tracks
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_deviation_browser_tracks
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --mem=8G
-#SBATCH --cpus-per-task=1
-#SBATCH --partition=queue0
-
-set -e -o pipefail
-
-TMPDIR=/tmp/\$SLURM_JOB_ID
-mkdir -p \${TMPDIR}
-
-cut -f1-3,4 ${output_dir}/interval.all.bedgraph | awk '\$4 > 0 { print; }' > \${TMPDIR}/interval.all.exp.bedgraph
-bedGraphToBigWig \${TMPDIR}/interval.all.exp.bedgraph ${chrom_sizes_file} ${output_dir}/interval.all.exp.bw
-
-cut -f1-3,5 ${output_dir}/interval.all.bedgraph | awk '\$4 > 0 { print; }' > \${TMPDIR}/interval.all.obs.bedgraph
-bedGraphToBigWig \${TMPDIR}/interval.all.obs.bedgraph ${chrom_sizes_file} ${output_dir}/interval.all.obs.bw
-
-cut -f1-3,6 ${output_dir}/interval.all.bedgraph > \${TMPDIR}/interval.all.lnpval.bedgraph
-bedGraphToBigWig \${TMPDIR}/interval.all.lnpval.bedgraph ${chrom_sizes_file} ${output_dir}/interval.all.lnpval.bw
-
-cut -f1-3,7 ${output_dir}/interval.all.bedgraph > \${TMPDIR}/interval.all.winlnpval.bedgraph
-bedGraphToBigWig \${TMPDIR}/interval.all.winlnpval.bedgraph ${chrom_sizes_file} ${output_dir}/interval.all.winlnpval.bw
-
-cut -f1-3,8 ${output_dir}/interval.all.bedgraph > \${TMPDIR}/interval.all.fpr.bedgraph
-bedGraphToBigWig \${TMPDIR}/interval.all.fpr.bedgraph ${chrom_sizes_file} ${output_dir}/interval.all.fpr.bw
-
-rm -rf \${TMPDIR}
-__SCRIPT__
-
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_deviation_cleanup
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --cpus-per-task=1
-#SBATCH --partition=queue1 
-
-set -e -o pipefail
-
-INPUT_FILES=(\`cat ${output_dir}/slurm.compute_deviation_chunk.params\`)
-OUTPUT_FILES=("\${INPUT_FILES[@]/%/.out}")
-
-rm -f \${INPUT_FILES[@]} \${OUTPUT_FILES[@]}
-
-rm -f ${output_dir}/interval.all.bedgraph
-__SCRIPT__
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_deviation_beta_prior
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --cpus-per-task=1
-#SBATCH --partition=queue1 
-
-ftd-learn-beta-prior \
-	--fdr-cutoff 0.05 --exp-cutoff 10 \
-	--processors \${SLURM_CPUS_PER_TASK} \
-	${output_dir}/interval.all.bedgraph.gz \
-> ${output_dir}/beta_prior.txt
-__SCRIPT__
-
-###
-
-JOB3=$(sbatch --export=ALL \
-	--job-name=${prefix}.working_tracks \
-	--depend=afterok:${JOB2##* } \
-	${output_dir}/slurm.compute_deviation_working_tracks)
-
-echo $JOB3
-
-JOB4=$(sbatch --export=ALL \
-	--job-name=${prefix}.footprints \
-	--array=1-${#fps_thresh[@]} \
-	--depend=afterok:${JOB2##* } \
-	${output_dir}/slurm.compute_deviation_footprints)
-
-echo $JOB4
-
-JOB5=$(sbatch --export=ALL \
-	--job-name=${prefix}.browser_tracks \
-	--depend=afterok:${JOB2##* } \
-	--partition=queue0 \
-	${output_dir}/slurm.compute_deviation_browser_tracks)
-
-echo $JOB5
-
-JOB6=$(sbatch --export=ALL \
-	--job-name=${prefix}.cleanup \
-	--depend=afterok:${JOB3##* }:${JOB4##* }:${JOB5##* } \
-	${output_dir}/slurm.compute_deviation_cleanup)
-
-echo $JOB6
-
-JOB7=$(sbatch --export=ALL \
-	--job-name=${prefix}.beta_prior \
-	--depend=afterok:${JOB3##* } \
-	${output_dir}/slurm.compute_deviation_beta_prior)
-
-echo $JOB7
diff --git a/examples/compute_posterior.slurm b/examples/compute_posterior.slurm
deleted file mode 100644
index 9092da1..0000000
--- a/examples/compute_posterior.slurm
+++ /dev/null
@@ -1,75 +0,0 @@
-#!/bin/bash
-
-# 2018 Jeff Vierstra
-# jvierstra@altius.org
-
-set -o pipefail
-
-metadata_file=$1 # Metadata file
-interval_file=$2 # File (BED-format) containing regions where footprint detection will occur (i.e., hotspots)
-output_dir=$3 # Path to directory where output will be stored (created by script if non-existent)
-
-chunksize=2500
-
-prefix="fps"
-
-####
-
-mkdir -p ${output_dir}/logs
-
-# load modules
-module load bedops/2.4.3
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_posterior_chunk
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --mem=32G
-#SBATCH --cpus-per-task=8
-
-set -e -o pipefail
-
-INPUT_FILES=(\`cat ${output_dir}/slurm.compute_posterior_chunk.params | head -n \${SLURM_ARRAY_TASK_ID} | tail -n 1\`)
-
-ftd-compute-posterior \
-	${base_dir}/metadata.txt \${INPUT_FILES[0]} \
-| sort --buffer-size=16G -k1,1 -k2,2n > \${INPUT_FILES[0]}.out
-__SCRIPT__
-
-cat ${interval_file} | split -l ${chunksize} -a 4 -d - ${output_dir}/interval.
-ls ${output_dir}/interval.* > ${output_dir}/slurm.compute_posterior_chunk.params
-
-JOB0=$(sbatch --export=ALL \
-	--job-name=${prefix}.compute_posterior_chunk \
-	--array=1-$(wc -l < ${output_dir}/slurm.compute_posterior_chunk.params) \
-	${output_dir}/slurm.compute_posterior_chunk)
-
-echo $JOB0
-
-cat <<__SCRIPT__ > ${output_dir}/slurm.compute_posterior_merge_sort
-#!/bin/bash
-#
-#SBATCH --output=${output_dir}/logs/%J.out
-#SBATCH --mem=32G
-#SBATCH --cpus-per-task=1
-
-set -e -o pipefail
-
-TMPDIR=/tmp/\$SLURM_JOB_ID
-mkdir -p \${TMPDIR}
-
-INPUT_FILES=(\`cat ${output_dir}/slurm.compute_posterior_chunk.params\`)
-OUTPUT_FILES=("\${INPUT_FILES[@]/%/.out}")
-
-sort -k1,1 -k2,2n -S 32G -m \${OUTPUT_FILES[@]} > \${TMPDIR}/interval.all.bed
-
-bgzip -c \${TMPDIR}/interval.all.bed > ${output_dir}/interval.all.bed.gz
-tabix -0 -p bed ${output_dir}/interval.all.bed.gz
-
-rm -rf \${TMPDIR}
-__SCRIPT__
-
-JOB1=$(sbatch --export=ALL \
-	--job-name=${prefix}.compute_posterior_merge_sort \
-	--depend=afterok:${JOB0##* }\
-	${output_dir}/slurm.compute_posterior_merge_sort)
diff --git a/footprint_tools/__main__.py b/footprint_tools/__main__.py
new file mode 100644
index 0000000..3e6d95d
--- /dev/null
+++ b/footprint_tools/__main__.py
@@ -0,0 +1,45 @@
+import click
+
+# Set up console logging from package file
+# This configure root logger which is inherited 
+# by all modules/submodules
+import pkg_resources
+import logging, logging.config
+logging.config.fileConfig(pkg_resources.resource_filename(__name__, "logging.conf"))
+
+logger = logging.getLogger(__name__)
+
+import footprint_tools
+import footprint_tools.cli.learn_bm as learn_bm
+import footprint_tools.cli.learn_dm as learn_dm
+import footprint_tools.cli.detect as detect
+import footprint_tools.cli.learn_beta as learn_beta
+import footprint_tools.cli.plot_dm as plot_dm
+import footprint_tools.cli.post as post
+
+epilog = """See http://github.com/jvierstra/footprint-tools for extended documentation.
+
+Citation: Vierstra, J., Lazar, J., Sandstrom, R. et al. Global reference mapping of human transcription factor footprints. Nature 583, 729–736 (2020)
+
+Written by Jeff Vierstra (jvierstra@altius.org) (2015-2021). Software licensed under GNU General Public License version 3."""
+
+@click.group(epilog=epilog)
+@click.version_option(version=footprint_tools.__version__)
+def main():
+    """footprint_tools: analysis of digital genomic footprints 
+    
+    This software packages enables de novo detection and analysis of genomic footprints from DNase 
+    I data. The underying model simulates expected cleavage rates using a 6-mer DNase I cleavage
+    preference model combined with density smoothing. Statistical significance of per-nucleotide 
+    cleavages are computed from a series emperically fit negative binomial distribution. In addition
+    to footprint detection in a single, isolated dataset, this packages has a stastical framework
+    to jointly analyze 100s to 1000s of datasets in unison.
+    """
+    pass
+
+main.add_command(learn_bm.run)
+main.add_command(learn_dm.run)
+main.add_command(detect.run)
+main.add_command(learn_beta.run)
+main.add_command(plot_dm.run)
+main.add_command(post.run)
\ No newline at end of file
diff --git a/footprint_tools/cli/__init__.py b/footprint_tools/cli/__init__.py
new file mode 100644
index 0000000..5a4829a
--- /dev/null
+++ b/footprint_tools/cli/__init__.py
@@ -0,0 +1,2 @@
+"""
+"""
\ No newline at end of file
diff --git a/footprint_tools/cli/autocorr.py b/footprint_tools/cli/autocorr.py
new file mode 100644
index 0000000..5f2d4fc
--- /dev/null
+++ b/footprint_tools/cli/autocorr.py
@@ -0,0 +1,17 @@
+import click
+from click_option_group import optgroup
+
+from footprint_tools.cli.utils import tuple_args
+
+import logging
+logger = logging.getLogger(__name__)
+
+@click.command('autocorr')
+@click.argument('tabix_file')
+@click.option('--range', type=str, 
+    default="-25,25" callback=tuple_args(int),
+    help='')
+def run(tabix_file, range):
+    """Calculate nucleotide-level signal autocorrelation for genomic intervals
+    """
+    raise NotImplementedError
diff --git a/footprint_tools/cli/detect.py b/footprint_tools/cli/detect.py
new file mode 100644
index 0000000..8fef838
--- /dev/null
+++ b/footprint_tools/cli/detect.py
@@ -0,0 +1,298 @@
+import click
+from click_option_group import optgroup
+
+from multiprocessing import cpu_count
+
+import numpy as np
+import scipy as sp
+import pandas as pd
+import pysam
+
+from genome_tools import genomic_interval
+from footprint_tools import cutcounts
+from footprint_tools.modeling import bias, predict, dispersion
+from footprint_tools.stats import fdr, windowing
+from footprint_tools.cli.utils import (list_args, tuple_args, get_kwargs, 
+                                        verify_bam_file, verify_fasta_file, write_output_header,
+                                        write_stats_to_output, write_segments_to_output)
+from footprint_tools.data.process import process
+
+from tqdm import tqdm
+from tqdm.contrib.logging import logging_redirect_tqdm
+
+import logging
+logger = logging.getLogger(__name__)
+
+# kill numpy warnings
+np.seterr(all="ignore")
+
+class deviation_stats(process):
+    """Class that computes per-nucleotide cleavage
+    deviation statistics
+    """
+    def __init__(self, interval_file, bam_file, fasta_file, bm, dm, **kwargs):
+
+        self.intervals = pd.read_table(interval_file, header=None)
+
+        logger.info(f"BED file contains {len(self.intervals):,} regions")
+
+        self.bam_file = bam_file
+        self.fasta_file = fasta_file
+        self.bm = bm
+        self.dm = dm
+
+        # kwargs for cutcounts.bamfile
+        self.counts_reader_kwargs = get_kwargs(
+            [
+                'min_qual',
+                'remove_dups',
+                'remove_qcfail',
+                'offset'
+            ], 
+            kwargs)
+
+        self.fasta_reader_kwargs = {}
+
+        self.counts_predictor_kwargs = get_kwargs(
+            [
+                'half_win_width',
+                'smoothing_half_win_width',
+                'smoothing_clip'
+            ],
+            kwargs)
+
+        self.fdr_shuffle_n = kwargs['fdr_shuffle_n']
+        self.seed = kwargs['seed']
+
+        self.counts_reader = None
+        self.fasta_reader = None
+        self.count_predictor = None
+
+        self.win_pval_fn = lambda z: windowing.stouffers_z(np.ascontiguousarray(z), 3)
+
+    def cleanup(self):
+        """Needs to be implemented"""
+        pass
+
+    def __len__(self):
+        return len(self.intervals)
+
+    def __getitem__(self, index):
+        """Process data for a single interval"""
+        
+        # Open file handlers on first call. This avoids problems when
+        # parallel processing data with non-thread safe code (i.e., pysam)
+        if not self.counts_reader:
+            self.counts_reader = cutcounts.bamfile(self.bam_file, **self.counts_reader_kwargs)
+            self.fasta_reader = pysam.FastaFile(self.fasta_file, **self.fasta_reader_kwargs)
+            self.count_predictor = predict.prediction(self.counts_reader, self.fasta_reader, 
+                                                        self.bm, **self.counts_predictor_kwargs)
+
+        chrom, start, end = (self.intervals.iat[index, 0], 
+                             self.intervals.iat[index, 1], 
+                             self.intervals.iat[index, 2])
+
+        interval = genomic_interval(chrom, start, end)
+
+        obs, exp, _ = self.count_predictor.compute(interval)
+        obs = obs['+'][1:] + obs['-'][:-1]
+        exp = exp['+'][1:] + exp['-'][:-1]
+
+        assert len(obs) == len(exp)
+        n = len(obs)
+
+        if self.dm:
+            try:
+                pvals = self.dm.p_values(exp, obs)
+                win_pvals = self.win_pval_fn(pvals)
+
+                _, pvals_null = self.dm.sample(exp, self.fdr_shuffle_n)
+                win_pvals_null = np.apply_along_axis(self.win_pval_fn, 0, pvals_null)
+
+                efdr = fdr.emperical_fdr(win_pvals_null, win_pvals)
+            except Exception as e:
+                logger.warning(f"Error computing stats for '{interval.chrom}:{interval.start}-{interval.end}'")
+                pvals = win_pvals = efdr = np.ones(n) # should change to return 'nan'
+            finally:
+                stats = np.column_stack((exp, obs, -np.log(pvals), -np.log(win_pvals), efdr))
+        else:
+            stats = np.column_stack((exp, obs))
+
+        return {
+            'interval': interval,
+            'stats': stats
+        }
+
+@click.command(name='detect')
+@click.argument('interval_file')
+@click.argument('bam_file')
+@click.argument('fasta_file')
+@optgroup.group('Bias modeling options')
+@optgroup.option('--bias_model_file', type=click.STRING,
+    help='Use a k-mer model for sequence bias (supplied by file). '
+        'If argument is not provided the model defaults to uniform '
+        'sequence bias.')
+@optgroup.option('--half_win_width', type=click.INT,
+    default=5, show_default=True,
+    help='Half window width to apply bias model')
+@optgroup.group('Smoothing options')
+@optgroup.option('--smooth_half_win_width', type=click.INT,
+    default=50, show_default=True,
+    help='Half window width to apply smoothing model. When set to '
+        '0, no smoothing is applied.')
+@optgroup.option('--smooth_clip', type=click.FLOAT,
+    default=0.01, show_default=True,
+    help='Fraction of bases to clip when computing trimmed mean in the smoothing window')
+@optgroup.group("Statistics options")
+@optgroup.option('--dispersion_model_file', type=click.STRING,
+    help='Dispersion model for negative binomial tests. If argument '
+        'is not provided then no stastical output is provided. File is in '
+        'JSON format and generated using the command learn_dm')
+@optgroup.option('--fdr_shuffle_n', type=click.INT,
+    default=100, show_default=True,
+    help='Number of times to shuffle data for FDR calculation')
+@optgroup.group('Read filtering options')
+@optgroup.option('--min_qual', type=click.INT, 
+    default=1, show_default=True,
+    help='Ignore reads with mapping quality lower than this threshold')
+@optgroup.option('--keep_dups', type=click.BOOL,
+    default=False, show_default=True,
+    help='Keep duplicate reads')
+@optgroup.option('--keep_qcfail', type=click.BOOL,
+    default=False, show_default=True,
+    help='Keep QC-failed reads')
+@optgroup.group('Output options')
+@optgroup.option('--outprefix', type=click.STRING,
+    default='out', show_default=True,
+    help='Output prefix')
+@optgroup.option('--write_footprints', type=click.STRING, 
+    default="0.001,0.01,0.05", show_default=True, callback=list_args(float),
+    help='Output footprints at specified FDRs')
+@optgroup.group('Other options')
+@optgroup.option('--bam_offset', type=click.STRING,
+    default="0,-1", show_default=True, callback=tuple_args(int),
+    help='BAM file offset (enables support for other datatypes -- e.g. Tn5/ATAC)')
+@optgroup.option('--seed', type=click.INT,
+    help='Seed for random number generation (not currently used)')
+@optgroup.option('--n_threads', type=click.IntRange(1, cpu_count()),
+    default=cpu_count(), show_default=True,
+    help='Number of processors to use')
+@optgroup.option('--batch_size', type=click.INT,
+    default=100, show_default=True,
+    help='Batch size of intervals to process')
+def run(interval_file,
+        bam_file,
+        fasta_file,
+        bias_model_file=None,
+        dispersion_model_file=None,
+        min_qual=1,
+        keep_dups=False,
+        keep_qcfail=False,
+        bam_offset=(0, -1),
+        half_win_width=5,
+        smooth_half_win_width=50,
+        smooth_clip=0.01,
+        fdr_shuffle_n=50,
+        seed=None,
+        n_threads=cpu_count(),
+        batch_size=100,
+        write_footprints=[0.001, 0.01, 0.05],
+        outprefix='out'):
+    """Compute per-nucleotide cleavage deviation statistics	
+
+    \b
+    Inputs:
+    INTERVAL_FILE  Path to BED-formatted file contained genomic regions to be analyzed
+    BAM_FILE       Path to BAM-format tag alignment file
+    FASTA_FILE     Path to genome FASTA file (requires associated FASTA index in same folder
+                    (see documentation on how to create an index)
+    
+    \b
+    Outputs:
+    {outprefix}.bedgraph    bedGraph file written to 'outprefix'.bedgraph
+                            Columns: contig start start+1 obs exp -log(pval) -log(winpval) fdr
+    {outprefix}.fdr{t}.bed  BED3-format file with FDR thresholded footprints.
+                            FDR levels specified by '--write_footprints'.
+    """
+    proc_kwargs = {
+        "min_qual": min_qual,
+        "remove_dups": ~keep_dups,
+        "remove_qcfail": ~keep_qcfail,
+        "offset": bam_offset,
+        "half_win_width": half_win_width,
+        "smoothing_half_win_width": smooth_half_win_width,
+        "smoothing_clip": smooth_clip,
+        "fdr_shuffle_n": fdr_shuffle_n,
+        "seed": seed, # not used...yet
+    }
+
+    # Validate and load inputs
+    logger.info("Validating input files")
+    try:
+        verify_bam_file(bam_file)
+        verify_fasta_file(fasta_file)
+
+        # Load bias model (if specified), otherwise use the uniform model
+        if bias_model_file:
+            logger.info(f"Loading bias model from file {bias_model_file}")
+            bm = bias.kmer_model(bias_model_file)
+        else:
+            logger.info(f"No bias model file specified -- using uniform model")
+            bm = bias.uniform_model()
+
+        # Load dispersion model (if specified)
+        if dispersion_model_file:
+            logger.info(f"Loading dispersion model from file {dispersion_model_file}")
+            dm = dispersion.load_dispersion_model(dispersion_model_file)
+        else:
+            logger.info(f"No dispersion model file specified -- not be reporting base-level statistics and footprints")
+            dm = None
+            write_footprints = []
+
+    except IOError as e:
+        logger.critical(e)
+        raise click.Abort()
+
+    # Open output files
+    try:
+        # Output stats file
+        output_bedgraph_file = outprefix + '.bedgraph'
+        logger.info(f"Writing per-nucleotide stats to {output_bedgraph_file}")
+        output_bedgraph_filehandle = open(output_bedgraph_file , 'w')
+        write_output_header(["exp", "obs", "-log(pval)", "-log(winpval)", "fdr"], file=output_bedgraph_filehandle)
+
+        # Output footprints filex
+        output_bed_file_template = outprefix + '.fdr{0}.bed'
+        output_bed_filehandles = {}
+
+        if len(write_footprints) > 0:
+            logger.info(f"Writing FDR thresholded footprints to {output_bed_file_template.format('{threshold}')} for threshold \u22f2 {write_footprints}")
+            for t in write_footprints:
+                fh =  open(output_bed_file_template.format(t), 'w')
+                write_output_header(["name", "fdr"], file=fh, extra=f"thresholded @ FDR {t}")
+                output_bed_filehandles.update({t:fh})
+
+    except IOError as e:
+        logger.critical(e)
+        raise click.Abort()
+    
+    # Create data processor and iterator
+    dp = deviation_stats(interval_file, bam_file, fasta_file, bm, dm, **proc_kwargs)
+    dp_iter = dp.batch_iter(batch_size=batch_size, num_workers=n_threads)
+
+    with logging_redirect_tqdm():
+        
+        for batch in tqdm(dp_iter, colour='#cc951d'):
+
+            for interval, stats in zip(batch["interval"], batch["stats"]):
+                # write stats
+                write_stats_to_output(interval, stats, output_bedgraph_filehandle)
+
+                # write footprints
+                for thresh, fh in output_bed_filehandles.items():
+                    # fdr is last column in stats array
+
+                    write_segments_to_output(interval, stats[:,-1], thresh, file=fh, decreasing=True)
+
+    output_bedgraph_filehandle.close()
+    [f.close() for f in output_bed_filehandles.values()]
\ No newline at end of file
diff --git a/footprint_tools/cli/diff.py b/footprint_tools/cli/diff.py
new file mode 100644
index 0000000..9c11fc1
--- /dev/null
+++ b/footprint_tools/cli/diff.py
@@ -0,0 +1,10 @@
+import click
+from click_option_group import optgroup
+
+logger = logging.getLogger(__name__)    
+
+@click.command('diff_test')
+def run():
+    """Base-pair resolution test for differential DNase I cleavage
+    """
+    raise NotImplementedError
diff --git a/footprint_tools/cli/learn_beta.py b/footprint_tools/cli/learn_beta.py
new file mode 100644
index 0000000..5d6f742
--- /dev/null
+++ b/footprint_tools/cli/learn_beta.py
@@ -0,0 +1,90 @@
+import click
+from click_option_group import optgroup
+
+from yaspin import yaspin
+from yaspin.spinners import Spinners
+
+import numpy as np
+import scipy.stats
+
+import logging
+logger = logging.getLogger(__name__)
+
+@click.command(name='learn_beta')
+@click.argument('bedgraph_file')
+@click.option('--fdr_cutoff',
+    type=click.FLOAT, default=0.05,
+    help='Only consider nucleotides with FDR <= this value')
+@click.option('--exp_cutoff',
+    type=click.INT, default=10,
+    help='Only consider nucleotides with expected cleavages >= this value')
+@optgroup.group('Output options')
+@optgroup.option('--outfile',
+    type=click.STRING, default='beta.txt',
+    help='Output file path')
+def run(bedgraph_file, 
+        fdr_cutoff=0.05, 
+        exp_cutoff=10,
+        outfile='beta.txt'):
+    """Learn the parameters of a Beta distribution
+
+    \b
+    Input:
+    BEDGRAPH_FILE   Path to bedgraph file containg the per-nucleotide statistics
+                    generated from the 'detect' command.
+    
+    Outputs the fit Beta-distribution parameters to 'outfile'.
+
+    Note: This step is required to compute posterior footprint probabilities.
+    """
+
+    total_lines = total_passed = 0
+    obs_over_exp = []
+
+    logger.info(f"Filtering nucleotides with expected count >= {exp_cutoff} and within a FDR {fdr_cutoff} footprint")
+
+    with yaspin(Spinners.bouncingBar, text='Reading nucleotides') as sp:
+        
+        try:
+            filehandle = open(bedgraph_file, 'r')
+
+            for line in filehandle:
+                total_lines +=1
+                if total_lines % 500000 == 0:
+                    sp.text = "Reading nucleotides -- {:,} / {:,} passed filters".format(total_passed, total_lines)
+
+                fields = line.strip().split('\t')
+                exp = float(fields[3])
+                obs = float(fields[4])
+                fdr = float(fields[7])
+
+                if fdr <= fdr_cutoff and exp >= exp_cutoff:
+                    obs_over_exp.append( (obs+1)/(exp+1) )
+                    total_passed += 1
+        
+            filehandle.close()
+            
+        except IOError as e:
+            logger.critical(e)
+            click.Abort()
+        
+        sp.text = "Fitting Beta distribution"
+
+        obs_over_exp = np.array(obs_over_exp)
+
+        i = np.isfinite(obs_over_exp) & (obs_over_exp>0) & (obs_over_exp<1)
+        
+        prior = scipy.stats.beta.fit(obs_over_exp[i], floc = 0, fscale = 1)[0:2]
+
+    logger.info(f"{total_passed:,} positions used for fitting distribution")
+    logger.info(f"Beta distribution parameters: \u03b1 = {prior[0]}, \u03b2 = {prior[1]}")
+
+    logger.info(f"Writing parameters to file {outfile}")
+
+    try:
+        with open(outfile, 'w') as f:
+            print("%0.4f\t%0.4f" % (prior[0], prior[1]), file = f)
+    except IOError as e:
+        logger.critical(e)
+        click.Abort()
+
diff --git a/footprint_tools/cli/learn_bm.py b/footprint_tools/cli/learn_bm.py
new file mode 100644
index 0000000..3875c3c
--- /dev/null
+++ b/footprint_tools/cli/learn_bm.py
@@ -0,0 +1,34 @@
+import click
+
+from footprint_tools.cli.utils import (tuple_args, verify_bam_file, verify_fasta_file)
+
+from tqdm import tqdm
+from tqdm.contrib.logging import logging_redirect_tqdm
+
+import logging
+logger = logging.getLogger(__name__)
+
+@click.command(name='learn_bm')
+@click.argument('bam_file')
+@click.argument('fasta_file')
+@click.argument('mappability_file')
+@click.option('-k', type=click.INT,
+    default=6, show_default=True,
+    help='Width of k-mer model')
+@click.option('--mid', type=click.INT,
+    default=3, show_default=True,
+    help='Cleavage position of k-mer model relative to cut site')
+@click.option('--bam_offset', type=click.STRING,
+    default="0,-1", show_default=True, callback=tuple_args(int),
+    help='BAM file offset (enables support for other datatypes -- e.g. Tn5/ATAC)')
+def run(bam_file, fasta_file, mappability_file, k=6, mid=3, offset=(0, -1)):
+    """Learn a sequence bias model
+    """
+    # Validate and load inputs
+    logger.info("Validating input files")
+    try:
+        verify_bam_file(bam_file)
+        verify_fasta_file(fasta_file)
+    except IOError as e:
+        logger.critical(e)
+        raise click.Abort()
diff --git a/footprint_tools/cli/learn_dm.py b/footprint_tools/cli/learn_dm.py
new file mode 100644
index 0000000..2a02340
--- /dev/null
+++ b/footprint_tools/cli/learn_dm.py
@@ -0,0 +1,208 @@
+import click
+from click_option_group import optgroup
+
+from multiprocessing import cpu_count
+
+import numpy as np
+import scipy as sp
+import pandas as pd
+import pysam
+
+from genome_tools import genomic_interval
+from footprint_tools import cutcounts
+from footprint_tools.modeling import bias, predict, dispersion
+from footprint_tools.cli.utils import (tuple_args, get_kwargs, verify_bam_file, verify_fasta_file)
+from footprint_tools.data.process import process
+from footprint_tools.data.utils import numpy_collate_concat
+
+from tqdm import tqdm
+from tqdm.contrib.logging import logging_redirect_tqdm
+
+import logging
+logger = logging.getLogger(__name__)
+
+# kill numpy warnings
+np.seterr(all="ignore")
+
+class expected_counts(process):
+    """Computes observed and expected cleavage counts
+    """
+    def __init__(self, interval_file, bam_file, fasta_file, bm, **kwargs):
+       
+        self.intervals = pd.read_table(interval_file, header=None)
+
+        logger.info(f"BED file contains {len(self.intervals):,} regions")
+
+        self.bam_file = bam_file
+        self.fasta_file = fasta_file
+        self.bm = bm
+
+        # kwargs for cutcounts.bamfile
+        self.counts_reader_kwargs = get_kwargs(
+            [
+                'min_qual',
+                'remove_dups',
+                'remove_qcfail',
+                'offset'
+            ], 
+            kwargs)
+
+        self.fasta_reader_kwargs = {}
+
+        self.counts_predictor_kwargs = get_kwargs(
+            [
+                'half_win_width',
+            ],
+            kwargs)
+
+        self.counts_reader = None
+        self.fasta_reader = None
+        self.count_predictor = None
+
+    def __len__(self):
+        return len(self.intervals)
+
+    def __getitem__(self, index):
+        """Process data for a single interval"""
+        
+        # Open file handlers on first call. This avoids problems when
+        # parallel processing data with non-thread safe code (i.e., pysam)
+        if not self.counts_reader:
+            self.counts_reader = cutcounts.bamfile(self.bam_file, **self.counts_reader_kwargs)
+            self.fasta_reader = pysam.FastaFile(self.fasta_file, **self.fasta_reader_kwargs)
+            self.count_predictor = predict.prediction(self.counts_reader, self.fasta_reader, 
+                                                        self.bm, **self.counts_predictor_kwargs)
+
+        chrom, start, end = (self.intervals.iat[index, 0], 
+                             self.intervals.iat[index, 1], 
+                             self.intervals.iat[index, 2])
+
+        interval = genomic_interval(chrom, start, end)
+
+        obs, exp, _ = self.count_predictor.compute(interval)
+        obs = obs['+'][1:] + obs['-'][:-1]
+        exp = exp['+'][1:] + exp['-'][:-1]
+
+        return np.column_stack((exp, obs))
+
+@click.command(name='learn_dm')
+@click.argument('interval_file')
+@click.argument('bam_file')
+@click.argument('fasta_file')
+@optgroup.group('Bias modeling options')
+@optgroup.option('--bias_model_file', type=click.STRING,
+    help='Use a k-mer model for sequence bias (supplied by file). '
+        'If argument is not provided the model defaults to uniform '
+        'sequence bias.')
+@optgroup.option('--half_win_width', type=click.INT,
+    default=5, show_default=True,
+    help='Half window width to apply bias model')
+@optgroup.group('Read filtering options')
+@optgroup.option('--min_qual', type=click.INT, 
+    default=1, show_default=True,
+    help='Ignore reads with mapping quality lower than this threshold')
+@optgroup.option('--keep_dups', type=click.BOOL,
+    default=False, show_default=True,
+    help='Keep duplicate reads')
+@optgroup.option('--keep_qcfail', type=click.BOOL,
+    default=False, show_default=True,
+    help='Keep QC-failed reads')
+@optgroup.group('Output options')
+@optgroup.option('--outfile', type=click.STRING,
+    default='out', show_default=True,
+    help='Output prefix')
+@optgroup.group('Other options')
+@optgroup.option('--bam_offset', type=click.STRING,
+    default="0,-1", show_default=True, callback=tuple_args(int),
+    help='BAM file offset (enables support for other datatypes -- e.g. Tn5/ATAC)')
+@optgroup.option('--n_threads', type=click.IntRange(1, cpu_count()),
+    default=cpu_count(), show_default=True,
+    help='Number of processors to use')
+@optgroup.option('--batch_size', type=click.INT,
+    default=100, show_default=True,
+    help='Batch size of intervals to process')
+def run(interval_file,
+        bam_file,
+        fasta_file,
+        bias_model_file=None,
+        min_qual=1,
+        keep_dups=False,
+        keep_qcfail=False,
+        bam_offset=(0, -1),
+        half_win_width=5,
+        n_threads=cpu_count(),
+        batch_size=100,
+        outfile='dm.json'):
+    """Learn a negative binomial dispersion model 
+    
+    This creates expected cleavage counts from data corrected for intrinsic 
+    sequence preferences and then builds a dispersion model used for footprint
+    detection and analysis.
+
+    \b
+    Inputs:
+    INTERVAL_FILE  Path to BED-formatted file contained genomic regions to be analyzed
+    BAM_FILE       Path to BAM-format tag alignment file
+    FASTA_FILE     Path to genome FASTA file (requires associated FASTA index in same folder
+                    (see documentation on how to create an index)
+
+    Outputs a JSON-formated dispersion model
+    """
+
+    proc_kwargs = {
+        "min_qual": min_qual,
+        "remove_dups": ~keep_dups,
+        "remove_qcfail": ~keep_qcfail,
+        "offset": bam_offset,
+        "half_win_width": half_win_width,
+    }
+
+    # Validate and load input files
+    logger.info("Validating input files")
+    try:
+        verify_bam_file(bam_file)
+        verify_fasta_file(fasta_file)
+
+        # Load bias model (if specified), otherwise use the uniform model
+        if bias_model_file:
+            logger.info(f"Loading bias model from file {bias_model_file}")
+            bm = bias.kmer_model(bias_model_file)
+        else:
+            logger.info(f"No bias model file specified -- using uniform model")
+            bm = bias.uniform_model()
+
+    except IOError as e:
+        logger.critical(e)
+        raise click.Abort()
+
+    # Initiate data processor
+    dp = expected_counts(interval_file, bam_file, fasta_file, bm, **proc_kwargs)
+    dp_iter = dp.batch_iter(batch_size=batch_size, collate_fn=numpy_collate_concat, num_workers=n_threads)
+
+    hist_size = (200, 1000)
+    hist = np.zeros(hist_size, dtype=int)
+
+    with logging_redirect_tqdm():
+
+        for cnts in tqdm(dp_iter, colour='#cc951d'):
+            for i in range(cnts.shape[0]):
+                try:
+                    # expected, observed
+                    hist[int(cnts[i,0]),int(cnts[i, 1])] += 1
+                # ignore counts bigger than histogram bounds
+                except IndexError:
+                    pass
+
+    logger.info("Learning dispersion model")
+    
+    model = dispersion.learn_dispersion_model(hist)
+
+    logger.info("Writing dispersion model to {}".format(outfile))
+
+    try:
+        with open(outfile, "w") as output_filehandle:
+            print(dispersion.write_dispersion_model(model), file=output_filehandle)
+
+    except IOError as e:
+        logger.critical(e)
+        raise click.Abort()
\ No newline at end of file
diff --git a/footprint_tools/cli/plot_dm.py b/footprint_tools/cli/plot_dm.py
new file mode 100644
index 0000000..8186311
--- /dev/null
+++ b/footprint_tools/cli/plot_dm.py
@@ -0,0 +1,171 @@
+import click
+from click_option_group import optgroup
+
+import math
+import numpy as np
+import scipy.stats
+
+from matplotlib.pylab import rcParams
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib.ticker import MaxNLocator
+
+from footprint_tools.modeling import dispersion
+
+from footprint_tools.cli.utils import list_args
+
+import logging
+logger = logging.getLogger(__name__)
+
+def plot_model_mu(dm, ax=None, xlim=(0, 100)):
+    """
+    Plot model mu parameters
+    """
+    x = np.arange(xlim[0], xlim[1])
+
+    # Raw parameters
+    r = np.array([dm.r[i] for i in x])
+    p = np.array([dm.p[i] for i in x])
+    mu = p*r/(1.0-p)
+
+    # Smoothed parameters
+    fit_mu = np.array([dm.fit_mu(i) for i in x])
+
+    # Plot functions & appearance
+    ax.plot(x, mu, label='MLE neg. binomial fit')
+    ax.plot(x, fit_mu, label='Smoothed fit', ls='dashed')
+    ax.plot(xlim, xlim, label='y=x', color='grey', ls='dashed', zorder=-10)
+
+    ax.set_xlabel('Expected cleavage count')
+    ax.set_ylabel('Observed cleavages')
+
+    [ax.spines[loc].set_color('none') for loc in ['top', 'right']]
+
+    ax.xaxis.set_ticks_position('bottom')
+    ax.xaxis.set_tick_params(direction='out')
+    ax.xaxis.set(major_locator = MaxNLocator(4))
+
+    ax.yaxis.set_ticks_position('left')
+    ax.yaxis.set_tick_params(direction = 'out')
+    ax.yaxis.set(major_locator = MaxNLocator(4))
+
+    ax.legend()
+
+def plot_model_r(dm, ax=None, xlim=(1, 100)):
+    """
+    Plot model dispersion parameters
+    """
+    x = np.arange(xlim[0], xlim[1])
+
+    # Raw parameters
+    r = np.array([dm.r[i] for i in x])
+
+    # Smoothed parameters
+    fit_r = np.array([dm.fit_r(i) for i in x])
+
+    ax.plot(x, 1/r, label='MLE neg. binomial fit')
+    ax.plot(x, 1/fit_r, label='Smoothed fit', ls='dashed')
+
+    ax.set_xlabel("Expected cleavage count")
+    ax.set_ylabel("1/r")
+
+    [ax.spines[loc].set_color('none') for loc in ['top', 'right']]
+
+    ax.xaxis.set_ticks_position('bottom')
+    ax.xaxis.set_tick_params(direction='out')
+    ax.xaxis.set(major_locator = MaxNLocator(4))
+
+    ax.yaxis.set_ticks_position('left')
+    ax.yaxis.set_tick_params(direction = 'out')
+    ax.yaxis.set(major_locator = MaxNLocator(4))
+
+    ax.legend()
+
+def plot_histogram(dm, n=25, show_poisson=True, ax=None, xlim=(0, 125)):
+    """
+    Plots a density histogram of the observed cleavage counts
+    at an expected cleavage rate (n).
+    """
+    x = np.arange(xlim[0], xlim[1])
+
+    mu = dm.fit_mu(n)
+    r = dm.fit_r(n)
+
+    # Raw observed counts
+    ax.bar(x, dm.h[n,x[0]:(x[-1]+1)]/np.sum(dm.h[n,:]), width=1, color='lightgrey', label="Observed")
+
+    # NB fit
+    y_nbinom=scipy.stats.nbinom.pmf(x, r, r/(r+mu))
+    ax.plot(x, y_nbinom, color="red", label="Negative binomial")
+
+    # Poisson
+    if show_poisson:
+        y_pois=scipy.stats.poisson.pmf(x, mu=n)
+        ax.plot(x, y_pois, color="blue", label="Poisson")
+
+    ax.set_xlim(x[0], x[-1])
+
+    [ax.spines[loc].set_visible(False) for loc in ["top", "right"]]
+    ax.set_xlabel("Observed cleavage count")
+    ax.set_ylabel("Density")
+
+    ax.set_title(f"{n} expected cleavages")
+
+    ax.legend()
+
+@click.command(name='plot_dm')
+@click.argument('dispersion_model_file')
+@optgroup.group('Output options')
+@optgroup.option('--histograms',
+    type=click.STRING, default="5,25,50,75", callback=list_args(int),
+    help='Plot histograms of observed counts at site with expected counts (comma-seperated list)')
+@optgroup.option('--outfile',
+    type=click.STRING, default='dm.pdf',
+    help='Output file path for plot (suffix determines image format)')
+def run(dispersion_model_file, histograms=[15,25,50,75], outfile='dm.pdf'):
+    """Diagnostic plotting of a dispersion model
+    
+    Output:
+        dm.pdf - a PDF file with plots
+    """
+    
+    dm = dispersion.load_dispersion_model(dispersion_model_file)
+
+    plt_params = {
+        'legend.fontsize': 'xx-small',
+        'axes.labelsize': 'small',
+        'axes.titlesize': 'small',
+        'xtick.labelsize': 'x-small',
+        'ytick.labelsize': 'x-small'}
+    rcParams.update(plt_params)
+
+    npanels = len(histograms)+2
+    ncols = 2
+    nrows = math.ceil(npanels/ncols)
+
+    fig = plt.figure()
+    gs = gridspec.GridSpec(nrows, ncols,  wspace=0.75, hspace=0.75)
+
+    logger.info("Plotting model parameters")
+
+    ax = fig.add_subplot(gs[0,0])
+    plot_model_mu(dm, ax)
+
+    ax = fig.add_subplot(gs[0,1])
+    plot_model_r(dm, ax)
+
+    logger.info(f"Plotting histograms - {histograms}")
+
+    for i, n in enumerate(histograms):
+        row_index = (i // ncols) + 1
+        col_index = i % ncols
+        ax = fig.add_subplot(gs[row_index, col_index])
+        plot_histogram(dm, n=n, ax=ax)
+
+    fig.set_size_inches(2.5*ncols, 2*nrows)
+    
+    logger.info(f"Saving plots to {outfile}")	
+
+    plt.savefig(outfile, transparent=True)
+
+    return 0
diff --git a/footprint_tools/cli/post.py b/footprint_tools/cli/post.py
new file mode 100644
index 0000000..aeb368e
--- /dev/null
+++ b/footprint_tools/cli/post.py
@@ -0,0 +1,210 @@
+import click
+from click_option_group import optgroup
+
+from multiprocessing import cpu_count
+
+import numpy as np
+import scipy as sp
+import pandas as pd
+import pysam
+
+from genome_tools import genomic_interval
+import footprint_tools
+from footprint_tools.modeling import dispersion
+from footprint_tools.stats import posterior
+from footprint_tools.cli.utils import (verify_tabix_file, write_stats_to_output, 
+                                        write_output_header)
+from footprint_tools.data.process import process
+from footprint_tools.data.utils import numpy_collate_concat
+
+from tqdm import tqdm
+from tqdm.contrib.logging import logging_redirect_tqdm
+
+import logging
+logger = logging.getLogger(__name__)
+
+# kill numpy warnings
+np.seterr(all="ignore")
+
+# columns required in sample data file
+required_sample_data_cols = ["id", "tabix_file", "dm_file", "beta_a", "beta_b"]
+
+class posterior_stats(process):
+    def __init__(self, interval_file, samples_data, fdr_cutoff):
+
+        self.intervals = pd.read_table(interval_file)
+        self.samples_data = samples_data
+        self.fdr_cutoff = fdr_cutoff
+
+        self.tabix_files = None # these get loaded on first call of __getitem__
+        self.disp_models = [dispersion.load_dispersion_model(fn) for fn in self.samples_data["dm_file"]]
+        self.betas = self.samples_data[["beta_a", "beta_b"]].to_numpy()
+
+    def _open_tabix_files(self):
+        self.tabix_files = [pysam.TabixFile(fn) for fn in self.samples_data["tabix_file"]]
+
+    def _load_data(self, interval):
+        assert self.tabix_files or len(self.tabix_files) > 0
+        
+        n = len(self.tabix_files)
+        m = len(interval)
+
+        obs = exp = fdr = w = np.zeros((n, m), dtype=np.float)
+        i = j = 0
+
+        for i, tbf in enumerate(self.tabix_files):
+            try:
+                for row in tbf.fetch(interval.chrom, interval.start,interval.end, parser=pysam.asTuple()):
+                    j = int(row[1])-interval.start
+                    obs[i, j] = np.float(row[3])
+                    exp[i, j] = np.float(row[4])
+                    fdr[i, j] = np.float(row[7])
+                    w[i, j] = 1
+            except:
+                pass
+        
+        return obs, exp, fdr, w
+
+    def cleanup(self):
+        """Clean-up filehandlers"""
+        if self.tabix_files and len(self.tabix_files) > 0:
+            [tbf.close() for tbf in self.tabix_files]
+        self.tabix_files = None
+
+    def __len__(self):
+        return len(self.intervals)
+    
+    def __getitem__(self, index):
+        """Process data for a single interval"""
+
+        if len(self.tabix_files)==0:
+            self._open_tabix_files()
+
+        chrom, start, end = (self.intervals.iat[index, 0], 
+                        self.intervals.iat[index, 1], 
+                        self.intervals.iat[index, 2])
+
+        interval = genomic_interval(chrom, start, end)
+        
+        obs, exp, fdr, w = self._load_data(interval)
+
+        prior = posterior.compute_prior_weighted(fdr, w, cutoff=self.fdr_cutoff)   
+        delta = posterior.compute_delta_prior(obs, exp, fdr, self.beta, cutoff=self.fdr_cutoff)
+        ll_on = posterior.log_likelihood(obs, exp, self.dms, delta=delta, w=3) 
+        ll_off = posterior.log_likelihood(obs, exp, self.disp_models, w=3)
+
+        post = -posterior.posterior(prior, ll_on, ll_off)
+        post[post<=0] = 0.0
+
+        return {
+            'interval': interval,
+            'post': post,
+        }
+
+@click.command(name='posterior')
+@click.argument('sample_data_file')
+@click.argument('interval_file')
+@click.option('--fdr_cutoff', type=click.FLOAT,
+    default=0.05, show_default=True,
+    help='FDR cutoff to use when computing priors')
+@optgroup.group('Output options')
+@optgroup.option('--post_cutoff', type=click.FLOAT,
+    default=0.2, show_default=True,
+    help='Print only positions where that maximum posterior across '
+        'all samples meets this threshold. Used to control for '
+        'file size.')
+@optgroup.option('--outprefix', type=click.STRING,
+    default='out', show_default=True,
+    help='Output prefix')
+@optgroup.group('Other options')
+@optgroup.option('--n_threads', type=click.IntRange(1, cpu_count()),
+    default=cpu_count(), show_default=True,
+    help='Number of processors to use')
+@optgroup.option('--batch_size', type=click.INT,
+    default=100, show_default=True,
+    help='Batch size of intervals to process')
+def run(sample_data_file,
+        interval_file, 
+        fdr_cutoff=0.05,
+        post_cutoff=0.2,
+        batch_size=100,
+        n_threads=cpu_count(),
+        outprefix='out'):
+    """Compute footprint posterior probabilities
+
+    Applies an emperical Bayesian approach to compute the posterior probability a
+    nucleotide is protected by jointly analyzing many datasets.
+
+    \b
+    Inputs:
+    INTERVAL_FILE       Path to BED-formatted file contained genomic regions to be analyzed
+    SAMPLE_DATA_FILE    Path to file that contains sample data. File is tab-delimited with the
+                        following columns:
+
+    \b
+                        id          Sample identifier (unique)
+                        tabix_file  Path to TABIX-format cleavage statistics file
+                        dm_file     Path to dataset JSON-encoded dispersion model file
+                        beta_a      \u03b1 parameter (see 'learn_beta' command)
+                        beta_b      \u03b2 parameter
+    \b
+                        Note: File must contain a header row. Lines ignored when '#' is first
+                        character.
+
+    \b
+    Output:
+    {outprefix}.bedgraph    bedGraph file written to 'outprefix'.bedgraph .Columns are
+                                'contig start start+1 -log(1-p)_1 ... -log(1-p)_n'
+                            where n is the total number samples. Note that values are
+                            1-posterior probability. Samples (columns) are in the same 
+                            order as the sample data file.
+    """
+
+    logger.info(f"Loading sample data file {sample_data_file}")
+    try:
+        sample_data = pd.read_table(sample_data_file, header=0, comment='#')
+        if not all([col in sample_data.columns for col in required_sample_data_cols]):
+            raise ValueError(f"Sample data file does not contain the required columns: {required_sample_data_cols}")
+    except Exception as e:
+        logger.critical(e)
+        raise click.Abort()
+
+    # Validate and load inputs
+    logger.info("Validating input files")
+    try:
+        [verify_tabix_file(fn) for fn in sample_data["tabix_file"]]
+        [open(fn).close() for fn in sample_data["dm_file"]]
+
+    except IOError as e:
+        logger.critical(e)
+        raise click.Abort()
+    
+    # Open output files
+    try:
+        # Output stats file
+        output_bedgraph_file = outprefix + '.bedgraph'
+        logger.info(f"Writing per-nucleotide stats to {output_bedgraph_file}")
+        output_bedgraph_filehandle = open(output_bedgraph_file , 'w')
+
+        #write header lines
+        write_output_header(sample_data["id"], output_bedgraph_filehandle)
+
+    except IOError as e:
+        logger.critical(e)
+        raise click.Abort()
+
+    # Create data processor and iterator
+    dp = posterior_stats(interval_file, sample_data, fdr_cutoff)
+    dp_iter = dp.batch_iter(batch_size=batch_size, num_workers=n_threads)
+
+    # filter function to apply when writing posteriors to file
+    filter_fn = lambda x: np.min(x, axis=1) > post_cutoff
+
+    with logging_redirect_tqdm():
+
+        for batch in tqdm(dp_iter, colour='#cc951d'):
+            
+            for interval, stats in zip(batch["interval"], batch["stats"]):
+                write_stats_to_output(interval, stats, output_bedgraph_filehandle, filter_fn=filter_fn)
+    
+    output_bedgraph_filehandle.close()
diff --git a/footprint_tools/cli/utils.py b/footprint_tools/cli/utils.py
new file mode 100644
index 0000000..d8a3f02
--- /dev/null
+++ b/footprint_tools/cli/utils.py
@@ -0,0 +1,174 @@
+"""
+Argument validation functions
+"""
+import click
+import click_option_group
+
+def tuple_args(value_type=int):
+    def _parse_tuple(ctx, params, value):
+        """Function to parse a tuple from command line"""
+        try:
+            items = tuple(map(value_type, value.split(',')))
+            assert len(items) == 2
+        except:
+            raise click.BadOptionUsage(f'needs to be a comma-delimited tuple of type {value_type.__name__}')
+        return items
+    return _parse_tuple
+
+def list_args(value_type=int):
+    def _parse_list(ctx, params, value):
+        """Function to parse a list from command line"""
+        try:
+            items = list(map(value_type, value.split(',')))
+        except:
+            raise click.BadOptionUsage(f'needs to be a comma-delimited list of type {value_type.__name__}')
+        return items
+    return _parse_list
+
+def get_kwargs(keys, kwargs):
+    return {k:kwargs[k] for k in keys if k in kwargs}
+
+"""
+File validation functions
+"""
+import pysam
+
+def verify_bam_file(fn):
+    """Tries to open a file, raises IOError with problems"""
+    try:
+        pysam.AlignmentFile(fn).close()
+    except IOError:
+        raise IOError(f"No such file: {fn}")
+    except ValueError:
+        raise IOError(f"BAM-index not found for {fn}")
+
+def verify_tabix_file(fn):
+    """Tries to open a file, raises IOError with problems"""
+    try:
+        pysam.TabixFile(fn).close()
+    except IOError:
+        raise IOError(f"No such file: {fn}")
+    except ValueError:
+        raise IOError(f"TABIX-index not found for {fn}")
+
+def verify_fasta_file(fn):
+    """Tries to open a file, raises IOError with problems"""
+    try:
+        pysam.FastaFile(fn).close()
+    except IOError:
+        raise IOError(f"No such file: {fn}")
+    except ValueError:
+        raise IOError(f"FASTA-index not found for {fn}")
+
+"""
+Writer functions
+"""
+import sys
+import numpy as np
+
+import footprint_tools
+from footprint_tools.stats import utils
+
+def write_output_header(columns, file=sys.stdout, delim='\t', include_name=True, extra=None):
+    """Write header to ouptput file
+    
+    Parameters
+    ----------
+    columns : list
+        Columns in file (apart from standard chrom, start, end)
+    file : filehandle
+        Filehandle to write header to
+    delim : str
+        Delimiter to use when printing output
+    include_name : bool
+        Include name column before writting additional columns. Makes
+        the file consistent with BED3+2 format
+    """
+    file.write(f"# generated by {footprint_tools.__name__} version {footprint_tools.__version__}\n")
+    if extra:
+        file.write('# '+extra+'\n')
+    file.write(f"# chrom{delim}start{delim}end")
+    if include_name:
+        file.write(f"{delim}name")
+    file.write(delim + delim.join(columns)+'\n')
+
+
+
+def write_stats_to_output(interval, stats, file=sys.stdout, delim='\t', filter_fn=None, fmt_string='0.4f'):
+    """Write per-nucleotide statistics to file
+    
+    The input stats array can be of any width the function just writes all
+    columns
+
+    Parameters
+    ----------
+    interval : genomic_interval
+        Genomic regions corresponding to stats
+    stats : ndarray
+        2-D array
+    file : filehandle
+        File handle to write output
+    delim : str
+        Delimiter to use when printing output
+    filter_fn: callable
+        Function to apply that determines which positions get
+        written to file
+    fmt_string: str
+        Format string to apply when writing scores
+
+    Notes
+    -----
+    The filter_fn is a function that returns an boolean array of same length of stats
+    specifying whether a position should be written. For example:
+
+        lambda x: np.max(x, axis=1) >= 0.2
+
+    If filter_fn is not specified all rows are outputted
+    """    
+    chrom = interval.chrom
+    start = interval.start
+
+    idxs = np.nonzero(filter_fn(stats)) if filter_fn else range(stats.shape[0])
+
+    for i in idxs:
+        out = f'{chrom}{delim}{start+i}{delim}{start+i+1}{delim}'
+        out += delim.join([('{0:'+fmt_string+'}').format(val) for val in stats[i,:]])
+        file.write(out+'\n')
+
+def write_segments_to_output(interval, stats, threshold, name='.', file=sys.stdout,
+                            delim='\t', score_fn=np.min, decreasing=False, fmt_string='0.4f'):
+    """Write footprints to file
+    
+    Parameters
+    ----------
+    interval : genomic_interval
+        Genomic regions corresponding to stats
+    stats : ndarray
+        1D array
+    threshold : int or float
+        Threshold to use when performing segmentations
+    file : filehandle
+        File handle to write output
+    delim : str
+        Delimiter to use when printing output
+    score_fn: callable
+        Function to apply when computing score within the segments
+    decreasing: bool
+        Values below threshold are outputted
+    name: str
+        Name for BED entry
+    fmt_string: str
+        Format string to apply when writing scores
+    """
+    chrom = interval.chrom
+    start = interval.start
+
+    assert stats.ndim == 1
+
+    segments = utils.segment(stats, threshold, 3, decreasing=decreasing)
+
+    for s, e in segments:
+        score = score_fn(stats[s:e])        
+        out = f"{chrom}{delim}{start+s}{delim}{start+e}{delim}{name}{delim}"
+        out += ('{0:'+fmt_string+'}').format(score)
+        file.write(out+'\n')
diff --git a/footprint_tools/cutcounts.py b/footprint_tools/cutcounts.py
index ff2b18f..f852694 100644
--- a/footprint_tools/cutcounts.py
+++ b/footprint_tools/cutcounts.py
@@ -1,410 +1,467 @@
-# Copyright 2015 Jeff Vierstra
-
+"""
+This modules contains classes and functions to compute cleavage counts
+directly from an alignment file.
+"""
 import genome_tools
 import pysam
 
 import numpy as np
 from collections import defaultdict
 
+import logging
+logger = logging.getLogger(__name__)
+
 class ReadError(Exception):
-	
-	ERROR_ALIGNMENT = (0, "Read alignment problematic (QC fail, duplicate, or MAPQ < 1)")
-	ERROR_5PROXIMITY = (1, "Variant too close to 5' end of tag") 
-	ERROR_BASEQ = (2, "Base quality < 20")
-	ERROR_GENOTYPE = (3, "Base does not match reference or expected alternate allele")
-	ERROR_MISMATCH = (4, "Read contains too many mismatches")
-	
-	def __init__(self, e):
-		self.value = e[0]
-		self.message = e[1]
+    
+    ERROR_ALIGNMENT = (0, "Read alignment problematic (QC fail, duplicate, or MAPQ < 1)")
+    ERROR_5PROXIMITY = (1, "Variant too close to 5' end of tag") 
+    ERROR_BASEQ = (2, "Base quality < 20")
+    ERROR_GENOTYPE = (3, "Base does not match reference or expected alternate allele")
+    ERROR_MISMATCH = (4, "Read contains too many mismatches")
+    
+    def __init__(self, e):
+        self.value = e[0]
+        self.message = e[1]
 
 
 class GenotypeError(Exception):
-	pass
+    pass
 
 class bamfile(object):
-	"""Class to access BAM files
-	
-	Attributes
-	----------
-	min_qual : int
-	    Filter reads by minimim mapping quality (MAPQ)
-	offset : tuple
-	    Position offsets to apply to the `+` and `-` strands,. DNase I data (0, -1).
+    """Class to access BAM files
+    
+    Attributes
+    ----------
+    min_qual : int
+        Filter reads by minimim mapping quality (MAPQ)
+    offset : tuple
+        Position offsets to apply to the `+` and `-` strands,. DNase I data (0, -1).
         Tn5-derived data would use (4,-5). (default = (0, -1))
-	remove_dups : bool
-	    Remove reads with duplicate flag (512) set
-	remove_qcfail : bool
-	    Remove reads with QC fail flag (1024) set
-	samfile : pysam.Samfile
-	    SAM/BAM file object
-	"""
-	
-	def __init__(self, filepath, min_qual = 1, remove_dups = False, remove_qcfail = True, offset = (0, -1), return_type="counts"):
-		"""Constructor
-		
-		Parameters
-		----------
-		filepath : str
-		    Fileath to SAM/BAM file
-		min_qual : int, optional
-		    Filter reads by minimim mapping quality (MAPQ)
-		remove_dups : bool, optional
-		    Remove reads with duplicate flag (512) set
-		remove_qcfail : bool, optional
-		    Remove reads with QC fail flag (1024) set
-		offset : tuple, optional
-		    Position offsets to apply to the `+` and `-` strands (default =(0, -1))
-		return_type: string, optional
-            Type of data to return; currently supports "counts" (stranded cutcounts) or "fragments" (mapped fragments).
-		Raises
-		------
-		IOError
-		    If file is not found, has no index, or something else bad happened while trying to open BAM file
-		"""
-
-		try:
-			self.samfile = pysam.Samfile(filepath, "rb")
-		except:
-			raise IOError("Cannot open BAM file: %s" % filepath)
-
-		self.offset = offset #-1 # a hack for the mis-aligned data from 2010
-		self.min_qual = min_qual
-		self.remove_dups = remove_dups
-		self.remove_qcfail = remove_qcfail
+    remove_dups : bool
+        Remove reads with duplicate flag (512) set
+    remove_qcfail : bool
+        Remove reads with QC fail flag (1024) set
+    samfile : pysam.Samfile
+        SAM/BAM file object
+    """
+    
+    def __init__(self, 
+                filepath, 
+                min_qual = 1, 
+                remove_dups = False,
+                remove_qcfail = True,
+                offset = (0, -1)):
+        """Constructor
+        
+        Parameters
+        ----------
+        filepath : str
+            Fileath to SAM/BAM file
+        min_qual : int, optional
+            Filter reads by minimim mapping quality (MAPQ)
+        remove_dups : bool, optional
+            Remove reads with duplicate flag (512) set
+        remove_qcfail : bool, optional
+            Remove reads with QC fail flag (1024) set
+        offset : tuple, optional
+            Position offsets to apply to the `+` and `-` strands (default =(0, -1))
+   
+        Raises
+        ------
+        IOError
+            If file is not found
+        ValueError
+            If file has not associated index
+        """
+
+        try:
+            self.samfile = pysam.Samfile(filepath, "rb")
+        except:
+            raise IOError("Cannot open BAM file: %s" % filepath)
+
+        self.offset = offset #-1 # a hack for the mis-aligned data from 2010
+        self.min_qual = min_qual
+        self.remove_dups = remove_dups
+        self.remove_qcfail = remove_qcfail
+        
+    def close(self):
+        """Closes BAM file
+        """
+        if self.samfile:
+            self.samfile.close()
+        return True
+
+    def __del__(self):
+        return self.close()
+
+    def validate_read(self, read):
+        """Validate BAM tag
+        
+        Parameters
+        ----------
+        read : pysam.AlignedSegment
+            Read from BAM/SAM file
+        
+        Returns
+        -------
+        read: pysam.AlignedSegment
+            Same read as input
         
-	def close(self):
-		if self.samfile:
-			self.samfile.close()
-		return True
-
-	def __del__(self):
-		return self.close()
-
-	def validate_read(self, read):
-		"""
-		Validate BAM tag
-		
-		Parameters
-		----------
-		read : pysam.AlignedSegment
-		    Read from BAM/SAM file
-		
-		Returns
-		-------
-		pysam.AlignedSegment
-		    Same read as input
-		
-		Raises
-		------
-		ReadError
-		    Raises error if read fails QC flag, is a duplicate or MAPQ < minimum
-		"""
-
-		if self.remove_qcfail and read.is_qcfail:
-			raise ReadError(ReadError.ERROR_ALIGNMENT)
-		if self.remove_dups and read.is_duplicate:
-			raise ReadError(ReadError.ERROR_ALIGNMENT)
-		if read.mapping_quality < self.min_qual:
-			raise ReadError(ReadError.ERROR_ALIGNMENT)
-
-		return read
-
-	def get_read_mate(self, read):
-		"""
-		Fetch the mate pair for paired-end reads
-		
-		Parameters
-		----------
-		read : pysam.AlignedSegment
-			One  of the read-pairs
-		
-		Returns
-		-------
-		pysam.AlignedSegment
-		    The corresponding mate read
-		
-		"""
-
-		fpos = self.samfile.tell()
-		try:
-			mate = self.samfile.mate(read)
-		except ValueError:
-			mate = None
-		finally:
-			self.samfile.seek(fpos)
-		return mate
-
-	def read_pair_generator(self, chrom, start, end):
-		"""
-		Generator function that returns sequencing tags within a given region
-		
-		Parameters
-		----------
-		chrom : str
-		    Chromosome
-		start : int
-		    Start coordinate
-		end : int
-		    End coordinate
-		
-		Yields
-		------
-		tuple
-		    A tuple of :class:`pysam.AlignedSegments`. Elements may be `NoneType` if single-end sequencing or one of pairs falls outisde of query range.
-		
-		"""
-
-		read_dict = defaultdict(lambda: [None, None])
-
-		for read in self.samfile.fetch(chrom, max(start-10, 0), end+10):
-
-			try:
-
-				self.validate_read(read)
-
-				# Single-end:,just pass the read
-				if not read.is_paired:
-					yield read, None 
-				
-				#  Pair-end; do some validation
-				else:
-
-					if not read.is_proper_pair or read.is_secondary or read.is_supplementary:
-						continue
-
-					qname = read.query_name
-
-					if qname not in read_dict:
-						if read.is_read1:
-							read_dict[qname][0] = read
-						else:
-							read_dict[qname][1] = read
-					else:
-						if read.is_read1:
-							yield read, read_dict[qname][1]
-						else:
-							yield read_dict[qname][0], read
-						del read_dict[qname]
-
-			except ReadError as e:
-				continue
-
-		"""Flush out the rest of dictionary. (for example if one the 
-		mates wasn't in the original region). It might be reasonable to 
-		manually grab the remaining read using `__get_read_mate`"""
-		for k, reads in read_dict.items():
-			yield reads[0], reads[1]
+        Raises
+        ------
+        ReadError
+            Raises error if read fails QC flag, is a duplicate or MAPQ < minimum
+        """
+
+        if self.remove_qcfail and read.is_qcfail:
+            raise ReadError(ReadError.ERROR_ALIGNMENT)
+        if self.remove_dups and read.is_duplicate:
+            raise ReadError(ReadError.ERROR_ALIGNMENT)
+        if read.mapping_quality < self.min_qual:
+            raise ReadError(ReadError.ERROR_ALIGNMENT)
+
+        return read
+
+    def get_read_mate(self, read):
+        """Fetch the mate pair for paired-end reads
+        
+        Parameters
+        ----------
+        read : pysam.AlignedSegment
+            One  of the read-pairs
+        
+        Returns
+        -------
+        mate: pysam.AlignedSegment
+            The corresponding mate read two input read
+        """
+
+        fpos = self.samfile.tell()
+        try:
+            mate = self.samfile.mate(read)
+        except ValueError:
+            mate = None
+        finally:
+            self.samfile.seek(fpos)
+        return mate
+
+    def read_pair_generator(self, chrom, start, end):
+        """Generator function that returns sequencing tags within a given region
+        
+        Parameters
+        ----------
+        chrom : str
+            Chromosome
+        start : int
+            Start coordinate
+        end : int
+            End coordinate
+        
+        Yields
+        ------
+        reads: tuple
+            A tuple of :class:`pysam.AlignedSegment`. Elements may be `NoneType` if 
+            single-end sequencing or one of pairs falls outisde of query range.
+        """
+
+        read_dict = defaultdict(lambda: [None, None])
+
+        for read in self.samfile.fetch(chrom, max(start-10, 0), end+10):
+
+            try:
+
+                self.validate_read(read)
+
+                # Single-end:,just pass the read
+                if not read.is_paired:
+                    yield read, None 
+                
+                #  Pair-end; do some validation
+                else:
+
+                    if not read.is_proper_pair or read.is_secondary or read.is_supplementary:
+                        continue
+
+                    qname = read.query_name
+
+                    if qname not in read_dict:
+                        if read.is_read1:
+                            read_dict[qname][0] = read
+                        else:
+                            read_dict[qname][1] = read
+                    else:
+                        if read.is_read1:
+                            yield read, read_dict[qname][1]
+                        else:
+                            yield read_dict[qname][0], read
+                        del read_dict[qname]
+
+            except ReadError as e:
+                continue
+
+        """Flush out the rest of dictionary. (for example if one the 
+        mates wasn't in the original region). It might be reasonable to 
+        manually grab the remaining read using `__get_read_mate`"""
+        for k, reads in read_dict.items():
+            yield reads[0], reads[1]
             
-	def __add_read(self, read, fw, rev):
-		"""
-
-		"""
-		if read.is_reverse:
-			a = int(read.reference_end)+self.offset[1]
-			rev[a] = rev.get(a, 0.0) + 1.0
-		else:
-			a = int(read.reference_start)+self.offset[0]
-			fw[a] = fw.get(a, 0.0) + 1.0
-
-	def __get_fragment(self, read):
-		"""
-		Returns a fragement from a mapped read
-		"""
-		tlen = read.template_length
-		if read.is_reverse:
-			end = int(read.reference_end)+self.offset[1]
-			start = end + tlen
-		else:
-			start = int(read.reference_start)+self.offset[0]
-			end = start + tlen
-
-		return genome_tools.genomic_interval(read.reference_name, start, end)   
-
-	def lookup(self, interval):
-		"""
-
-		"""	
-		chrom = interval.chrom
-		start = interval.start
-		end = interval.end
-		flip = True if interval.strand == '-' else False
-
-		tmp_fw = {}
-		tmp_rev = {}
-		reads = []
-
-		for read1, read2 in self.read_pair_generator(chrom, max(start-10, 0), end+10):
-			if read1:
-				self.__add_read(read1, tmp_fw, tmp_rev)
-				reads.append(self.__get_fragment(read1))
-			if read2:
-				self.__add_read(read2, tmp_fw, tmp_rev)
+    def __add_read(self, read, fw, rev):
+        """Add read to output dictionary
+
+        Parameters
+        ----------
+        read: pysam.AlignedSegment
+            Aligned read segment to add
+        fw: dict
+            Dictionary of read counts on forward strand
+        rev: dict
+            Dictionary of read counts on reverse strand
+        """
+        if read.is_reverse:
+            a = int(read.reference_end)+self.offset[1]
+            rev[a] = rev.get(a, 0.0) + 1.0
+        else:
+            a = int(read.reference_start)+self.offset[0]
+            fw[a] = fw.get(a, 0.0) + 1.0
+
+    def __get_fragment(self, read):
+        """Returns a fragment from a mapped read
+        
+        Parameters
+        ----------
+        read : pysam.AlignedSegment
+            An aligned read
+
+        Returns
+        -------
+        fragment: genomic_interval
+            Fragment that corresponds to `read`. Start and end coordinates
+            correspond to the 5' ends of the paired reads.
+        """
+        tlen = read.template_length
+        if read.is_reverse:
+            end = int(read.reference_end)+self.offset[1]
+            start = end + tlen
+        else:
+            start = int(read.reference_start)+self.offset[0]
+            end = start + tlen
+
+        return genome_tools.genomic_interval(read.reference_name, start, end)   
+
+    def lookup(self, interval):
+        """Lookup reads in a defined genomic region
+
+        Parameters
+        ----------
+        intervals: iterable (genomic_interval)
+
+        Returns
+        -------
+        counts: dict
+            Dictionary of read counts (keys: '+' or '-), which contain arrays
+            with counts on each strand
+        """
+
+        chrom = interval.chrom
+        start = interval.start
+        end = interval.end
+        flip = True if interval.strand == '-' else False
+
+        tmp_fw = {}
+        tmp_rev = {}
+        reads = []
+
+        for read1, read2 in self.read_pair_generator(chrom, max(start-10, 0), end+10):
+            if read1:
+                self.__add_read(read1, tmp_fw, tmp_rev)
+                reads.append(self.__get_fragment(read1))
+            if read2:
+                self.__add_read(read2, tmp_fw, tmp_rev)
             
-		fw_cutarray = np.array([tmp_fw.get(i, 0.0) for i in range(start, end)])
-		rev_cutarray = np.array([tmp_rev.get(i, 0.0) for i in range(start, end)])
-
-		return {
- 			"+": rev_cutarray[::-1] if flip else fw_cutarray, 
-			"-": fw_cutarray[::-1] if flip else rev_cutarray,
-			"fragments": reads
-			}
-
-	def __validate_genotype(self, read, pos, ref, alt):
-		"""
-
-		Returns
-		-------
-		bool
-			True if read contains reference allele and passes filters
-			False if read contains alternate allele and passes filters
-
-		Raises
-		------
-		IndexError
-		"""
-		if not read:
-			return 0
-
-		try:
-
-			offset = pos-read.reference_start
-			if offset<0:
-				raise IndexError
-
-			base_call = read.query_sequence[offset]
-			qual = read.query_qualities[offset]
-
-			if qual<20:
-				raise ReadError(ReadError.ERROR_BASEQ)
-
-			offset_5p = read.reference_end-1-pos if read.is_reverse else pos-read.reference_start
-			if offset_5p<=3:
-				raise ReadError(ReadError.ERROR_5PROXIMITY)
-
-			mm = int(read.get_tag("XM", with_value_type = False))
-
-			if base_call==ref: 
-				if mm>1:
-					raise ReadError(ReadError.ERROR_MISMATCH)
-				return ref
-			elif base_call==alt:
-				if mm>2:
-					raise ReadError(ReadError.ERROR_MISMATCH)
-				return alt
-			else:
-				raise ReadError(ReadError.ERROR_GENOTYPE)
-		
-		# throws when variant pos is outside of read, handle these differently
-		except IndexError as e:
-			return None
-
-	def lookup_allelic(self, chrom, start, end, pos, ref, alt, flip=False):
-		"""
-
-		"""
-		#visited_read_pairs = set()
-
-		tmp_ref_fw = {}
-		tmp_ref_rev = {}
-		tmp_alt_fw = {}
-		tmp_alt_rev = {}
-		tmp_non_fw = {}
-		tmp_non_rev = {}
-
-		reads_ref = []
-		reads_alt = []
-		reads_non = []
+        fw_cutarray = np.array([tmp_fw.get(i, 0.0) for i in range(start, end)])
+        rev_cutarray = np.array([tmp_rev.get(i, 0.0) for i in range(start, end)])
+
+        return {
+            "+": rev_cutarray[::-1] if flip else fw_cutarray, 
+            "-": fw_cutarray[::-1] if flip else rev_cutarray,
+            "fragments": reads
+        }
+
+    def __validate_genotype(self, read, pos, ref, alt):
+        """Validate read genotype
+
+        Parameters
+        ----------
+        read : pysam.AlignedSegment
+            Read to validate
+        pos: int
+            Genomic postion of variant (0-based)
+        ref : str
+            Reference allele
+        alt : str
+            Alternate allele
+
+        Returns
+        -------
+        out : str or bool 
+
+        Raises
+        ------
+        ReadError
+        """
+
+        if not read:
+            return None
+
+        try:
+
+            offset = pos-read.reference_start
+            if offset<0:
+                raise IndexError
+
+            base_call = read.query_sequence[offset]
+            qual = read.query_qualities[offset]
+
+            if qual<20:
+                raise ReadError(ReadError.ERROR_BASEQ)
+
+            offset_5p = read.reference_end-1-pos if read.is_reverse else pos-read.reference_start
+            if offset_5p<=3:
+                raise ReadError(ReadError.ERROR_5PROXIMITY)
+
+            mm = int(read.get_tag("XM", with_value_type = False))
+
+            if base_call==ref: 
+                if mm>1:
+                    raise ReadError(ReadError.ERROR_MISMATCH)
+                return ref
+            elif base_call==alt:
+                if mm>2:
+                    raise ReadError(ReadError.ERROR_MISMATCH)
+                return alt
+            else:
+                raise ReadError(ReadError.ERROR_GENOTYPE)
+        
+        # throws when variant pos is outside of read, handle these differently
+        except IndexError as e:
+            return None
+
+    def lookup_allelic(self, chrom, start, end, pos, ref, alt, flip=False):
+        """Lookup function for allelically resolved counts
+
+        Parameters
+        ----------
+
+        Returns
+        -------
+
+        """
+        #visited_read_pairs = set()
+
+        tmp_ref_fw = {}
+        tmp_ref_rev = {}
+        tmp_alt_fw = {}
+        tmp_alt_rev = {}
+        tmp_non_fw = {}
+        tmp_non_rev = {}
+
+        reads_ref = []
+        reads_alt = []
+        reads_non = []
+        
+        for read1, read2 in self.read_pair_generator(chrom, max(start-10, 0), end+10):
+
+            #try read1 for genotype
+            try:
+                # Check reads for proper genotypes; raise Genotype error
+                # if some thing is problematic, if a read pair doesn't exist or
+                # doesn't overlap the SNV then returns 0
+                read1_gt = self.__validate_genotype(read1, pos, ref, alt)
+                read2_gt = self.__validate_genotype(read2, pos, ref, alt)
+
+                if (not read1_gt) and (not read2_gt): # neither read overlaps SNV
+                    fw = tmp_non_fw
+                    rev = tmp_non_rev
+                    reads = reads_non
+                elif (read1_gt and read2_gt) and read1_gt != read2_gt: # Both reads have a GT, but discordant genotypes
+                    raise ReadError(ReadError.ERROR_GENOTYPE)
+                elif read1_gt == ref or read2_gt == ref: # Matches REF allele
+                    fw = tmp_ref_fw
+                    rev = tmp_ref_rev
+                    reads = reads_ref
+                elif read1_gt == alt or read2_gt == alt: # Matches ALT allele
+                    fw = tmp_alt_fw
+                    rev = tmp_alt_rev
+                    reads = reads_alt
+                else:
+                    raise ReadError()
+
+                if read1:
+                    self.__add_read(read1, fw, rev)
+                if read2:
+                    self.__add_read(read2, fw, rev)					
+
+                if read1:
+                    reads.append(self.__get_fragment(read1))
+                else:
+                    reads.append(self.__get_fragment(read2))
+    
+            except ReadError as e:
+                continue
+    
+        ref_fw_cutarray = np.array([tmp_ref_fw.get(i, 0.0) for i in range(start, end)])
+        ref_rev_cutarray = np.array([tmp_ref_rev.get(i, 0.0) for i in range(start, end)])
+
+        alt_fw_cutarray = np.array([tmp_alt_fw.get(i, 0.0) for i in range(start, end)])
+        alt_rev_cutarray = np.array([tmp_alt_rev.get(i, 0.0) for i in range(start, end)])
+        
+        non_fw_cutarray = np.array([tmp_non_fw.get(i, 0.0) for i in range(start, end)])
+        non_rev_cutarray = np.array([tmp_non_rev.get(i, 0.0) for i in range(start, end)])
+
+        return {
+            ref: {
+                "+": ref_rev_cutarray[::-1] if flip else ref_fw_cutarray, 
+                "-": ref_fw_cutarray[::-1] if flip else ref_rev_cutarray,
+                "fragments": reads_ref
+                },
+            alt: {
+                "+": alt_rev_cutarray[::-1] if flip else alt_fw_cutarray, 
+                "-": alt_fw_cutarray[::-1] if flip else alt_rev_cutarray,
+                "fragments": reads_alt
+            },
+            "other": {
+                "+": non_rev_cutarray[::-1] if flip else non_fw_cutarray, 
+                "-": non_fw_cutarray[::-1] if flip else  non_rev_cutarray,
+                "fragments": reads_non
+            }
+        }
+
+
+    def __getitem__(self, x):
+        """General function to retrieve cutcounts. Currently supports
+        only genomic_intervals and pysam.VariantRecord as input.
+        
+        Parameters
+        ----------
+        x : genome_tools.genomic_interval or pysam.VariantRecord
+            Retrieve cleavages over a windowed region or resolve allelically 
+            over a known variant
+        
+        Returns
+        -------
+        counts : dict
+            A dictionary of read counts resolved to each strand, and 
+            each allele (if approriate)
         
-		for read1, read2 in self.read_pair_generator(chrom, max(start-10, 0), end+10):
-
-			#try read1 for genotype
-			try:
-				# Check reads for proper genotypes; raise Genotype error
-				# if some thing is problematic, if a read pair doesn't exist or
-				# doesn't overlap the SNV then returns 0
-				read1_gt = self.__validate_genotype(read1, pos, ref, alt)
-				read2_gt = self.__validate_genotype(read2, pos, ref, alt)
-
-				if (not read1_gt) and (not read2_gt): # neither read overlaps SNV
-					fw = tmp_non_fw
-					rev = tmp_non_rev
-					reads = reads_non
-				elif (read1_gt and read2_gt) and read1_gt != read2_gt: # Both reads have a GT, but discordant genotypes
-					raise ReadError(ReadError.ERROR_GENOTYPE)
-				elif read1_gt == ref or read2_gt == ref: # Matches REF allele
-					fw = tmp_ref_fw
-					rev = tmp_ref_rev
-					reads = reads_ref
-				elif read1_gt == alt or read2_gt == alt: # Matches ALT allele
-					fw = tmp_alt_fw
-					rev = tmp_alt_rev
-					reads = reads_alt
-				else:
-					raise ReadError()
-
-				if read1:
-					self.__add_read(read1, fw, rev)
-				if read2:
-					self.__add_read(read2, fw, rev)					
-
-				if read1:
-					reads.append(self.__get_fragment(read1))
-				else:
-					reads.append(self.__get_fragment(read2))
-	
-			except ReadError as e:
-				continue
-	
-		ref_fw_cutarray = np.array([tmp_ref_fw.get(i, 0.0) for i in range(start, end)])
-		ref_rev_cutarray = np.array([tmp_ref_rev.get(i, 0.0) for i in range(start, end)])
-
-		alt_fw_cutarray = np.array([tmp_alt_fw.get(i, 0.0) for i in range(start, end)])
-		alt_rev_cutarray = np.array([tmp_alt_rev.get(i, 0.0) for i in range(start, end)])
-		
-		non_fw_cutarray = np.array([tmp_non_fw.get(i, 0.0) for i in range(start, end)])
-		non_rev_cutarray = np.array([tmp_non_rev.get(i, 0.0) for i in range(start, end)])
-
-		return {
-			ref: {
-				"+": ref_rev_cutarray[::-1] if flip else ref_fw_cutarray, 
-				"-": ref_fw_cutarray[::-1] if flip else ref_rev_cutarray,
-				"fragments": reads_ref
-				},
-			alt: {
-				"+": alt_rev_cutarray[::-1] if flip else alt_fw_cutarray, 
-				"-": alt_fw_cutarray[::-1] if flip else alt_rev_cutarray,
-				"fragments": reads_alt
-			},
-			"other": {
-				"+": non_rev_cutarray[::-1] if flip else non_fw_cutarray, 
-				"-": non_fw_cutarray[::-1] if flip else  non_rev_cutarray,
-				"fragments": reads_non
-			}
-		}
-
-
-	def __getitem__(self, x):
-		"""Getter function to retrieve cutcounts
-		
-		Parameters
-		----------
-		x : genome_tools.genomic_interval or pysam.VariantRecord
-		    Retrieve cleavages over a windowed region or resolve allelically over a known variant
-		
-		Returns
-		-------
-		dict
-		    A dictionary of read counts resolved to each strand, and each allele (if approriate)
-		
-		Raises
-		------
-		TypeError
-		    If input is neither a genome_tools.genomic_interval nor a pysam.VariantRecord
-		"""
-		if isinstance(x, genome_tools.genomic_interval):
-			return self.lookup(x)
-		elif isinstance(x, pysam.VariantRecord):
-			return self.lookup_allelic(x)
-		else:
-			raise TypeError()
+        Raises
+        ------
+        TypeError
+            Input is neither a genome_tools.genomic_interval nor 
+            pysam.VariantRecord
+        """
+        if isinstance(x, genome_tools.genomic_interval):
+            return self.lookup(x)
+        elif isinstance(x, pysam.VariantRecord):
+            return self.lookup_allelic(x)
+        else:
+            raise TypeError(f"Query type not supported: {type(x)}")
diff --git a/footprint_tools/data/__init__.py b/footprint_tools/data/__init__.py
new file mode 100644
index 0000000..7c86dc8
--- /dev/null
+++ b/footprint_tools/data/__init__.py
@@ -0,0 +1,5 @@
+"""Module to facilitate the loading a procesing of data. Largely
+inspired/copied from Kipoi (http://github.com/kipoi) and pyTorch
+data loaders"""
+
+__all__ = ["dataset", "loader", "sampler", "utils"]
\ No newline at end of file
diff --git a/footprint_tools/data/process.py b/footprint_tools/data/process.py
new file mode 100644
index 0000000..609c9a1
--- /dev/null
+++ b/footprint_tools/data/process.py
@@ -0,0 +1,46 @@
+from footprint_tools.data.processor import data_processor
+
+import logging
+logger = logging.getLogger(__name__)
+
+class base_process(object):
+    def batch_iter(self, **kwargs):
+        raise  NotImplementedError
+
+    def load_all(self, **kwargs):
+        """Load an entire dataset"""
+        return [x for x in self.batch_iter(**kwargs)]
+        
+class process(base_process):
+    """All datasets should subclass this class. All subclases should
+    override `__len__` and `__getitem__` to support integer indexing"""
+
+    def __getitem__(self, index):
+        """Return one sample"""
+        raise NotImplementedError
+    
+    def __len__(self):
+        """Return number of all elements"""
+        raise NotImplementedError
+
+    def _batch_iterable(self, batch_size=1, num_workers=0, **kwargs):
+        
+        dp = data_processor(self, 
+                        batch_size=batch_size,
+                        num_workers=num_workers,
+                        **kwargs)
+
+        return dp
+
+    def batch_iter(self, batch_size=1, num_workers=0, **kwargs):
+        
+        dp =  self._batch_iterable(batch_size=batch_size,
+                                    num_workers=num_workers,
+                                    **kwargs)
+        
+        return iter(dp)
+
+    # Example of overloaded `load_all` function
+    # def load_all(self, batch_size=1, **kwargs):
+    #     """Load all data"""
+    #     return numpy_collate_concat([x for x in tqdm(self.batch_iter(batch_size, **kwargs))])
diff --git a/footprint_tools/data/processor.py b/footprint_tools/data/processor.py
new file mode 100644
index 0000000..caf6743
--- /dev/null
+++ b/footprint_tools/data/processor.py
@@ -0,0 +1,161 @@
+import sys
+import multiprocessing
+import threading
+import traceback
+
+from footprint_tools.data.sampler import sequential_sampler, minibatch_sampler
+from footprint_tools.data.utils import list_collate
+
+import logging
+logger = logging.getLogger(__name__)
+
+default_collate = list_collate
+
+class exception_wrapper(object):
+    """Wraps an exception plus traceback to communicate across threads"""
+
+    def __init__(self, exc_info):
+        self.exc_type = exc_info[0]
+        self.exc_msg = "".join(traceback.format_exception(*exc_info))
+
+def _worker_loop(process, index_queue, data_queue, collate_fn):
+    """Worker loop to load data"""
+    while True:
+        r = index_queue.get()
+        if r is None:
+            break
+        idx, batch_indices = r
+
+        try:
+            batch = collate_fn([process[i] for i in batch_indices])
+        except Exception:
+            data_queue.put((idx, exception_wrapper(sys.exc_info())))
+        else:
+            data_queue.put((idx, batch))
+
+    # clean up function that closes filehandlers, etc.
+    if hasattr(process, 'cleanup'):
+        process.cleanup()
+
+class data_processor_iter(object):
+    def __init__(self, processor):
+        self.process = processor.process
+        self.collate_fn = processor.collate_fn
+        self.batch_sampler = processor.batch_sampler
+        self.num_workers = processor.num_workers
+        self.done_event = threading.Event()
+
+        self.sample_iter = iter(self.batch_sampler)
+                
+        if self.num_workers > 0:
+
+            logger.info(f"Using {self.num_workers} threads to process data")
+
+            self.index_queue = multiprocessing.SimpleQueue()
+            self.data_queue = multiprocessing.SimpleQueue()
+            self.batches_outstanding = 0
+            self.shutdown = False
+            self.send_idx = 0
+            self.rcvf_idx = 0
+            self.reorder_dict = {}
+
+            self.workers = [
+                multiprocessing.Process(
+                    target=_worker_loop,
+                    args=(self.process, self.index_queue, self.data_queue, self.collate_fn))
+                for _ in range(self.num_workers)
+            ]
+
+            for w in self.workers:
+                w.daemon = True
+                w.start()
+
+            for _ in range(2 * self.num_workers):
+                self._put_indices()
+        else:
+
+            logger.info(f"Using a single threads to process data")
+
+    
+    def __len__(self):
+        return len(self.batch_sampler)
+
+    def __next__(self):
+        if self.num_workers == 0:
+            indices = next(self.sample_iter)
+            batch = self.collate_fn([self.process[i] for i in indices])
+            return batch
+        
+        if self.rcvf_idx in self.reorder_dict:
+            batch = self.reorder_dict.pop(self.rcvf_idx)
+            return self._process_next_batch(batch)
+        
+        if self.batches_outstanding == 0:
+            self._shutdown_workers()
+            raise StopIteration
+
+        while True:
+            assert (not self.shutdown and self.batches_outstanding > 0)
+            idx, batch = self.data_queue.get()
+            self.batches_outstanding -= 1
+            if idx != self.rcvf_idx:
+                # store out of order samples
+                self.reorder_dict[idx] = batch
+                continue
+            return self._process_next_batch(batch)
+            
+    def __iter__(self):
+        return self
+
+    def _put_indices(self):
+        assert self.batches_outstanding < 2 * self.num_workers
+        indices = next(self.sample_iter, None)
+        if indices == None:
+            return
+        self.index_queue.put((self.send_idx, indices))
+        self.batches_outstanding += 1
+        self.send_idx += 1
+
+    def _process_next_batch(self, batch):
+        self.rcvf_idx += 1
+        self._put_indices()
+        if isinstance(batch, exception_wrapper):
+            raise batch.exc_type(batch.exc_msg)
+        return batch
+
+    def __getstate__(self):
+        raise NotImplementedError("Class cannot be pickled")
+
+    def _shutdown_workers(self):
+        if not self.shutdown:
+            self.shutdown = True
+            self.done_event.set()
+            for _ in self.workers:
+                self.index_queue.put(None)
+                # wait for worker to join?
+
+    def __del__(self):
+        if self.num_workers > 0:
+            self._shutdown_workers()
+
+class data_processor(object):
+    def __init__(self, process, batch_size=1, num_workers=0, collate_fn=default_collate):
+        self.process = process
+        self.batch_size = batch_size
+        self.num_workers = num_workers
+        self.collate_fn = collate_fn
+
+        sampler = sequential_sampler(process)
+        batch_sampler = minibatch_sampler(sampler, batch_size)
+
+        self.sampler = sampler
+        self.batch_sampler = batch_sampler
+
+        logger.info(f"Processing data with batch_size = {self.batch_size:,} resulting in {len(self):,} batches" )
+        logger.info(f"Using '{self.collate_fn.__doc__}' to collate batch chunks")
+
+    def __iter__(self):
+        return data_processor_iter(self)
+
+    def __len__(self):
+        return len(self.batch_sampler)
diff --git a/footprint_tools/data/sampler.py b/footprint_tools/data/sampler.py
new file mode 100644
index 0000000..6ede369
--- /dev/null
+++ b/footprint_tools/data/sampler.py
@@ -0,0 +1,39 @@
+class base_sampler(object):
+    def __init__(self, data_source):
+        pass
+
+    def __iter__(self):
+        raise NotImplementedError
+
+    def __len__(self):
+        raise NotImplementedError
+
+class sequential_sampler(base_sampler):
+    def __init__(self, data_source):
+        self.data_source = data_source
+    
+    def __iter__(self):
+        return iter(range(len(self.data_source)))
+
+    def __len__(self):
+        return len(self.data_source)
+
+class minibatch_sampler(base_sampler):
+    def __init__(self, sampler, batch_size):
+        self.sampler = sampler
+        self.batch_size = batch_size
+
+    def __iter__(self):
+        batch = []
+        for idx in self.sampler:
+            batch.append(idx)
+            if len(batch) == self.batch_size:
+                yield batch
+                batch = []
+        if len(batch) > 0:
+            yield batch
+    
+    def __len__(self):
+        return (len(self.sampler) + self.batch_size - 1) // self.batch_size
+        
+
diff --git a/footprint_tools/data/utils.py b/footprint_tools/data/utils.py
new file mode 100644
index 0000000..cf00434
--- /dev/null
+++ b/footprint_tools/data/utils.py
@@ -0,0 +1,59 @@
+import numpy as np
+import collections
+import itertools
+
+string_classes = (str, bytes)
+
+def _list_collate(stack_fn=list):
+    def list_collate_fn(batch):
+        if isinstance(batch[0], collections.Mapping):
+            return {key: list_collate_fn([d[key] for d in batch]) for key in batch[0]}
+        elif isinstance(batch[0], collections.Sequence):
+            transposed = zip(*batch)
+            return [list_collate_fn(samples) for samples in transposed]
+        else:
+            return stack_fn(batch)
+
+    return list_collate_fn
+
+list_collate = _list_collate(list)
+list_collate.__doc__ = __name__ + '.list_collate'
+
+list_collate_concat = _list_collate(lambda x: list(itertools.chain(*x)))
+list_collate_concat.__doc__ = __name__ + '.list_collate_concat'
+
+def _numpy_collate(stack_fn=np.stack):
+    def numpy_collate_fn(batch):
+        "Puts each data field into a tensor with outer dimension batch size"
+        if type(batch[0]).__module__ == 'numpy':
+            elem = batch[0]
+            if type(elem).__name__ == 'ndarray':
+                return stack_fn(batch, 0)
+            if elem.shape == ():  # scalars
+                return np.array(batch)
+        elif isinstance(batch[0], int):
+            return np.asarray(batch)
+        elif isinstance(batch[0], float):
+            return np.asarray(batch)
+        elif batch[0] is None:
+            return np.asarray(batch)
+        elif isinstance(batch[0], string_classes):
+            # Also convert to a numpy array
+            return np.asarray(batch)
+            # return batch
+        elif isinstance(batch[0], collections.Mapping):
+            return {key: numpy_collate_fn([d[key] for d in batch]) for key in batch[0]}
+        elif isinstance(batch[0], collections.Sequence):
+            transposed = zip(*batch)
+            return [numpy_collate_fn(samples) for samples in transposed]
+
+        raise TypeError(("batch must contain tensors, numbers, dicts or lists; found {}"
+                         .format(type(batch[0]))))
+
+    return numpy_collate_fn
+
+numpy_collate = _numpy_collate(np.stack)
+numpy_collate.__doc__ = __name__ + '.numpy_collate'
+
+numpy_collate_concat = _numpy_collate(np.concatenate)
+numpy_collate_concat.__doc__ = __name__ + '.numpy_collate_concat'
\ No newline at end of file
diff --git a/footprint_tools/logging.conf b/footprint_tools/logging.conf
new file mode 100644
index 0000000..df960d9
--- /dev/null
+++ b/footprint_tools/logging.conf
@@ -0,0 +1,21 @@
+[loggers]
+keys=root
+
+[logger_root]
+handlers=stream
+level=INFO
+
+[formatters]
+keys=color
+
+[formatter_color]
+class=colorlog.ColoredFormatter
+format=%(bg_black)s%(asctime)s%(reset)s %(log_color)s%(levelname)-8s%(reset)s %(bold_black)s%(name)s:%(reset)s%(message)s
+datefmt=%m-%d %H:%M:%S
+
+[handlers]
+keys=stream
+
+[handler_stream]
+class=StreamHandler
+formatter=color
diff --git a/footprint_tools/modeling/bias.py b/footprint_tools/modeling/bias.py
index 4ce9ab0..96aaf84 100644
--- a/footprint_tools/modeling/bias.py
+++ b/footprint_tools/modeling/bias.py
@@ -1,101 +1,102 @@
-# Copyright 2015 Jeff Vierstra
+"""
+"""
+
 import numpy as np
 import itertools
 import random
 
-
 class bias_model(object):
 
-	def __init__(self):
+    def __init__(self):
 
-		self.model = {}
-		
-		self.k = 6
-		self.mid = 3
+        self.model = {}
+        
+        self.k = 6
+        self.mid = 3
 
-	def __getitem__(self, key):
-		
-		return self.model.get(key, 1e-6)
+    def __getitem__(self, key):
+        
+        return self.model.get(key, 1e-6)
 
-	def __setitem__(self, key, value):
+    def __setitem__(self, key, value):
 
-		self.model[key] = value
+        self.model[key] = value
 
-	def offset(self):
+    def offset(self):
 
-		return max(self.k-self.mid, self.mid)
+        return max(self.k-self.mid, self.mid)
 
-	def shuffle(self):
-		"""Randomly shuffle the bias model
+    def shuffle(self):
+        """Randomly shuffle the bias model
 
-		:returns: 
-		"""
+        :returns: 
+        """
 
-		ret = bias_model()
-		ret.model = { x: y for (x, y) in zip(list(self.model.keys()), sorted(list(self.model.values()), key = lambda k: random.random())) }
-		ret.offset = self.offset
-		return ret
+        ret = bias_model()
+        ret.model = { x: y for (x, y) in zip(list(self.model.keys()), sorted(list(self.model.values()), key = lambda k: random.random())) }
+        ret.offset = self.offset
+        return ret
 
-	def predict(self, probs, n = 100):
-		"""Compute cleavage propensities from sequence
-		
-		:param probs: an array of probilities (relative values)
-		:type probs: numpy array (float64)
-		:param n: number tags to distrbute
-		:type n: integer
+    def predict(self, probs, n = 100):
+        """Compute cleavage propensities from sequence
+        
+        :param probs: an array of probilities (relative values)
+        :type probs: numpy array (float64)
+        :param n: number tags to distrbute
+        :type n: integer
 
-		:returns: integer array of relative cleavage counts
-		"""
+        :returns: integer array of relative cleavage counts
+        """
 
-		return np.around( probs / np.sum(probs) * n )
+        return np.around( probs / np.sum(probs) * n )
 
 class kmer_model(bias_model):
 
-	def __init__(self, filepath):
-		
-		bias_model.__init__(self)		
-		self.read_model(filepath)
-	
-	def read_model(self, filepath):
-		"""Read the k-mer model from a file.
+    def __init__(self, filepath):
+        
+        bias_model.__init__(self)		
+        self.read_model(filepath)
+    
+    def read_model(self, filepath):
+        """Read the k-mer model from a file.
 
-		:param filepath: the path to a 6-model
-		:type filepath: string
+        :param filepath: the path to a 6-model
+        :type filepath: string
 
-		:returns:
-		"""
+        :returns:
+        """
 
-		try:
+        try:
 
-			for line in open(filepath, 'r'):
-				(seq, prob) = line.strip().split('\t')
-				self.model[seq.upper()] = float(prob)
+            for line in open(filepath, 'r'):
+                (seq, prob) = line.strip().split('\t')
+                self.model[seq.upper()] = float(prob)
 
-		except IOError:
-			
-			raise IOError("Cannot open file: %s" % filepath)
-	
-	def probs(self, seq):
-		"""Generate cleavage preference array from DNA sequence
+        except IOError:
+            
+            raise IOError("Cannot open file: %s" % filepath)
+    
+    def probs(self, seq):
+        """Generate cleavage preference array from DNA sequence
 
-		:param seq: DNA sequence
+        :param seq: DNA sequence
 
-		:returns: array (float)
-		"""
-		k = self.k
-		mid = self.mid
-		ofst = self.offset()
+        :returns: array (float)
+        """
+        k = self.k
+        mid = self.mid
+        offset = self.offset()
 
-		return np.array([ self.__getitem__( seq[(i-mid):(i+(k-mid))] ) for i in range(ofst, len(seq)-ofst) ], dtype = np.float64)
+        return np.array([ self.__getitem__( seq[(i-mid):(i+(k-mid))] ) for i in range(offset, len(seq)-offset) ], dtype = np.float64)
 
 class uniform_model(bias_model):
 
-	def __init__(self):
+    def __init__(self):
 
-		bias_model.__init__(self)
+        bias_model.__init__(self)
 
-		for seq in itertools.product('ATCG', repeat = self.k):
-			self.model[''.join(seq)] = 1.0
+        for seq in itertools.product('ATCG', repeat = self.k):
+            self.model[''.join(seq)] = 1.0
 
-	def probs(self, seq):
-		return np.ones( len(seq) )
+    def probs(self, seq):
+        return np.ones( len(seq) )
diff --git a/footprint_tools/modeling/dispersion.pxd b/footprint_tools/modeling/dispersion.pxd
index add2d57..650c2ed 100644
--- a/footprint_tools/modeling/dispersion.pxd
+++ b/footprint_tools/modeling/dispersion.pxd
@@ -1,37 +1,35 @@
+# Copyright 2015-2021 Jeff Vierstra
+
 cimport numpy as np
 
 ctypedef np.float64_t data_type_t
 
-
 cpdef data_type_t piecewise_four(x, data_type_t x0, data_type_t x1, data_type_t x2, data_type_t x3, 
-	 								data_type_t y0, data_type_t y1, data_type_t y2, data_type_t y3,
-	 								data_type_t k0, data_type_t k1, data_type_t k2, data_type_t k3)
+                                     data_type_t y0, data_type_t y1, data_type_t y2, data_type_t y3,
+                                     data_type_t k0, data_type_t k1, data_type_t k2, data_type_t k3)
 
 
 cpdef data_type_t piecewise_three(x, data_type_t x0, data_type_t x1, data_type_t x2,
-	 								data_type_t y0, data_type_t y1, data_type_t y2,
-	 								data_type_t k0, data_type_t k1, data_type_t k2)
-
-# cpdef data_type_t piecewise_three(x, data_type_t x0, data_type_t y0, data_type_t x1, data_type_t k1, data_type_t k2, data_type_t k3)
-# cpdef data_type_t piecewise_three_fast(data_type_t x, data_type_t x0, data_type_t y0, data_type_t x1, data_type_t k1, data_type_t k2, data_type_t k3) nogil
-
+                                     data_type_t y0, data_type_t y1, data_type_t y2,
+                                     data_type_t k0, data_type_t k1, data_type_t k2)
 
 cdef class dispersion_model:
-	
-	cdef np.ndarray _h
-	cdef np.ndarray _p
-	cdef np.ndarray _r
-	cdef np.ndarray _r_params
-	cdef np.ndarray _mu_params
+    
+    cdef np.ndarray _h
+    cdef np.ndarray _p
+    cdef np.ndarray _r
+    cdef np.ndarray _r_params
+    cdef np.ndarray _mu_params
 
-	cpdef data_type_t fit_mu(self, data_type_t x)
-	cpdef data_type_t fit_r(self, data_type_t x)
+    cpdef data_type_t fit_mu(self, data_type_t x)
+    cpdef data_type_t fit_r(self, data_type_t x)
 
-	cpdef data_type_t [:] log_pmf_values(self, data_type_t [:] exp, data_type_t [:] obs)
-	cpdef data_type_t [:] pmf_values(self, data_type_t [:] exp, data_type_t [:] obs)
+    cpdef data_type_t [:] log_pmf_values(self, data_type_t [:] exp, data_type_t [:] obs)
+    cpdef data_type_t [:] pmf_values(self, data_type_t [:] exp, data_type_t [:] obs)
 
-	cpdef data_type_t [:] log_pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res)
-	cpdef data_type_t [:] pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res)
+    cpdef data_type_t [:] log_pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res)
+    cpdef data_type_t [:] pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res)
 
-	cpdef data_type_t [:] p_values(self, data_type_t [:] exp, data_type_t [:] obs)
-	cpdef resample_p_values(self, data_type_t [:] x, int times)
\ No newline at end of file
+    cpdef data_type_t [:] p_values(self, data_type_t [:] exp, data_type_t [:] obs)
+    cpdef sample(self, data_type_t [:] x, int times)
+    
\ No newline at end of file
diff --git a/footprint_tools/modeling/dispersion.pyx b/footprint_tools/modeling/dispersion.pyx
index d1e7d9b..23b1caa 100644
--- a/footprint_tools/modeling/dispersion.pyx
+++ b/footprint_tools/modeling/dispersion.pyx
@@ -1,39 +1,41 @@
-# Copyright 2015-2021 Jeff Vierstra
+"""
+This module contains classess and functions that implement a dispersion model.
+"""
 
 # cython: boundscheck=False
 # cython: wraparound=False
 # cython: nonecheck=False
 # cython: embedsignature=True
 
-import warnings
-
-from footprint_tools.stats.distributions import nbinom
-import numpy as np
-from scipy import optimize
-
-import pwlf
-
 cimport cython
 cimport numpy as np
+import numpy as np
 
 from footprint_tools.stats.distributions cimport nbinom
+from footprint_tools.stats.distributions import nbinom
 
+import warnings
+from scipy import optimize
+import pwlf
+
+import logging
+logger = logging.getLogger(__name__)
 
 cpdef data_type_t piecewise_three(x, data_type_t x0, data_type_t x1, data_type_t x2,
-	 								data_type_t y0, data_type_t y1, data_type_t y2,
-	 								data_type_t k0, data_type_t k1, data_type_t k2):
-	
-	return (
+                                     data_type_t y0, data_type_t y1, data_type_t y2,
+                                     data_type_t k0, data_type_t k1, data_type_t k2):
+    
+    return (
         (x<x0) * (y0 + k0*(x)) +
         ((x>=x0) & (x<x1)) * (y1 + k1*(x)) +
         (x>=x1) * (y2 + k2*(x))
     )
 
 cpdef data_type_t piecewise_four(x, data_type_t x0, data_type_t x1, data_type_t x2, data_type_t x3, 
-	 								data_type_t y0, data_type_t y1, data_type_t y2, data_type_t y3,
-	 								data_type_t k0, data_type_t k1, data_type_t k2, data_type_t k3):
-	
-	return (
+                                     data_type_t y0, data_type_t y1, data_type_t y2, data_type_t y3,
+                                     data_type_t k0, data_type_t k1, data_type_t k2, data_type_t k3):
+    
+    return (
         (x<x0) * (y0 + k0*(x)) +
         ((x>=x0) & (x<x1)) * (y1 + k1*(x)) +
         ((x>=x1) & (x<x2)) * (y2 + k2*(x)) +
@@ -41,10 +43,10 @@ cpdef data_type_t piecewise_four(x, data_type_t x0, data_type_t x1, data_type_t
     )
 
 cpdef data_type_t piecewise_five(x, data_type_t x0, data_type_t x1, data_type_t x2, data_type_t x3, data_type_t x4, 
-	 								data_type_t y0, data_type_t y1, data_type_t y2, data_type_t y3, data_type_t y4,
-	 								data_type_t k0, data_type_t k1, data_type_t k2, data_type_t k3, data_type_t k4):
-	
-	return (
+                                     data_type_t y0, data_type_t y1, data_type_t y2, data_type_t y3, data_type_t y4,
+                                     data_type_t k0, data_type_t k1, data_type_t k2, data_type_t k3, data_type_t k4):
+    
+    return (
         (x<x0) * (y0 + k0*(x)) +
         ((x>=x0) & (x<x1)) * (y1 + k1*(x)) +
         ((x>=x1) & (x<x2)) * (y2 + k2*(x)) +
@@ -53,367 +55,470 @@ cpdef data_type_t piecewise_five(x, data_type_t x0, data_type_t x1, data_type_t
     )
 
 cdef class dispersion_model:
-	"""
-	Dispersion model base class
-	"""
-
-	def __init__(self):
-		
-		# Histograms
-		self._h = None
-		
-		# Emperical fits
-		self._p = None
-		self._r = None
-
-		self._mu_params = self._r_params = None
-
-		pass
-
-	# Pickling function
-	def __reduce__(self):
-		x = {}
-		x['mu_params'] = self.mu_params
-		x['r_params'] = self.r_params
-		return (dispersion_model, (), x)
-
-	# Pickling function
-	def __setstate__(self, x):
-		self.mu_params = x['mu_params']
-		self.r_params = x['r_params']
-
-	property h:
-		"""Histrogram of observed cleavages at each predicted cleavage rate"""
-		def __get__(self):
-			return self._h
-		def __set__(self, x):
-			self._h = x
-
-	property p:
-		"""Array of the negative binomial ML fit parameter `p`"""
-		def __get__(self):
-			return self._p
-		def __set__(self, x):
-			self._p = x
-
-	property r:
-		"""Array of the negative binomial ML fit parameter `r`"""
-		def __get__(self):
-			return self._r
-		def __set__(self, x):
-			self._r = x
-
-	property mu_params:
-		
-		def __get__(self):
-			return self._mu_params
-		def __set__(self, x):
-			self._mu_params = np.array(x, order = 'c')
-
-	property r_params:
-		def __get__(self):
-			return self._r_params
-		def __set__(self, x):
-			self._r_params = np.array(x, order = 'c')
-
-
-	cpdef data_type_t fit_mu(self, data_type_t x):
-		"""Computes the mu term for the negative binomial.
-		
-		:param x: predicted values to be converted to observed means
-		
-		:returns mu: (float)
-		"""
-
-		cdef data_type_t [:] par = self._mu_params
-		cdef data_type_t res = piecewise_three(x, *par)
-
-		return res if res > 0.0 else 0.1
-
-	cpdef data_type_t fit_r(self, data_type_t x):
-		"""Computes the dispersion term for the negative binomial.
-		Note that the model parameters estimate the inverse.
-
-		:param x: predicted values to be converted to dispersion parameter r
-
-		:returns: A fitted paramter for the negative bionomial distribution
-		:rtype: float
-		"""
-
-		cdef data_type_t [:] par = self._r_params
-		cdef data_type_t res = 1.0/piecewise_five(x, *par)
-
-		return res if res > 0.0 else 1e-6
-
-	def __str__(self):
-		"""Print model to string"""
-		res = "Oops...the functionality has not yet been implemented!"
-		return res
-
-	@cython.cdivision(True)
-	cpdef data_type_t [:] log_pmf_values(self, data_type_t [:] exp, data_type_t [:] obs):
-		"""Computing the probability mass function"""
-
-		cdef int i, n = exp.shape[0]
-		cdef double r, mu
-		
-		cdef data_type_t [:] res = np.zeros(n, dtype = np.float64, order = 'c')
-	
-		for i in range(n):
-			r = self.fit_r(exp[i])
-			mu = self.fit_mu(exp[i])
-			res[i] = nbinom.logpmf(<int>obs[i], (r/(r+mu)), r)
-
-		return res
-
-	@cython.cdivision(True)
-	cpdef data_type_t [:] pmf_values(self, data_type_t [:] exp, data_type_t [:] obs):
-		"""Computing the probability mass function"""
-
-		cdef int i, n = exp.shape[0]
-		cdef double r, mu
-		
-		cdef data_type_t [:] res = np.zeros(n, dtype = np.float64, order = 'c')
-
-		for i in range(n):
-			r = self.fit_r(exp[i])
-			mu = self.fit_mu(exp[i])
-			res[i] = nbinom.pmf(<int>obs[i], (r/(r+mu)), r)
-
-		return res
-
-
-	@cython.cdivision(True)
-	cpdef data_type_t [:] pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res):
-		"""Computing the probability mass function (but writes to a matrix pointer)"""
-
-		cdef int i, n = exp.shape[0]
-		cdef double r, mu
-		
-		#cdef data_type_t [:] res = np.zeros(n, dtype = np.float64, order = 'c')
-
-		for i in range(n):
-			r = self.fit_r(exp[i])
-			mu = self.fit_mu(exp[i])
-			res[i] = nbinom.pmf(<int>obs[i], (r/(r+mu)), r)
-
-		return res
-
-	@cython.cdivision(True)
-	cpdef data_type_t [:] log_pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res):
-		"""Computing the probability mass function (but writes to a matrix pointer)"""
-
-		cdef int i, n = exp.shape[0]
-		cdef double r, mu
-		
-		#cdef data_type_t [:] res = np.zeros(n, dtype = np.float64, order = 'c')
-
-		for i in range(n):
-			r = self.fit_r(exp[i])
-			mu = self.fit_mu(exp[i])
-			res[i] = nbinom.logpmf(<int>obs[i], (r/(r+mu)), r)
-
-		return res
-
-	cpdef data_type_t [:] p_values(self, data_type_t [:] exp, data_type_t [:] obs):
-		"""Computes CDF p-value from negative binomial
-
-		:param exp: 
-		:param obs:
-
-		:return p:
-		"""
-		
-		cdef int i, n = exp.shape[0]
-		cdef double r, mu
-		cdef data_type_t [:] res = np.ones(n, dtype = np.float64, order = 'c')
-
-		for i in range(n):
-			r = self.fit_r(exp[i])
-			mu = self.fit_mu(exp[i])
-			res[i] = nbinom.cdf(<int>obs[i], r/(r+mu), r)
-
-		return res
-
-	cpdef resample_p_values(self, data_type_t [:] x, int times):
-
-		cdef int i, j, n = x.shape[0]
-		cdef data_type_t r, mu, k
-
-		cdef long [:] vals
-
-		cdef data_type_t [:,:] res_pvals = np.ones((n, times), dtype = np.float64, order = 'c')
-		cdef long [:,:] res_vals = np.zeros((n, times), dtype = np.int_, order = 'c')
-
-		for i in range(n):
-			r = self.fit_r(x[i])
-			mu = self.fit_mu(x[i])
-		
-			vals = np.random.negative_binomial(r, r/(r+mu), times)
-			
-			for j in range(times):
-				res_pvals[i, j] = nbinom.cdf(<int>vals[j], r/(r+mu), r)
-				res_vals[i, j] = vals[j]
-
-		return res_vals, res_pvals
-
-
-@cython.boundscheck(True)
-@cython.wraparound(True)
-
+    """
+    Dispersion model class
+    """
+    def __init__(self):
+        # Histograms
+        self._h = None
+        
+        # Emperical fits
+        self._p = None
+        self._r = None
+
+        self._mu_params = self._r_params = None
+
+    # Pickling function
+    def __reduce__(self):
+        x = {}
+        x['mu_params'] = self.mu_params
+        x['r_params'] = self.r_params
+        return (dispersion_model, (), x)
+
+    # Pickling function
+    def __setstate__(self, x):
+        self.mu_params = x['mu_params']
+        self.r_params = x['r_params']
+
+    property h:
+        """Histrogram of observed cleavages at each predicted cleavage rate"""
+        def __get__(self):
+            return self._h
+        def __set__(self, x):
+            self._h = x
+
+    property p:
+        """Array of the negative binomial MLE fit parameters `p`"""
+        def __get__(self):
+            return self._p
+        def __set__(self, x):
+            self._p = x
+
+    property r:
+        """Array of the negative binomial MLE fit parameters `r`"""
+        def __get__(self):
+            return self._r
+        def __set__(self, x):
+            self._r = x
+
+    property mu_params:
+        def __get__(self):
+            return self._mu_params
+        def __set__(self, x):
+            self._mu_params = np.array(x, order = 'c')
+
+    property r_params:
+        def __get__(self):
+            return self._r_params
+        def __set__(self, x):
+            self._r_params = np.array(x, order = 'c')
+
+
+    cpdef data_type_t fit_mu(self, data_type_t x):
+        """Computes the fitted mu term for the negative binomial
+        from a piece-wise linear fit.
+        
+        Parameters
+        ----------
+        x : float
+
+        Returns
+        -------
+        mu : float
+            mu computed from the regression fit
+        """
+
+        cdef data_type_t [:] par = self._mu_params
+        cdef data_type_t res = piecewise_three(x, *par)
+
+        return res if res > 0.0 else 0.1
+
+    cpdef data_type_t fit_r(self, data_type_t x):
+        """Computes the dispersion term for the negative binomial
+        from a piece-wise linear fit. Note that the model parameters 
+        estimate the inverse.
+        
+        Parameters
+        ----------
+        x : float
+        
+        Returns
+        -------
+        r : float
+            r computed from the regression fit
+        """
+        cdef data_type_t [:] par = self._r_params
+        cdef data_type_t res = 1.0/piecewise_five(x, *par)
+
+        return res if res > 0.0 else 1e-6
+
+    def __str__(self):
+        """Print model to string"""
+        raise NotImplementedError
+
+    @cython.cdivision(True)
+    cpdef data_type_t [:] log_pmf_values(self, data_type_t [:] exp, data_type_t [:] obs):
+        """Compute the log probability mass function
+
+        Parameters
+        ----------
+        exp: ndarray
+            Expected cleavage counts
+        obs: ndarray
+            Observed cleavage counts
+
+        Returns
+        -------
+        logp : ndarray
+            Array of log probability mass function values computed from 
+            the expected cleavage distributions
+        """
+        cdef int i, n = exp.shape[0]
+        cdef double r, mu
+        
+        cdef data_type_t [:] res = np.zeros(n, dtype = np.float64, order = 'c')
+    
+        for i in range(n):
+            r = self.fit_r(exp[i])
+            mu = self.fit_mu(exp[i])
+            res[i] = nbinom.logpmf(<int>obs[i], (r/(r+mu)), r)
+        
+        return res
+
+    @cython.cdivision(True)
+    cpdef data_type_t [:] pmf_values(self, data_type_t [:] exp, data_type_t [:] obs):
+        """Compute the probability mass function
+
+        Parameters
+        ----------
+        exp: ndarray
+            Expected cleavage counts
+        obs: ndarray
+            Observed cleavage counts
+
+        Returns
+        -------
+        p : ndarray
+            Array of probability mass function values computed from 
+            the expected cleavage distributions
+        """
+        cdef int i, n = exp.shape[0]
+        cdef double r, mu
+        
+        cdef data_type_t [:] res = np.zeros(n, dtype = np.float64, order = 'c')
+
+        for i in range(n):
+            r = self.fit_r(exp[i])
+            mu = self.fit_mu(exp[i])
+            res[i] = nbinom.pmf(<int>obs[i], (r/(r+mu)), r)
+
+        return res
+
+    @cython.cdivision(True)
+    cpdef data_type_t [:] log_pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res):
+        """Computing the log probability mass function
+
+        Parameters
+        ----------
+        exp: ndarray
+            Expected cleavage counts
+        obs: ndarray
+            Observed cleavage counts
+
+        Returns
+        -------
+        logp : ndarray (memoryview)
+            Array pointer to log probability mass function values computed from 
+            the expected cleavage distributions
+
+        Notes
+        -----
+        This function is equivalent to `log_pmf_values`, except it stores values to
+        a matrix pointer
+        """
+        cdef int i, n = exp.shape[0]
+        cdef double r, mu
+        
+        for i in range(n):
+            r = self.fit_r(exp[i])
+            mu = self.fit_mu(exp[i])
+            res[i] = nbinom.logpmf(<int>obs[i], (r/(r+mu)), r)
+
+        return res
+    
+    @cython.cdivision(True)
+    cpdef data_type_t [:] pmf_values_0(self, data_type_t [:] exp, data_type_t [:] obs, data_type_t [:] res):
+        """Compute the probability mass function
+
+        Parameters
+        ----------
+        exp: ndarray
+            Expected cleavage counts
+        obs: ndarray
+            Observed cleavage counts
+
+        Returns
+        -------
+        p : ndarray (memoryview)
+            Array pointer to probability mass function values computed from 
+            the expected cleavage distributions
+
+        Notes
+        -----
+        This function is equivalent to `pmf_values`, except it stores values to
+        a matrix pointer
+        """
+        cdef int i, n = exp.shape[0]
+        cdef double r, mu
+        
+        for i in range(n):
+            r = self.fit_r(exp[i])
+            mu = self.fit_mu(exp[i])
+            res[i] = nbinom.pmf(<int>obs[i], (r/(r+mu)), r)
+
+        return res
+
+    cpdef data_type_t [:] p_values(self, data_type_t [:] exp, data_type_t [:] obs):
+        """Compute cumulative distribution (lower-tail p-value)
+        from negative binomial
+
+        Parameters
+        ----------
+        exp: ndarray
+            Expected cleavage counts
+        obs: ndarray
+            Observed cleavage counts
+
+        Returns
+        -------
+        pvals : ndarray
+            Array of p-values
+        """     
+        cdef int i, n = exp.shape[0]
+        cdef double r, mu
+        cdef data_type_t [:] res = np.ones(n, dtype = np.float64, order = 'c')
+
+        for i in range(n):
+            r = self.fit_r(exp[i])
+            mu = self.fit_mu(exp[i])
+            res[i] = nbinom.cdf(<int>obs[i], r/(r+mu), r)
+
+        return res
+
+    cpdef sample(self, data_type_t [:] x, int times):
+        """Sample counts from negative binomial distribution and
+        compute p-values
+
+        Parameters
+        ----------
+        x : ndarray
+            Count values to specifying from which distribution
+            to resample. This typically expected count values.
+        times : int
+            Number of times to sample (per element)
+
+        Returns
+        -------
+        sampled_counts : ndarray
+            Array of sample counts (2-D array - positions by number of samples)
+        sampled_pvals : ndarray
+            Array of sample counts (2-D array - positions by number of samples)
+        """
+        cdef int i, j, n = x.shape[0]
+        cdef data_type_t r, mu, k
+
+        cdef long [:,:] sampled_vals = np.zeros((n, times), dtype = np.int_, order = 'c')
+        cdef data_type_t [:,:] sampled_pvals = np.ones((n, times), dtype = np.float64, order = 'c')
+
+        cdef long [:] vals
+
+        for i in range(n):
+            r = self.fit_r(x[i])
+            mu = self.fit_mu(x[i])
+        
+            vals = np.random.negative_binomial(r, r/(r+mu), times)
+            
+            for j in range(times):
+                sampled_vals[i, j] = vals[j]
+                sampled_pvals[i, j] = nbinom.cdf(<int>vals[j], r/(r+mu), r)
+
+        return sampled_vals, sampled_pvals
 
 def learn_dispersion_model(h, cutoff = 250, trim = (2.5, 97.5)):
-	"""Learn a dispersion model from the expected vs. observed histogram
-	
-	:param h: 2-D :class:`numpy.array` containing the distribution of observerd cleavages at each expected cleavage rate
-	:type h: numpy.array
-	:param cutoff: Mininum number of observed cleavages tp perform ML negative binomial fit
-	:type cutoff: int
-	:param trim: Percent of data to trim  from the observed cleavage count (to mitigate outlier effects)
-	:type trim: tuple of int
-
-	:return: A dispersion model class :class:`dispersion_model`
-	:rtype: dispersion_model
-
-	:todo: Add exceptions for failure to fit, etc.
-	"""
-
-	size = int(h.shape[0])
-	p = np.zeros(size)
-	r = np.zeros(size)
-
-	# Make an initial negative binomial fit
-	for i in range(0, size):
-	
-		pos = 0
-		x = np.zeros(int(np.sum(h[i,:])))
-		for j in range(len(h[i,:])):
-			num = int(h[i, j])
-			x[pos:pos+num] = j
-			pos += num
-
-		# If more than 500k points downsample to
-		# make curve-fitting tractable
-		if len(x)>1e5:
-			x=np.random.choice(x, size=int(1e5))
-			x=np.sort(x)
-
-		if len(x) >= cutoff:		
-
-			lower = int(np.floor(x.shape[0]*(trim[0]/100.0)))
-			upper = int(np.ceil(x.shape[0]*(trim[1]/100.0)))
-
-			mu = np.mean(x[lower:upper])
-			var = np.var(x[lower:upper])
-
-			est_r = (mu * mu) / (var - mu)
-			if est_r <= 0.0: 
-				est_r = 10.0
-			est_p = est_r / (est_r + mu)
-			
-			(p[i], r[i]) = nbinom.fit(x[lower:upper], p = est_p, r = est_r)
-
-		else:
-			p[i] = r[i] = np.nan
-
-	# Back-compute the mean values from the negative binomial parameters
-	mus = p*r/(1-p)
-
-	# hack
-	r[r>200]=200.
-
-	# Extrapolate using polynomial fit
-	x = np.arange(size)
-	sele = np.isfinite(mus)
-
-	first_x = np.min(x[sele])
-	last_x = np.max(x[sele])*0.75
-
-	# fit mu with a 3 segments
-	fit_mu=pwlf.PiecewiseLinFit(x[sele], mus[sele])
-	res=fit_mu.fit_with_breaks_force_points(np.linspace(first_x, last_x, 4), [0], [mus[0]])
-
-	# fit r with 5 s segments
-	fit_r=pwlf.PiecewiseLinFit(x[sele], 1.0/r[sele])
-	res = optimize.minimize(fit_r.fit_with_breaks_opt, [3.0, 7.0, 15.0, 25.0])
-	
-	x0=np.zeros(6)
-	x0[0]=first_x
-	x0[-1]=last_x
-	x0[1:-1]=res.x
-
-	#res=fit_r.fit_with_breaks_force_points([first_x, 2, 7, 15, 25, last_x], [1], [1.0/r[1]])
-	res=fit_r.fit_with_breaks_force_points(x0, [1], [1.0/r[1]])
-
-
-	# Create a dispersion model class
-	res = dispersion_model()
-	res.h = h
-	res.p = p
-	res.r = r
-
-	res.mu_params = list(fit_mu.fit_breaks[1:]) + list(fit_mu.intercepts) + list(fit_mu.slopes)
-	res.r_params = list(fit_r.fit_breaks[1:]) + list(fit_r.intercepts) + list(fit_r.slopes)
-	
-	return res
-
-# Read and write functions for making a portable dispersion model
-
-import simplejson as json
+    """Learn a dispersion model from the expected 
+    vs. observed histogram
+
+    Parameters
+    ----------
+    h : ndarray
+        A 2-dimemsional array containing the distribution of
+        observerd cleavages at each expected cleavage rate
+    cutoff : int
+        Mininum number of observed cleavages to perform ML 
+        negative binomial fit at each value of expected cleavages
+    trim : tuple (float)
+        Percent of data to trim  from the observed cleavage count 
+        (to mitigate outlier effects)
+
+    Returns
+    -------
+    model : dispersion_model
+        A dispersion model learned from observed and expected counts
+
+    Todo
+    ----
+    Add exceptions for failure to fit, etc.
+    """
+    size = int(h.shape[0])
+    p = np.zeros(size)
+    r = np.zeros(size)
+
+    logger.info("Analyzing {:,} cleavages".format(np.sum(h)))
+
+    # Make an initial negative binomial fit
+    for i in range(0, size):
+    
+        # Unpack histogram
+        pos = 0
+        x = np.zeros(int(np.sum(h[i,:])))
+        for j in range(len(h[i,:])):
+            num = int(h[i, j])
+            x[pos:pos+num] = j
+            pos += num
+
+        # If more than 500k points downsample to
+        # make curve-fitting tractable
+        if len(x)>1e5:
+            x=np.random.choice(x, size=int(1e5))
+            x=np.sort(x)
+
+        if len(x) >= cutoff:		
+            # Find data points to trim
+            lower = int(np.floor(x.shape[0]*(trim[0]/100.0)))
+            upper = int(np.ceil(x.shape[0]*(trim[1]/100.0)))
+
+            # Compute intial estimates of mean and variance
+            mu = np.mean(x[lower:upper])
+            var = np.var(x[lower:upper])
+
+            est_r = (mu * mu) / (var - mu)
+            if est_r <= 0.0: 
+                est_r = 10.0
+            est_p = est_r / (est_r + mu)
+            
+            # Fit negative binomial
+            (p[i], r[i]) = nbinom.fit(x[lower:upper], p = est_p, r = est_r)
+
+        else:
+            # Not enough data points to fit - will exclude 
+            # from further analysis
+            p[i] = r[i] = np.nan
+
+    # Back-compute the mean values from the negative binomial parameters
+    mus = p*r/(1-p)
+
+    # Hack to control for problematic fits
+    r[r>200] = 200.
+
+    # Extrapolate using polynomial fit
+    x = np.arange(size)
+    sele = np.isfinite(mus)
+
+    first_x = np.min(x[sele])
+    last_x = np.max(x[sele])*0.75
+
+    # fit mu with a 3 segments
+    fit_mu = pwlf.PiecewiseLinFit(x[sele], mus[sele])
+    res = fit_mu.fit_with_breaks_force_points(np.linspace(first_x, last_x, 4), [0], [mus[0]])
+
+    # fit r with 5 s segments
+    fit_r = pwlf.PiecewiseLinFit(x[sele], 1.0/r[sele])
+    res = optimize.minimize(fit_r.fit_with_breaks_opt, [3.0, 7.0, 15.0, 25.0])
+    
+    x0 = np.zeros(6)
+    x0[0] = first_x
+    x0[-1] = last_x
+    x0[1:-1] = res.x
+
+    res = fit_r.fit_with_breaks_force_points(x0, [1], [1.0/r[1]])
+    
+    # Create a dispersion model class
+    res = dispersion_model()
+    res.h = h
+    res.p = p
+    res.r = r
+
+    res.mu_params = list(fit_mu.fit_breaks[1:]) + list(fit_mu.intercepts) + list(fit_mu.slopes)
+    res.r_params = list(fit_r.fit_breaks[1:]) + list(fit_r.intercepts) + list(fit_r.slopes)
+    
+    return res
+
 import base64
 
 def base64encode(x):
-	
-	return [str(x.dtype), base64.b64encode(x), x.shape]
+    return [str(x.dtype), base64.b64encode(x), x.shape]
 
 def base64decode(x):
+    dtype = np.dtype(x[0])
+    arr = np.frombuffer(base64.b64decode(x[1]), dtype)
+    if len(x) > 2:
+        return arr.reshape(x[2])
+    return arr
 
-	dtype = np.dtype(x[0])
-	arr = np.frombuffer(base64.b64decode(x[1]), dtype)
-	if len(x) > 2:
-		return arr.reshape(x[2])
-	return arr
+def load_dispersion_model(filename):
+    """Load a dispersion model encoded in JSON format
 
-def read_dispersion_model(filename):
-	
-	import urllib.request as request
+    Parameters
+    ----------
+    filename : str
+        Path to JSON-format dispersion model
+    """
+    import simplejson as json
+    import urllib.request as request
 
-	if filename.startswith('http'):
-		file = request.urlopen(filename)
-	else:
-		file = open(filename, 'r')
+    if filename.startswith('http'):
+        file = request.urlopen(filename)
+    else:
+        file = open(filename, 'r')
 
-	params = json.load(file)
+    params = json.load(file)
 
-	file.close()
+    file.close()
 
-	model = dispersion_model()
-	model.mu_params = base64decode(params["mu_params"])
-	model.r_params = base64decode(params["r_params"])
+    model = dispersion_model()
+    model.mu_params = base64decode(params["mu_params"])
+    model.r_params = base64decode(params["r_params"])
 
-	if "h" in params:
-		model.h = base64decode(params["h"])
-	if "p" in params:
-		model.p = base64decode(params["p"])
-	if "r" in params:
-		model.r = base64decode(params["r"])
+    if "h" in params:
+        model.h = base64decode(params["h"])
+    if "p" in params:
+        model.p = base64decode(params["p"])
+    if "r" in params:
+        model.r = base64decode(params["r"])
 
-	return model
+    return model
 
 def write_dispersion_model(model):
-	"""Write a JSON format dispersion model
-
-	:param model: An instance of :class:`dispersion_model`
-	:return: JSON-formatted dump of dispersion model
-	:rtype: str
-	"""
-
-	out = { "mu_params": base64encode(np.asarray(model.mu_params, order = 'C')), 
-			"r_params": base64encode(np.asarray(model.r_params, order = 'C')),
-			"h": base64encode(np.asarray(model.h, order = 'C')),
-			"p": base64encode(np.asarray(model.p, order = 'C')),
-			"r": base64encode(np.asarray(model.r, order = 'C'))
-		}
-
-	return json.dumps(out, indent = 4)
+    """Write a JSON format dispersion model
+
+    Parameters
+    ----------
+    model : dispersion_model
+        An instance of dispersion_model
+
+    Returns
+    -------
+    out : str
+        JSON-formatted dump of dispersion model
+    """
+    import simplejson as json
+
+    out = { "mu_params": base64encode(np.asarray(model.mu_params, order = 'C')), 
+            "r_params": base64encode(np.asarray(model.r_params, order = 'C')),
+            "h": base64encode(np.asarray(model.h, order = 'C')),
+            "p": base64encode(np.asarray(model.p, order = 'C')),
+            "r": base64encode(np.asarray(model.r, order = 'C')),
+    }
+
+    return json.dumps(out, indent = 4)
 
diff --git a/footprint_tools/modeling/predict.h b/footprint_tools/modeling/predict.h
index 6e99dff..41a9a02 100644
--- a/footprint_tools/modeling/predict.h
+++ b/footprint_tools/modeling/predict.h
@@ -1,4 +1,4 @@
-// Copyright 2015 Jeff Vierstra
+// Fast predict functions implemented in native C
 
 #ifndef __PREDICT_H_
 #define __PREDICT_H_
@@ -8,62 +8,68 @@
 
 typedef struct result
 {
-	double *exp;
-	double *win;
+    double *exp;
+    double *win;
 } result_t;
 
 void free_result_t(result_t* res)
 {
-	free(res->exp);
-	free(res->win);
-	free(res);
+    free(res->exp);
+    free(res->win);
+    free(res);
 }
 
-result_t* fast_predict(double const *obs, double const *probs, int l, int half_window_width, int smoothing_half_window_width, double smoothing_clip)
+result_t* fast_predict(
+        double const *obs, 
+        double const *probs, 
+        int l, 
+        int half_win_width, 
+        int smoothing_half_win_width, 
+        double smoothing_clip)
 {
-	//probs = l
-	int i, j;
-	//int w = (half_window_width * 2);
+    //probs = l
+    int i, j;
+    //int w = (half_window_width * 2);
 
-	// Create memory space for windows
-	double *win_counts	= (double*) calloc(l, sizeof(double));
-	double *win_probs 	= (double*) calloc(l, sizeof(double));
-	double *exp_counts	= (double*) calloc(l, sizeof(double));
+    // Create memory space for windows
+    double *win_counts	= (double*) calloc(l, sizeof(double));
+    double *win_probs 	= (double*) calloc(l, sizeof(double));
+    double *exp_counts	= (double*) calloc(l, sizeof(double));
 
-	// Make windowed counts and probabilities
-	for (i = half_window_width; i < l-half_window_width; i++)
-	{
-		for (j = -half_window_width; j < half_window_width; j++)
-		{
-			win_counts[i] += obs[i+j];
-			win_probs[i] += probs[i+j]; 
-		}
-	}
+    // Make windowed counts and probabilities
+    for (i = half_win_width; i < l-half_win_width; i++)
+    {
+        for (j = -half_win_width; j < half_win_width; j++)
+        {
+            win_counts[i] += obs[i+j];
+            win_probs[i] += probs[i+j]; 
+        }
+    }
 
-	if (smoothing_half_window_width > 0)
-	{
-		double *smoothed_win_counts = windowed_trimmed_mean(&win_counts[0], l, smoothing_half_window_width, smoothing_clip);
-		// free
-		free(win_counts);
-		// set
-		win_counts = smoothed_win_counts;
-	}
+    if (smoothing_half_win_width > 0)
+    {
+        double *smoothed_win_counts = windowed_trimmed_mean(&win_counts[0], l, smoothing_half_win_width, smoothing_clip);
+        // free
+        free(win_counts);
+        // set
+        win_counts = smoothed_win_counts;
+    }
 
-	// Compute expected
-	for (i = half_window_width; i < l-half_window_width; i++)
-	{
-		exp_counts[i] = c_round( (probs[i] / win_probs[i]) * win_counts[i] );
-	}
+    // Compute expected
+    for (i = half_win_width; i < l-half_win_width; i++)
+    {
+        exp_counts[i] = c_round( (probs[i] / win_probs[i]) * win_counts[i] );
+    }
 
-	result_t *res = (result_t*) malloc(sizeof(result_t));
-	res->exp = exp_counts;
-	res->win = win_counts;
+    result_t *res = (result_t*) malloc(sizeof(result_t));
+    res->exp = exp_counts;
+    res->win = win_counts;
 
-	// Clean up
-	free(win_probs);
+    // Clean up
+    free(win_probs);
 
-	// Return
-	return res;
+    // Return
+    return res;
 }
 
 #endif
\ No newline at end of file
diff --git a/footprint_tools/modeling/predict.pyx b/footprint_tools/modeling/predict.pyx
index ee2569b..58102ca 100644
--- a/footprint_tools/modeling/predict.pyx
+++ b/footprint_tools/modeling/predict.pyx
@@ -1,177 +1,161 @@
-# Copyright 2015 Jeff Vierstra
+"""
 
+"""
 # cython: boundscheck=False
 # cython: wraparound=False
 # cython: nonecheck=False
 # cython: embedsignature=True
 
 cimport numpy as np
-
 import numpy as np
 
-from genome_tools import genomic_interval, genomic_interval_set
+from genome_tools import genomic_interval
 
 cdef extern from "predict.h":
-	struct result:
-		double* exp
-		double* win	
-	void free_result_t(result*)
-	result* fast_predict(double*, double*, int, int, int, double)
+    struct result:
+        double* exp
+        double* win	
+    void free_result_t(result*)
+    result* fast_predict(double*, double*, int, int, int, double)
 
 ctypedef np.float64_t data_type_t
 
-cdef predict(data_type_t [:] obs, data_type_t [:] probs, int half_window_width, int smoothing_half_window_width, double smoothing_clip):
+cdef predict(data_type_t [:] obs, data_type_t [:] probs, int half_win_width, int smoothing_half_win_width, double smoothing_clip):
 
-	cdef int i
-	cdef int l = obs.shape[0]
+    cdef int i
+    cdef int l = obs.shape[0]
 
-	# Predict
-	cdef result* res = fast_predict(&obs[0], &probs[0], l, half_window_width, smoothing_half_window_width, smoothing_clip)
+    # Predict
+    cdef result* res = fast_predict(&obs[0], &probs[0], l, half_win_width, smoothing_half_win_width, smoothing_clip)
 
-	# Copy back into a python object
-	cdef np.ndarray[data_type_t, ndim = 1, mode = 'c'] exp = np.zeros(l, dtype = np.float64, order = 'c')
-	cdef data_type_t [:] exp_view = exp
+    # Copy back into a python object
+    cdef np.ndarray[data_type_t, ndim = 1, mode = 'c'] exp = np.zeros(l, dtype = np.float64, order = 'c')
+    cdef data_type_t [:] exp_view = exp
 
-	cdef np.ndarray[data_type_t, ndim = 1, mode = 'c'] win = np.zeros(l, dtype = np.float64, order = 'c')
-	cdef data_type_t [:] win_view = win
+    cdef np.ndarray[data_type_t, ndim = 1, mode = 'c'] win = np.zeros(l, dtype = np.float64, order = 'c')
+    cdef data_type_t [:] win_view = win
 
-	for i in range(l): 
-		exp_view[i] = res.exp[i]
-		win_view[i] = res.win[i]
+    for i in range(l): 
+        exp_view[i] = res.exp[i]
+        win_view[i] = res.win[i]
 
-	# Free up the memory from C
-	free_result_t(res)
+    # Free up the memory from C
+    free_result_t(res)
 
-	return exp, win
+    return exp, win
 
 def reverse_complement(seq):
-	"""Computes reverse complement of a DNA sequence
-	
-	Parameters
-	----------
-	seq : str
-	    DNA sequence string
- 	
-	Returns
-	-------
-	str
-	    Reverse complement of ``seq``
-	"""
-	compl = { 'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N': 'N', 'a': 't', 'c': 'g', 'g': 'c', 't': 'a', 'n': 'n'}
-	return ''.join([ compl.get(base, 'N') for base in seq ])[::-1]
+    """Computes reverse complement of a DNA sequence
+    
+    Parameters
+    ----------
+    seq : str
+        DNA sequence string
+     
+    Returns
+    -------
+    out : str
+        Reverse complement of ``seq``
+    """
+    compl = { 'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N': 'N', 'a': 't', 'c': 'g', 'g': 'c', 't': 'a', 'n': 'n'}
+    return ''.join([ compl.get(base, 'N') for base in seq ])[::-1]
 
 class prediction(object):
-
-	"""Class that holds a wrapper function to 
-		compute the expected cleavage counts
-	
-	Attributes
-	----------
-	bm : bias.bias_model
-	    Instance of bias model class
-	fasta_func : pysam.FastaFile
-		Instance of FastaFile
-	half_window_width : int
-	    Window width to apply bias model (final windows size = 2W+1)
-	padding : int
-	    Padding applied to ``interval`` when retrieving per-nucleotide data
-	read_func : cutcounts.bamfile
-		Instance of counts reader
-	smoothing_clip : float
-	    Fraction of nucleotides to trim when computing smoothed mean
-	smoothing_half_window_width : int
-	    Desc
-	"""
-	
-	def __init__(self, read_func, fasta_func, bm, half_window_width = 5, smoothing_half_window_width = 0, smoothing_clip = 0.01):
-		"""
-		
-		Parameters
-		----------
-		read_func : TYPE
-		    Description
-		fasta_func : TYPE
-		    Description
-		bm : TYPE
-		    Description
-		half_window_width : int, optional
-		    Description
-		smoothing_half_window_width : int, optional
-		    Description
-		smoothing_clip : float, optional
-		    Description
-		"""
-		self.read_func = read_func
-		self.fasta_func = fasta_func
-		self.bm = bm
-
-		# Smoothing/windowing parameters
-		self.half_window_width = half_window_width
-		self.smoothing_half_window_width = smoothing_half_window_width
-		self.smoothing_clip = smoothing_clip
-
-		# Base padding
-		self.padding = self.half_window_width + smoothing_half_window_width
-
-		# self.interval = interval
-
-		# pad_interval = interval.widen(self.padding)
-		
-		# Note: We clip the first base when recombining the positive 
-		# and negative strand, so add an extra base upfront
-		# pad_interval.start -= 1
-
-		# self.counts = reads[pad_interval]
-		# self.seq = fasta[pad_interval.chrom][pad_interval.start-bm.offset():pad_interval.end+bm.offset()].seq.upper()
-
-	def compute(self, x):
-		"""Summary
-		
-		Parameters
-		----------
-		x : genome_toools.genomic_interval,  genome_toools.genomic_interval_set, or list
-		    Genomic region(s) to generate predicted cleavages
-		
-		Returns
-		-------
-		tuple (or list of tuples)
-		    Observed, expected and windowed cleavage counts
-		"""
-
-		# if isinstance(x, (list, genomic_interval_set)):
-		# 	print("LIST")
-		# 	for i in x:
-		# 		yield self.compute(i)
-	
-		# Note: We clip the first base when recombining the positive 
-		# and negative strand, so add an extra base upfront
-		pad_interval = x.widen(self.padding)
-		pad_interval.start -= 1
-
-		# Get the raw cleavage counts and FASTA sequence
-		raw_counts = self.read_func[pad_interval]
-		raw_seq = self.fasta_func.fetch(pad_interval.chrom, pad_interval.start-self.bm.offset(), pad_interval.end+self.bm.offset()).upper()
-
-		obs_counts = {'+': None, '-': None}
-		exp_counts = {'+': None, '-': None}
-		win_counts = {'+': None, '-': None}
-
-		for strand in ['+', '-']:
-
-			# Pre-calculate the sequence bias propensity table from bias model
-			if strand == '+':
-				probs = self.bm.probs(raw_seq)
-			else:
-				probs = self.bm.probs(reverse_complement(raw_seq))[::-1]
-
-			exp, win = predict(np.ascontiguousarray(raw_counts[strand]), np.ascontiguousarray(probs), self.half_window_width, self.smoothing_half_window_width, self.smoothing_clip)
-
-			w = raw_counts[strand].shape[0] - self.padding
-
-			obs_counts[strand] = raw_counts[strand][self.padding:w]
-			exp_counts[strand] = exp[self.padding:w]
-			win_counts[strand] = win[self.padding:w]
-
-		return obs_counts, exp_counts, win_counts
-
-
+    """Class that holds a wrapper function to 
+        compute the expected cleavage counts
+    
+    Attributes
+    ----------
+    bm : bias.bias_model
+        Sequence bias model to apply
+    read_func : cutcounts.bamfile
+       Cut-counts reader
+    fasta_func : pysam.FastaFile
+        FASTA-file reader
+    half_win_width : int
+        Window width to apply bias model (final windows size = 2W+1)
+    padding : int
+        Padding applied to region when retrieving per-nucleotide data
+    smoothing_clip : float
+        Fraction of nucleotides to trim when computing smoothed mean
+    smoothing_half_win_width : int
+        Half width of window used to compute windowed tag counts 
+    """
+    
+    def __init__(self, read_func, fasta_func, bm, half_win_width=5, smoothing_half_win_width=0, smoothing_clip=0.01):
+        """
+        
+        Parameters
+        ----------
+        read_func : cutcounts.bamfile
+            Description
+        fasta_func : pysam.FastaFile
+            Description
+        bm : bias.bias_model
+            Description
+        half_win_width : int, optional
+            Description
+        smoothing_half_win_width : int, optional
+            Description
+        smoothing_clip : float, optional
+            Half width of window used to compute windowed tag counts 
+        """
+
+        self.read_func = read_func
+        self.fasta_func = fasta_func
+        self.bm = bm
+
+        # Smoothing/windowing parameters
+        self.half_win_width = half_win_width
+        self.smoothing_half_win_width = smoothing_half_win_width
+        self.smoothing_clip = smoothing_clip
+
+        # Base padding
+        self.padding = self.half_win_width + smoothing_half_win_width
+
+    def compute(self, x):
+        """Computed expected cleavage counts
+        
+        Parameters
+        ----------
+        x : genome_toools.genomic_interval
+            Genomic region to generate predicted cleavages
+        
+        Returns
+        -------
+        out: tuple of dict
+            Observed, expected and windowed cleavage counts
+        """
+    
+        # Note: We clip the first base when recombining the positive 
+        # and negative strand, so add an extra base upfront
+        pad_interval = x.widen(self.padding)
+        pad_interval.start -= 1
+
+        # Get the raw cleavage counts and FASTA sequence
+        raw_counts = self.read_func[pad_interval]
+        raw_seq = self.fasta_func.fetch(pad_interval.chrom, 
+                                        pad_interval.start-self.bm.offset(), 
+                                        pad_interval.end+self.bm.offset()).upper()
+
+        obs_counts = {'+': None, '-': None}
+        exp_counts = {'+': None, '-': None}
+        win_counts = {'+': None, '-': None}
+
+        for strand in ['+', '-']:
+
+            # Pre-calculate the sequence bias propensity table from bias model
+            if strand == '+':
+                probs = self.bm.probs(raw_seq)
+            else:
+                probs = self.bm.probs(reverse_complement(raw_seq))[::-1]
+
+            exp, win = predict(np.ascontiguousarray(raw_counts[strand]), np.ascontiguousarray(probs), self.half_win_width, self.smoothing_half_win_width, self.smoothing_clip)
+
+            w = raw_counts[strand].shape[0] - self.padding
+
+            obs_counts[strand] = raw_counts[strand][self.padding:w]
+            exp_counts[strand] = exp[self.padding:w]
+            win_counts[strand] = win[self.padding:w]
+
+        return obs_counts, exp_counts, win_counts
diff --git a/footprint_tools/modeling/smoothing.h b/footprint_tools/modeling/smoothing.h
index 8ad1bbf..b3c5fb6 100644
--- a/footprint_tools/modeling/smoothing.h
+++ b/footprint_tools/modeling/smoothing.h
@@ -1,4 +1,4 @@
-// Copyright 2015 Jeff Vierstra
+// Windowed smoothing functions implemented in native C
 
 #ifndef __SMOOTHING_H_
 #define __SMOOTHING_H_
@@ -104,20 +104,20 @@ double trimmed_mean(double x[], int n, int k)
 }
 
 
-double* windowed_trimmed_mean(double* win, int l, int half_window_width, double clip)
+double* windowed_trimmed_mean(double* win, int l, int half_win_width, double clip)
 {
     int i, j, k;
-    int w = (half_window_width * 2) + 1;
+    int w = (half_win_width * 2) + 1;
 
     k = (int)( (double)w * clip );
 
     double *tmp = malloc(w * sizeof(double));
     double *res = calloc(l, sizeof(double));
     
-    for (i = half_window_width; i < l-half_window_width; i++)
+    for (i = half_win_width; i < l-half_win_width; i++)
     {
         // index to start
-        j = i - half_window_width;
+        j = i - half_win_width;
 
         // copy data into new memory location
         memcpy(&tmp[0], &win[j], w * sizeof(double));
diff --git a/footprint_tools/stats/__init__.py b/footprint_tools/stats/__init__.py
index f1f1f5d..649f355 100644
--- a/footprint_tools/stats/__init__.py
+++ b/footprint_tools/stats/__init__.py
@@ -1,1 +1,5 @@
-__all__ = ["distributions", "windowing", "bayesian", "segment"]
+"""
+This module contains classes and functions for statistical analysis
+"""
+
+__all__ = ["distributions", "windowing", "posterior", "utils"]
diff --git a/footprint_tools/stats/bayesian.py b/footprint_tools/stats/bayesian.py
deleted file mode 100644
index ffca20c..0000000
--- a/footprint_tools/stats/bayesian.py
+++ /dev/null
@@ -1,142 +0,0 @@
-# Copyright 2015 Jeff Vierstra
-
-import numpy as np
-import scipy.stats
-
-from footprint_tools.stats import windowing
-
-def compute_prior_weighted(fdr, w, cutoff = 0.05, pseudo = 0.5):
-	"""Returns prior of whether a TF is bound or not 
-	
-	Parameters
-	----------
-	fdr : numpy.array
-	    2D :class:`np.array` with FPR values for each datset and position
-	w : numpy.array
-	    Binary array whether samples is DHS
-	cutoff : float, optional
-	    FPR cutoff value to label as occupied when building prior
-	pseudo : float, optional
-	    Psuedocount to add as prior to Beta distribution
-	
-	Returns
-	-------
-	TYPE
-	    Description
-	"""
-
-	k = np.sum(fdr <= cutoff, axis = 0) # num of datasets w/ fp
-	n = np.sum(w, axis = 0) #np.sum(w) # num of datasets w/ HS
-	A = n-k + pseudo
-	B = k + pseudo
-	pr = A/(A+B)
-	
-	res = np.ones(fdr.shape)
-	res *= pr[np.newaxis,:]
-	res[w==0] = 1
-
-	return res
-
-def compute_delta_prior(obs, exp, fdr, beta_prior, cutoff = 0.05):
-	"""Returns a point estimate of exepected cleavage 
-	depletion with a footprint at each nucleotide 
-	
-	Parameters
-	----------
-	obs : np.array
-	     2D :class:`np.array` with observed cleavages counts
-	exp : np.array
-	     2D :class:`np.array` with expectd cleavages counts
-	fdr : np.array
-	    2D :class:`np.array` with adjusted p-values (FDR)
-	beta_prior : np.array
-	    Description
-	cutoff : float, optional
-	    FPR cutoff value to label as occupied when building prior
-	
-	Returns
-	-------
-	TYPE
-	    Description
-	"""
-
-	(n, w) = obs.shape
-
-	mus = np.ones((n, w))
-	ws = np.ones((n, w))
-	
-	for i in range(n):
-		
-		k = obs[i,:]
-		n = np.max(np.vstack([exp[i,:], obs[i,:]]), axis = 0)
-
-		mu, v = scipy.stats.beta.stats(k+beta_prior[i][0], n-k+beta_prior[i][1], 
-			loc = 0, scale = 1, moments = 'mv')
-
-		mus[i,:] = mu
-		ws[i,:] = 1/np.sqrt(v)
-
-	ws[fdr > cutoff] = 0
-
-	delta = np.sum(ws*mus, axis = 0) / np.sum(ws, axis = 0)
-	delta[np.isnan(delta)] = 1
-
-	return delta
-
-def log_likelihood(obs, exp, dm, delta = 1, w = 3):
-	"""
-	
-	Parameters
-	----------
-	obs : TYPE
-	    Description
-	exp : TYPE
-	    Description
-	dm : TYPE
-	    Description
-	delta : int, optional
-	    Description
-	w : int, optional
-	    Description
-	
-	Returns
-	-------
-	TYPE
-	    Description
-	"""
-	res = np.ones((obs.shape[0], obs.shape[1]), order = 'c')
-
-	n = obs.shape[0]
-	for i in range(n):
-		res[i,:] = windowing.sum(dm[i].log_pmf_values(exp[i,:] * delta, obs[i,:]), w)
-	
-	return res
-
-def posterior(prior, ll_on, ll_off):
-	"""
-	Parameters
-	----------
-	prior : TYPE
-	    Description
-	ll_on : TYPE
-	    Description
-	ll_off : TYPE
-	    Description
-	
-	Returns
-	-------
-	TYPE
-	    Description
-	"""
-	prior_on = np.log(1-prior)
-	prior_off = np.log(prior)
-
-	p_off = prior_off + ll_off
-	p_on = prior_on + ll_on
-	denom = np.logaddexp(p_on, p_off)
-
-	return (p_off - denom)
-
-
-
-
diff --git a/footprint_tools/stats/differential.pyx b/footprint_tools/stats/differential.pyx
index 382f142..36f66e9 100644
--- a/footprint_tools/stats/differential.pyx
+++ b/footprint_tools/stats/differential.pyx
@@ -1,12 +1,11 @@
 
-# Copyright 2017 Jeff Vierstra
+# Copyright 2017-2021 Jeff Vierstra
 
 # cython: boundscheck=False
 # cython: wraparound=False
 # cython: nonecheck=False
 # cython: embedsignature=True
 
-
 import numpy as np
 cimport numpy as np
 
@@ -16,79 +15,85 @@ import scipy.optimize
 ctypedef np.float64_t data_type_t
 
 from ..modeling cimport dispersion
+from .distributions cimport invchi2, t
 
 ctypedef dispersion.dispersion_model dm_t
 
 cdef extern from "cephes.h":
-	double c_lgamma(double) nogil
-	double c_log(double) nogil
-	double c_sqrt(double) nogil
-	double PI
+    double c_lgamma(double) nogil
+    double c_log(double) nogil
+    double c_sqrt(double) nogil
+    double PI
 
 cpdef np.ndarray[data_type_t, ndim = 3, mode = 'c'] compute_logpmf_values(dm, np.ndarray[data_type_t, ndim = 2, mode = 'c'] obs, np.ndarray[data_type_t, ndim = 2, mode = 'c'] exp,  data_type_t lo, data_type_t hi, int nslices):
-	
-	cdef int i 		
-	cdef int n = obs.shape[0]
-	cdef int m = obs.shape[1]
-
-	bins = np.power(2, np.linspace(lo, hi, nslices))
+    """Compute log pmf
+    """
+    cdef int i 		
+    cdef int n = obs.shape[0]
+    cdef int m = obs.shape[1]
 
-	cdef np.ndarray[data_type_t, ndim = 3, mode = 'c'] e = np.ascontiguousarray(exp[:,:,np.newaxis] * bins)
-	cdef np.ndarray[data_type_t, ndim = 3, mode = 'c'] o = np.ascontiguousarray(np.repeat(obs[:,:,np.newaxis], nslices, axis = 2))
-	cdef np.ndarray[data_type_t, ndim = 3, mode = 'c'] res = np.zeros((n, m, nslices), order = 'c')
+    bins = np.power(2, np.linspace(lo, hi, nslices))
 
-	for i in range(n):
-		<dm_t>(dm[i]).log_pmf_values_0(np.ravel(e[i,:,:]), np.ravel(o[i,:,:]), np.ravel(res[i,:,:]))
+    cdef np.ndarray[data_type_t, ndim = 3, mode = 'c'] e = np.ascontiguousarray(exp[:,:,np.newaxis] * bins)
+    cdef np.ndarray[data_type_t, ndim = 3, mode = 'c'] o = np.ascontiguousarray(np.repeat(obs[:,:,np.newaxis], nslices, axis = 2))
+    cdef np.ndarray[data_type_t, ndim = 3, mode = 'c'] res = np.zeros((n, m, nslices), order = 'c')
 
-	return res.T
+    for i in range(n):
+        <dm_t>(dm[i]).log_pmf_values_0(np.ravel(e[i,:,:]), np.ravel(o[i,:,:]), np.ravel(res[i,:,:]))
 
+    return res.T
 
 #various stats functions (should be move to ..stats/distributions folder)
 
-def invchi2_logpdf(x, nu, tau2):
-	return (nu/2.0)*np.log(tau2*nu/2.0) + (-nu*tau2/2.0/x) - scipy.special.gammaln(nu/2.0) - (1.0+nu/2.0)*np.log(x)
+"""def invchi2_logpdf(x, nu, tau2):
+    return (nu/2.0)*np.log(tau2*nu/2.0) + (-nu*tau2/2.0/x) - scipy.special.gammaln(nu/2.0) - (1.0+nu/2.0)*np.log(x)
 
 def invchi2_likelihood(par, x):
-	return -np.sum(invchi2_logpdf(x, par[0], par[1]))
+    return -np.sum(invchi2_logpdf(x, par[0], par[1]))
 
 cdef double t_logpmf(double x, double nu, double mu, double sig2) nogil:
-	return c_lgamma((nu+1.0)/2.0) - c_lgamma(nu/2.0) - c_log(c_sqrt(PI*nu*sig2)) - ((nu+1.0)/2.0 * c_log(1.0+((x-mu)**2)/sig2/nu ))
-
+    return c_lgamma((nu+1.0)/2.0) - c_lgamma(nu/2.0) - c_log(c_sqrt(PI*nu*sig2)) - ((nu+1.0)/2.0 * c_log(1.0+((x-mu)**2)/sig2/nu ))
+"""
 
-cpdef np.ndarray[data_type_t, ndim = 2, mode = 'c'] compute_log_prior_t(np.ndarray[data_type_t, ndim = 2, mode = 'c'] obs_over_exp , data_type_t nu_0, data_type_t sig2_0,  data_type_t lo, data_type_t hi, int nslices):
-	
-	""" 
-	Computes the log posterior distribution of the hyperparameters (normal * inv chi2 = students t)
-	"""
+cpdef np.ndarray[data_type_t, ndim = 2, mode = 'c'] compute_log_prior_t(np.ndarray[data_type_t, ndim = 2, mode = 'c'] obs_over_exp , data_type_t nu_0, data_type_t sig2_0,  data_type_t lo, data_type_t hi, int nslices): 
+    """ Computes the log posterior distribution of the hyperparameters 
+    (normal * inv chi2 = students t)
+    
+    Parameters
+    ----------
 
-	cdef int i, j
-	cdef int n = obs_over_exp.shape[0]
-	cdef int m = obs_over_exp.shape[1]
+    Returns
+    -------
+    
+    """
+    cdef int i, j
+    cdef int n = obs_over_exp.shape[0]
+    cdef int m = obs_over_exp.shape[1]
 
-	cdef data_type_t mu, ssqdev, nu_1, sig2_1
+    cdef data_type_t mu, ssqdev, nu_1, sig2_1
 
-	cdef np.ndarray[data_type_t, ndim = 2, mode = 'c'] res = np.zeros((m, nslices), dtype = np.float64, order = 'c')
-	cdef data_type_t [:,:] res_view = res
+    cdef np.ndarray[data_type_t, ndim = 2, mode = 'c'] res = np.zeros((m, nslices), dtype = np.float64, order = 'c')
+    cdef data_type_t [:,:] res_view = res
 
-	cdef double h = (hi-lo)/<double>nslices
+    cdef double h = (hi-lo)/<double>nslices
 
-	for i in range(m):
-		
-		mu = np.nanmean(obs_over_exp[:,i])
-		ssqdev = np.nansum((obs_over_exp[:,i]-mu)**2)
+    for i in range(m):
+        
+        mu = np.nanmean(obs_over_exp[:,i])
+        ssqdev = np.nansum((obs_over_exp[:,i]-mu)**2)
 
-		nu_1 = nu_0 + n
-		sig2_1 = (nu_0*sig2_0 + ssqdev) / (nu_0 + n)
+        nu_1 = nu_0 + n
+        sig2_1 = (nu_0*sig2_0 + ssqdev) / (nu_0 + n)
 
-		with nogil:
+        with nogil:
 
-			for j in range(0, nslices):
-				#res_view[i, j] = scipy.stats.t.logpdf(bins[j], df = nu_1, loc = mu, scale = np.sqrt(sig2_1))
-				res_view[i, j] = t_logpmf(lo+(j*h), nu_1, mu, sig2_1)
+            for j in range(0, nslices):
+                #res_view[i, j] = scipy.stats.t.logpdf(bins[j], df = nu_1, loc = mu, scale = np.sqrt(sig2_1))
+                res_view[i, j] = t.logpmf(lo+(j*h), nu_1, mu, sig2_1)
 
-		res[i, :] += np.log(h)
+        res[i, :] += np.log(h)
 
-	return res.T # - res.sum(axis = 1)[:, np.newaxis]).T
+    return res.T # - res.sum(axis = 1)[:, np.newaxis]).T
 
 
 
diff --git a/footprint_tools/stats/distributions/__init__.py b/footprint_tools/stats/distributions/__init__.py
index 84a7553..4cc6362 100644
--- a/footprint_tools/stats/distributions/__init__.py
+++ b/footprint_tools/stats/distributions/__init__.py
@@ -1,1 +1,5 @@
-__all__ = ["nbinom"]
\ No newline at end of file
+"""
+This module contains imlpementations of statistical distributions.
+"""
+
+__all__ = ["nbinom", "beta", "invchi2"]
\ No newline at end of file
diff --git a/footprint_tools/stats/distributions/invchi2.pxd b/footprint_tools/stats/distributions/invchi2.pxd
new file mode 100644
index 0000000..5743de2
--- /dev/null
+++ b/footprint_tools/stats/distributions/invchi2.pxd
@@ -0,0 +1,7 @@
+cimport cython
+cimport numpy as np
+
+ctypedef np.float64_t data_type_t
+
+cpdef data_type_t logpmf(data_type_t x, data_type_t nu, data_type_t tau2) nogil
+cpdef data_type_t log_likelihood(data_type_t [:] x, data_type_t nu, data_type_t tau2) nogil
diff --git a/footprint_tools/stats/distributions/invchi2.pyx b/footprint_tools/stats/distributions/invchi2.pyx
new file mode 100644
index 0000000..46d4a31
--- /dev/null
+++ b/footprint_tools/stats/distributions/invchi2.pyx
@@ -0,0 +1,51 @@
+"""
+Implementation of inverse Chi-squared distribution
+"""
+
+cimport cython
+cimport numpy as np
+
+cdef extern from "cephes.h":
+    double c_log(double) nogil
+    double c_lgamma(double) nogil
+     
+cpdef data_type_t logpmf(data_type_t x, data_type_t nu, data_type_t tau2) nogil:
+    """Log probability mass function for inverse Chi-squared distribution
+
+    Parameters
+    ----------
+    x : float
+
+    nu : float
+
+    tau2 : float
+
+    Returns
+    -------
+    p : float
+
+    """
+    #return (nu/2.0)*np.log(tau2*nu/2.0) + (-nu*tau2/2.0/x) - scipy.special.gammaln(nu/2.0) - (1.0+nu/2.0)*np.log(x)
+    return nu/2.0 * c_log(tau2*nu/2.0) + (-nu*tau2/2.0/x) - c_lgamma(nu/2.0) - (1.0+nu/2.0)*c_log(x)
+
+cpdef data_type_t log_likelihood(data_type_t [:] x, data_type_t nu, data_type_t tau2) nogil:
+    """Log-likelihood for inverse Chi-squared distribution
+
+    Parameters
+    ----------
+    data : ndarray
+        Data values
+    nu : float
+
+    tau2 : float
+
+    Returns
+    -------
+    log_likelihood : float
+    """
+    cdef int n = x.shape[0]
+    cdef data_type_t res = 0
+
+    for i in range(n):
+        res += logpmf(x[i], nu, tau2)
+    return res
diff --git a/footprint_tools/stats/distributions/nbinom.pyx b/footprint_tools/stats/distributions/nbinom.pyx
index c22a5ef..1b5755d 100644
--- a/footprint_tools/stats/distributions/nbinom.pyx
+++ b/footprint_tools/stats/distributions/nbinom.pyx
@@ -1,8 +1,8 @@
-# Copyright 2015 Jeff Vierstra
-
+"""
+Implementation of negative binomial distribution
+"""
 # cython: embedsignature=True
 
-
 import scipy.special
 import scipy.optimize
 import numpy as np
@@ -11,92 +11,176 @@ cimport cython
 cimport numpy as np
 
 cdef extern from "cephes.h":
-	double c_exp(double) nogil
-	double c_log(double) nogil
-	double c_log1p(double) nogil
-	double c_lgamma(double) nogil
-	double c_incbet(double, double, double) nogil
+    double c_exp(double) nogil
+    double c_log(double) nogil
+    double c_log1p(double) nogil
+    double c_lgamma(double) nogil
+    double c_incbet(double, double, double) nogil
 
 import warnings
 
-
 def mle(par, data, sm):
-
-	"""Objective function for MLE estimate according to
-	https://en.wikipedia.org/wiki/Negative_binomial_distribution#Maximum_likelihood_estimation
-	
-	Parameters
-	----------
-	par : TYPE
-	    Description
-	data : TYPE
-	    Description
-	sm : TYPE
-	    Description
-	data: the points to be fit
-	sm: \sum data / len(data)
-	
-	Returns
-	-------
-	TYPE
-	    Description
-	"""
-
-	p = par[0]
-	r = par[1]
-	n = len(data)
-	f0 = sm / (r+sm) - p
-	f1 = np.sum( scipy.special.psi(data+r) ) - n * scipy.special.psi(r) + n * np.log( r / (r+sm) )
-	
-	return np.array([f0, f1])
-
+    """Objective function for MLE estimate according to
+    https://en.wikipedia.org/wiki/Negative_binomial_distribution#Maximum_likelihood_estimation
+    
+    Parameters
+    ----------
+    par : list
+        Distribution parameters that are iteratively updated
+    data : iterable list or :class:`np.array`
+        Data points to be fit
+    sm: \sum data / len(data)
+    
+    Returns
+    -------
+    out : ndarray
+        Likelihood of parameters given data
+    """
+
+    p = par[0]
+    r = par[1]
+    n = len(data)
+    f0 = sm / (r+sm) - p
+    f1 = np.sum( scipy.special.psi(data+r) ) - n * scipy.special.psi(r) + n * np.log( r / (r+sm) )
+    
+    return np.array([f0, f1])
 
 def fit(data, p = None, r = None):
-	"""Summary
-	
-	Parameters
-	----------
-	data : TYPE
-	    Description
-	p : None, optional
-	    Description
-	r : None, optional
-	    Description
-	
-	Returns
-	-------
-	TYPE
-	    Description
-	"""
-	if p is None or r is None:
-		av = np.average(data)
-		va = np.var(data)
-		r = (av*av) / (va-av)
-		p = (va-av) / (va)
-
-	sm = np.sum(data)/len(data)
-	
-	with warnings.catch_warnings():
-		warnings.simplefilter("ignore")
-		x = scipy.optimize.fsolve(mle, np.array([p, r]), args = (data, sm))
-
-	return (x[0], x[1])
+    """Function to fit the parameters of a negative binomial to data
+    
+    Parameters
+    ----------
+    data : :class:`np.array`
+        Data to fit (1-D array)
+    p : float, optional
+        Initial guess for `p`
+    r : float, optional
+        Initial guess for `r`
+    
+    Returns
+    -------
+    out: tuple
+       Maximum likelihood fit parameters (p, r)
+    """
+    if p is None or r is None:
+        av = np.average(data)
+        va = np.var(data)
+        r = (av*av) / (va-av)
+        p = (va-av) / (va)
+
+    sm = np.sum(data)/len(data)
+    
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        x = scipy.optimize.fsolve(mle, np.array([p, r]), args = (data, sm))
+
+    return (x[0], x[1])
 
 cpdef data_type_t logpmf(int k, data_type_t p, data_type_t r) nogil:
-	cdef double coeff = c_lgamma(k+r) - c_lgamma(k+1) - c_lgamma(r)
-	return coeff + r * c_log(p) + k * c_log1p(-p)
-
+    """log of probability mass function for negative binomial
+
+    Parameters
+    ----------
+    k : int
+        Number of successes
+    p : float
+        Probability of success
+    r: float
+        Number of failures
+
+    Returns
+    -------
+    p : float
+        log of probability mass at value `k`
+    """
+    cdef double coeff = c_lgamma(k+r) - c_lgamma(k+1) - c_lgamma(r)
+    return coeff + r * c_log(p) + k * c_log1p(-p)
 
 cpdef data_type_t pmf(int k, data_type_t p, data_type_t r) nogil:
-	return c_exp( logpmf(k, p, r) )
+    """Probability mass function for negative binomial
+
+    Parameters
+    ----------
+    k : int
+        Number of successes
+    p : float
+        Probability of success
+    r: float
+        Number of failures
+
+    Returns
+    -------
+    p : float
+        Probability mass at value `k`
+    """
+    return c_exp( logpmf(k, p, r) )
 
 cpdef data_type_t cdf(int k, data_type_t p, data_type_t r) nogil:
-	return c_incbet(r, k+1, p)
+    """Cumulative distribution function for negative binomial
+
+    Parameters
+    ----------
+    k : int
+        Number of successes
+    p : float
+        Probability of success
+    r: float
+        Number of failures
+
+    Returns
+    -------
+    p : float
+        Cumulative density at value `k`
+    """
+    return c_incbet(r, k+1, p)
 
 cpdef data_type_t mean(data_type_t p, data_type_t r) nogil:
-	return p*r/(1-p)
-
-"""Not yet implemented
-def rvs(data_type_t p, data_type_t r, int n = 1):
-	return np.random.negative_binomial(r, p, n)
-"""
\ No newline at end of file
+    """Mean of negative binomial distribution
+
+    Parameters
+    ----------
+    p : float
+        Probability of success
+    r: float
+        Number of failures
+
+    Returns
+    -------
+    mu : float
+        Mean of distribution
+    """
+    return p*r/(1-p)
+
+cpdef data_type_t var(data_type_t p, data_type_t r) nogil:
+    """Variance of negative binomial distribution
+
+    Parameters
+    ----------
+    p : float
+        Probability of success
+    r: float
+        Number of failures
+
+    Returns
+    -------
+    var : float
+        Variance of distribution
+    """
+    return (p*r)/((1-p)*(1-p))
+
+cpdef data_type_t rvs(data_type_t p, data_type_t r) nogil:
+    """Sample a random variate from negative binomial distribution
+    
+    Parameters
+    ----------
+    p : float
+        Probability of success
+    r: float
+        Number of failures
+
+    Returns
+    -------
+    x : float
+        Random variate
+    """
+    raise NotImplementedError
\ No newline at end of file
diff --git a/footprint_tools/stats/distributions/t.pxd b/footprint_tools/stats/distributions/t.pxd
new file mode 100644
index 0000000..63c18e2
--- /dev/null
+++ b/footprint_tools/stats/distributions/t.pxd
@@ -0,0 +1,6 @@
+cimport cython
+cimport numpy as np
+
+ctypedef np.float64_t data_type_t
+
+cpdef data_type_t logpmf(data_type_t x, data_type_t nu, data_type_t mu, data_type_t sig2) nogil
diff --git a/footprint_tools/stats/distributions/t.pyx b/footprint_tools/stats/distributions/t.pyx
new file mode 100644
index 0000000..a811107
--- /dev/null
+++ b/footprint_tools/stats/distributions/t.pyx
@@ -0,0 +1,32 @@
+"""
+Implementation of Student's t distribution
+"""
+# cython: embedsignature=True
+
+cimport cython
+cimport numpy as np
+
+cdef extern from "cephes.h":
+    double c_lgamma(double) nogil
+    double c_log(double) nogil
+    double c_sqrt(double) nogil
+    double PI
+
+cpdef data_type_t logpmf(data_type_t x, data_type_t nu, data_type_t mu, data_type_t sig2) nogil:
+    """Log probability mass function of Student's t distribution
+    
+    Parameters
+    ----------
+    x : float
+
+    nu : float
+
+    mu : float
+
+    sig2 : float
+
+    Returns
+    -------
+    p : float
+    """
+    return c_lgamma((nu+1.0)/2.0) - c_lgamma(nu/2.0) - c_log(c_sqrt(PI*nu*sig2)) - ((nu+1.0)/2.0 * c_log(1.0+((x-mu)**2)/sig2/nu ))
diff --git a/footprint_tools/stats/fdr/__init__.py b/footprint_tools/stats/fdr/__init__.py
index 90fe05d..71949b2 100644
--- a/footprint_tools/stats/fdr/__init__.py
+++ b/footprint_tools/stats/fdr/__init__.py
@@ -1,78 +1,113 @@
-# Copyright 2015 Jeff Vierstra
-
-from .bisect import *
-
+"""
+This module contains functions to perform multiple test correction.
+"""
 import numpy as np
 import scipy.stats
 
+import operator
 
-def emperical_fdr(pvals_null, pvals):
-	"""Computes emperical FDR from a distribution of null (sampled) p-values
-
-	"""
+from ..utils import bisect
 
-	sorted_pvals_null = np.sort(np.ravel(pvals_null))
-	sorted_pvals_idx = np.argsort(pvals)
+def emperical_fdr(pvals_null, pvals):
+    """Computes emperical FDR from a distribution of null (sampled) p-values
+
+    Parameters
+    ---------
+    pvals_null : array_like
+        Sampled p-values from null expectation
+    pvals : array_like
+        Observed p-values
+    
+    Returns
+    -------
+    out : ndarray
+        Empirically adjusted p-values to be compared directly with the desired FDR
+    """
+    sorted_pvals_null = np.sort(np.ravel(pvals_null))
+    sorted_pvals_idx = np.argsort(pvals)
 
-	counts = bisect(sorted_pvals_null, pvals[sorted_pvals_idx])
-	false_positive_rates = counts / len(sorted_pvals_null)
-	false_positive_rates[false_positive_rates > 1] = 1
-	return false_positive_rates[np.argsort(sorted_pvals_idx)]
+    counts = bisect(sorted_pvals_null, pvals[sorted_pvals_idx])
+    false_positive_rates = counts / len(sorted_pvals_null)
+    false_positive_rates[false_positive_rates > 1] = 1
+    return false_positive_rates[np.argsort(sorted_pvals_idx)]
 
 # Not sure if Storey FDR is best (appears to be too conservative)
 
 def pi0est(pvals, lamb = None):
-	""" """
+    """
+    
+    Parameters
+    ---------
+    pvals : :class:`np.array'
+        Observed p-values
+    lambda : :class:`np.array`, optional
+        Lambda 
 
-	lamb = np.arange(0.05, 1, 0.05) if lamb is None else lamb
-	
-	n = len(pvals)
+    Returns
+    -------
 
-	pi0 = np.array([np.mean(pvals>=l)/(1-l) for l in lamb])
+    """
+
+    n = len(pvals)
+    lamb = np.arange(0.05, 1, 0.05) if lamb is None else lamb
+    pi0 = np.array([np.mean(pvals>=l)/(1-l) for l in lamb])
 
-	#bootstrap method
-	minpi0 = np.percentile(pi0, q = 10)
-	W = np.array([np.sum(pvals >= l) for l in lamb])
-	mse = (W / (n**2 * (1-lamb)**2)) * (1 - W/n) + (pi0 - minpi0)**2
+    #bootstrap method
+    minpi0 = np.percentile(pi0, q = 10)
+    W = np.array([np.sum(pvals >= l) for l in lamb])
+    mse = (W / (n**2 * (1-lamb)**2)) * (1 - W/n) + (pi0 - minpi0)**2
 
-	return min(pi0[mse== np.min(mse)], 1)
-	#return lamb
+    return min(pi0[mse== np.min(mse)], 1)
 
 def qvalue(pvals):
-	""" """
+    """Storey q-values
+    
+    Parameters
+    ----------
 
-	pi0 = pi0est(pvals)
+    Returns
+    -------
+
+    """
 
-	n = len(pvals)
-	u = np.argsort(pvals)
-	v = scipy.stats.rankdata(pvals, method = 'max')
-	#qvals = (pi0*n*pvals) / v
-	qvals = (pi0*n*pvals) / (v * (1 - (1-pvals)**n))
-	qvals[u[n-1]] = min(qvals[u[n-1]], 1)
-	for i in range(0, n-1)[::-1]:
-		qvals[u[i]] = min(qvals[u[i]], qvals[u[i+1]])
+    pi0 = pi0est(pvals)
 
-	return qvals
+    n = len(pvals)
+    u = np.argsort(pvals)
+    v = scipy.stats.rankdata(pvals, method = 'max')
+    #qvals = (pi0*n*pvals) / v
+    qvals = (pi0*n*pvals) / (v * (1 - (1-pvals)**n))
+    qvals[u[n-1]] = min(qvals[u[n-1]], 1)
+    for i in range(0, n-1)[::-1]:
+        qvals[u[i]] = min(qvals[u[i]], qvals[u[i+1]])
+
+    return qvals
 
 def bh_qvalue(pvals):
-    """
-    Return Benjamini-Hochberg FDR q-values corresponding to p-values C{pv}.
+    """Return Benjamini-Hochberg FDR q-values corresponding to p-values C{pv}.
 
     This function implements an algorithm equivalent to L{bh_rejected} but
     yields a list of 'adjusted p-values', allowing for rejection decisions
     based on any given threshold.
 
-    @type pv: list
-    @param pv: p-values from a multiple statistical test
+    Parameters
+    ----------
+    pvals : array_like
+        P-values from a multiple statistical test
+
+    Raises
+    ------
+    ValueError
 
-    @rtype: list
-    @return: adjusted p-values to be compared directly with the desired FDR
-      level
+    Returns
+    -------
+    out : array_like
+        Adjusted p-values to be compared directly with the desired FDR
     """
     m = len(pvals)
     rank, sorted_pvals = list(zip(*sorted(enumerate(pvals), None, operator.itemgetter(1))))
     if pvals[0] < 0 or pvals[-1] > 1:
-        raise ValueError("p-values must be between 0 and 1")
+        raise ValueError("P-values must be between 0 and 1")
     qvalues = np.zeros(m)
     mincoeff = sorted_pvals[-1]
     qvalues[rank[-1]] = mincoeff
diff --git a/footprint_tools/stats/fdr/bisect.pyx b/footprint_tools/stats/fdr/bisect.pyx
deleted file mode 100644
index aa5007b..0000000
--- a/footprint_tools/stats/fdr/bisect.pyx
+++ /dev/null
@@ -1,33 +0,0 @@
-# Copyright 2015 Jeff Vierstra
-
-# cython: boundscheck=False
-# cython: wraparound=False
-# cython: nonecheck=False
-# cython: embedsignature=True
-
-
-import numpy as np
-
-cimport cython
-cimport numpy as np
-
-ctypedef np.float64_t data_type_t
-
-cpdef bisect(data_type_t [:] a, data_type_t [:] b):
-	"""
-	"""
-	cdef int i
-	cdef int n = b.shape[0]
-
-	cdef int lo = 0
-	cdef int hi = a.shape[0]
-
-	cdef np.ndarray[data_type_t, ndim = 1] ind = np.zeros(n, dtype = np.float64)
-
-	for i in range(n):
-		while lo < hi:
-			if b[i] < a[lo]: break
-			else: lo = lo + 1
-		ind[i] = lo
-
-	return ind
\ No newline at end of file
diff --git a/footprint_tools/stats/posterior.py b/footprint_tools/stats/posterior.py
new file mode 100644
index 0000000..995745a
--- /dev/null
+++ b/footprint_tools/stats/posterior.py
@@ -0,0 +1,145 @@
+"""
+This sub-module contains functions compute the posterior probability that
+a nucleotide is footprinted.
+"""
+
+import numpy as np
+import scipy.stats
+
+from footprint_tools.stats import windowing
+
+def compute_prior_weighted(fdr, w, cutoff = 0.05, pseudocount = 0.5):
+    """Returns prior of whether a nucleotide is footprinted
+    
+    Parameters
+    ----------
+    fdr : ndarray
+        FDR values for each datset and position
+    w : ndarray
+        Binary array whether nucleotide is within a DHS
+    cutoff : float, optional
+        FPR cutoff value to label as occupied when building prior
+    pseudocount : float, optional
+        Psuedocount to add as prior to Beta distribution
+    
+    Returns
+    -------
+    :class:`np.array`
+        Per-nucleotide occupancy prior
+    """
+
+    k = np.sum(fdr <= cutoff, axis = 0) # num of datasets w/ fp
+    n = np.sum(w, axis = 0) #np.sum(w) # num of datasets w/ HS
+    a = n-k + pseudocount
+    b = k + pseudocount
+    pr = a/(a+b)
+    
+    res = np.ones(fdr.shape)
+    res *= pr[np.newaxis,:]
+    res[w==0] = 1
+
+    return res
+
+def compute_delta_prior(obs, exp, fdr, beta_prior, cutoff = 0.05):
+    """Returns a point estimate of exepected cleavage 
+    depletion with a footprint at each nucleotide 
+    
+    Parameters
+    ----------
+    obs : ndarray
+        Array with observed cleavages counts (2D - nucleotides x samples)
+    exp : ndarray
+        Array with expectd cleavages counts (2D - nucleotides x samples)
+    fdr : ndarray
+        Array with adjusted p-values (2D - nucleotides x samples)
+    beta_prior : ndarray
+        Description
+    cutoff : float, optional
+        FDR cutoff value to label a nucleotide as occupied when building prior
+    
+    Returns
+    -------
+    delta: ndarray
+        Per-nucleotide priors of expected nucleotide protection at footprinted sites 
+        (values explicitly in [0, 1]). Arrays is 1-D corresponding to DHS
+    """
+
+    (n, w) = obs.shape
+
+    mus = np.ones((n, w))
+    ws = np.ones((n, w))
+    
+    for i in range(n):
+        
+        k = obs[i,:]
+        n = np.max(np.vstack([exp[i,:], obs[i,:]]), axis = 0)
+
+        mu, v = scipy.stats.beta.stats(k+beta_prior[i][0], n-k+beta_prior[i][1], 
+            loc = 0, scale = 1, moments = 'mv')
+
+        mus[i,:] = mu
+        ws[i,:] = 1/np.sqrt(v)
+
+    ws[fdr > cutoff] = 0
+
+    delta = np.sum(ws*mus, axis = 0) / np.sum(ws, axis = 0)
+    delta[np.isnan(delta)] = 1
+
+    return delta
+
+def log_likelihood(obs, exp, dm, delta = 1, w = 3):
+    """Likelihood function of observed counts given the bias corrected data
+    and an expected protection
+    
+    Parameters
+    ----------
+    obs : ndarray
+        Observed cleavage counts
+    exp : ndarray
+        Expected cleavage counts
+    dm : :class:`dispersion_model`
+        Dispersional model to compute cleavage statistics
+    delta : array_like, optional
+        Values to scale exp counts
+    w : int, optional
+        Half window width to compute log-likelihood. Default 3bp (combined 7bp window)
+    
+    Returns
+    -------
+    log_likelihood: ndarray
+        Array of log-likelihoods for each nucloetides and samples
+    """
+    res = np.ones((obs.shape[0], obs.shape[1]), order = 'c')
+
+    n = obs.shape[0]
+    for i in range(n):
+        res[i,:] = windowing.sum(dm[i].log_pmf_values(exp[i,:] * delta, obs[i,:]), w)
+    
+    return res
+
+def posterior(prior, ll_on, ll_off):
+    """Compute the posterior probability of a nucleotide
+
+    Parameters
+    ----------
+    prior : ndarray
+        Per-nucleotide prior that a site is footprinted
+    ll_on : ndarray
+        Log-likelihood that a nucleotide is occupied
+    ll_off : ndarray
+        Log-likelihood that a nucleotide is unoccupied
+    
+    Returns
+    -------
+    posteriors : ndarray
+        log posterior that a nucleotde is footprinted in sample. Array is
+        2-D corresponding to nucleotides and samples. 
+    """
+    prior_on = np.log(1-prior)
+    prior_off = np.log(prior)
+
+    p_off = prior_off + ll_off
+    p_on = prior_on + ll_on
+    denom = np.logaddexp(p_on, p_off)
+
+    return (p_off - denom)
diff --git a/footprint_tools/stats/segment.pyx b/footprint_tools/stats/segment.pyx
deleted file mode 100644
index b5ec81b..0000000
--- a/footprint_tools/stats/segment.pyx
+++ /dev/null
@@ -1,31 +0,0 @@
-# Copyright 2015 Jeff Vierstra
-
-# cython: embedsignature=True
-
-
-
-cimport numpy as np
-
-ctypedef np.float64_t data_type_t
-
-cpdef segment(data_type_t [:] x, data_type_t threshold, int w = 1, bint decreasing = 0):
-	"""Segment an array into continuous elements passing a threshhold
-
-	"""
-
-	cdef double dir = -1 if decreasing else 1
-
-	cdef list ret = []
-	cdef int i, curr_start = -1
-	for i in range(x.shape[0]):
-		if curr_start < 0:
-			if dir*x[i] >= dir*threshold:
-				curr_start = i-w+1
-		else:
-			if dir*x[i] < dir*threshold:
-				if len(ret) > 0 and curr_start <= ret[-1][1]:
-					ret[-1][1] = i-1+w
-				else:
-					ret.append( [curr_start, i-1+w] )
-				curr_start = -1
-	return ret
diff --git a/footprint_tools/stats/utils.pyx b/footprint_tools/stats/utils.pyx
new file mode 100644
index 0000000..29c3c00
--- /dev/null
+++ b/footprint_tools/stats/utils.pyx
@@ -0,0 +1,80 @@
+"""
+This sub-module contains utility functions for stats module 
+"""
+
+# cython: boundscheck=False
+# cython: wraparound=False
+# cython: nonecheck=False
+# cython: embedsignature=True
+
+cimport numpy as np
+import numpy as np
+
+ctypedef np.float64_t data_type_t
+
+cpdef segment(data_type_t [:] x, data_type_t threshold, int w = 1, bint decreasing = 0):
+    """Segment an array into continuous elements passing a threshold
+
+    Parameters
+    ---------
+    x : array_like
+        Array of values to segment
+    threshold : float
+        Threshold for grouping elements
+    w : int
+        Window size
+    decreasing : bool
+
+    Returns
+    -------
+    out : list
+        List of tuples which specific the start and end index 
+        of contiguous intervals that pass threshold
+
+    """
+    cdef double dir = -1 if decreasing else 1
+
+    cdef list ret = []
+    cdef int i, curr_start = -1
+    for i in range(x.shape[0]):
+        if curr_start < 0:
+            if dir*x[i] >= dir*threshold:
+                curr_start = i-w+1
+        else:
+            if dir*x[i] < dir*threshold:
+                if len(ret) > 0 and curr_start <= ret[-1][1]:
+                    ret[-1][1] = i-1+w
+                else:
+                    ret.append( [curr_start, i-1+w] )
+                curr_start = -1
+    return ret
+
+cpdef bisect(data_type_t [:] a, data_type_t [:] b):
+    """Bisect arrays
+
+    Parameters
+    ----------
+    a : array_like
+        
+    b : array_like
+
+    Returns
+    -------
+    out : ndarray
+        Bisected list
+    """
+    cdef int i
+    cdef int n = b.shape[0]
+
+    cdef int lo = 0
+    cdef int hi = a.shape[0]
+
+    cdef np.ndarray[data_type_t, ndim = 1] ind = np.zeros(n, dtype = np.float64)
+
+    for i in range(n):
+        while lo < hi:
+            if b[i] < a[lo]: break
+            else: lo = lo + 1
+        ind[i] = lo
+
+    return ind
\ No newline at end of file
diff --git a/footprint_tools/stats/windowing.h b/footprint_tools/stats/windowing.h
index 2965e33..7cfe3e7 100644
--- a/footprint_tools/stats/windowing.h
+++ b/footprint_tools/stats/windowing.h
@@ -1,5 +1,4 @@
-// Copyright 2015 Jeff Vierstra
-
+// Fast windowing functions implemented in native C
 #ifndef __WINDOWING_H__
 #define __WINDOWING_H__
 
@@ -16,115 +15,115 @@ typedef double (*weighted_func_t)(const double* const, const double* const, int)
 
 double fast_sum(const double* const x, int k)
 {
-	int i;
-	double s = 0.0;
+    int i;
+    double s = 0.0;
 
-	for (i = 0; i < k; i++)
-	{
-		s += x[i];
-	}
+    for (i = 0; i < k; i++)
+    {
+        s += x[i];
+    }
 
-	return s;
+    return s;
 }
 
 double fast_product(const double* const x, int k)
 {
-	int i;
-	double p = 1.0;
+    int i;
+    double p = 1.0;
 
-	for (i = 0; i < k; i++)
-	{
-		p *= x[i];
-	}
+    for (i = 0; i < k; i++)
+    {
+        p *= x[i];
+    }
 
-	return p;
+    return p;
 }
 
 
 double fast_fishers_combined(const double* const x, int k)
 {
-	int i;
-	double p, chi = 0.0;
-
-	for (i = 0; i < k; i++)
-	{
-		chi += c_log(x[i]);
-	}
-	chi *= -2.0;
-	
-	p = c_chdtrc( (double) 2.0 * k, chi);
-	return p;
+    int i;
+    double p, chi = 0.0;
+
+    for (i = 0; i < k; i++)
+    {
+        chi += c_log(x[i]);
+    }
+    chi *= -2.0;
+    
+    p = c_chdtrc( (double) 2.0 * k, chi);
+    return p;
 }
 
 double fast_stouffers_z(const double* const x, int k)
 {
-	int i;
-	double p, z;
-	double s = 0.0;
-	
-	for (i = 0; i < k; ++i)
-	{
-		s += c_ndtri(1.0 - x[i]);
-	}
-	z = s / c_sqrt((double)k);
-	p = c_ndtr(-z);
-
-	return p;
+    int i;
+    double p, z;
+    double s = 0.0;
+    
+    for (i = 0; i < k; ++i)
+    {
+        s += c_ndtri(1.0 - x[i]);
+    }
+    z = s / c_sqrt((double)k);
+    p = c_ndtr(-z);
+
+    return p;
 }
 
 double* fast_windowing_func(const double* const x, int n, int hw, func_t func_ptr)
 {
-	int i, k = 2 * hw + 1;
-	double* tmp  = (double*) malloc(k * sizeof(double));
-	double* res = (double*) calloc(n, sizeof(double));
+    int i, k = 2 * hw + 1;
+    double* tmp  = (double*) malloc(k * sizeof(double));
+    double* res = (double*) calloc(n, sizeof(double));
 
-	for (i = hw; i < n-hw; ++i)
-	{
-		memcpy(&tmp[0], &x[i-hw], k * sizeof(double));
-		res[i] = func_ptr(tmp, k);
-	}
+    for (i = hw; i < n-hw; ++i)
+    {
+        memcpy(&tmp[0], &x[i-hw], k * sizeof(double));
+        res[i] = func_ptr(tmp, k);
+    }
 
-	free(tmp);
+    free(tmp);
 
-	return res;
+    return res;
 }
 
 double fast_weighted_stouffers_z(const double* const x, const double* const w, int k)
 {
-	int i;
-	double p, z;
-	double s = 0.0;
-	double sw = 0.0;
-
-	for (i = 0; i < k; ++i)
-	{
-		s += w[i] * c_ndtri(1.0 - x[i]);
-		sw += w[i] * w[i];
-	}
-	z = s / c_sqrt(sw);
-	p = c_ndtr(-z);
-
-	return p;	
+    int i;
+    double p, z;
+    double s = 0.0;
+    double sw = 0.0;
+
+    for (i = 0; i < k; ++i)
+    {
+        s += w[i] * c_ndtri(1.0 - x[i]);
+        sw += w[i] * w[i];
+    }
+    z = s / c_sqrt(sw);
+    p = c_ndtr(-z);
+
+    return p;	
 }
 
 double* fast_weighted_windowing_func(const double* const x, const double* const w, int n, int hw, weighted_func_t func_ptr)
 {
-	int i, k = 2 * hw + 1;
-	double* xtmp = (double*) malloc(k * sizeof(double));
-	double* wtmp = (double*) malloc(k * sizeof(double));
-	double* res = (double*) calloc(n, sizeof(double));
-
-	for (i = hw; i < n-hw; ++i)
-	{
-		memcpy(&xtmp[0], &x[i-hw], k * sizeof(double));
-		memcpy(&wtmp[0], &w[i-hw], k * sizeof(double));
-		res[i] = func_ptr(xtmp, wtmp, k);
-	}
-
-	free(xtmp);
-	free(wtmp);
-
-	return res;
+    int i, k = 2 * hw + 1;
+    double* xtmp = (double*) malloc(k * sizeof(double));
+    double* wtmp = (double*) malloc(k * sizeof(double));
+    double* res = (double*) calloc(n, sizeof(double));
+
+    for (i = hw; i < n-hw; ++i)
+    {
+        memcpy(&xtmp[0], &x[i-hw], k * sizeof(double));
+        memcpy(&wtmp[0], &w[i-hw], k * sizeof(double));
+        res[i] = func_ptr(xtmp, wtmp, k);
+    }
+
+    free(xtmp);
+    free(wtmp);
+
+    return res;
 }
 
 #endif
\ No newline at end of file
diff --git a/footprint_tools/stats/windowing.pyx b/footprint_tools/stats/windowing.pyx
index dd8589d..ca9fb20 100644
--- a/footprint_tools/stats/windowing.pyx
+++ b/footprint_tools/stats/windowing.pyx
@@ -1,16 +1,18 @@
-# Copyright 2015 Jeff Vierstra
+"""
+This sub-module contains functions for performing sliding window functions
+written in native C. 
+"""
+
+# Copyright 2015-2021 Jeff Vierstra
 
 # cython: boundscheck=False
 # cython: wraparound=False
 # cython: nonecheck=False
 # cython: embedsignature=True
 
-
-from libc.stdlib cimport free
-
-
 import numpy as np
 
+from libc.stdlib cimport free
 cimport numpy as np
 
 ctypedef np.float64_t data_type_t
@@ -30,7 +32,18 @@ cdef extern from "windowing.h":
     double* fast_weighted_windowing_func(const double* const, const double* const, int, int, weighted_func_t)
 
 cdef windowing_func(data_type_t [:] x, int hw, func_t func_ptr):
-    
+    """Windowing function wrapper that passes ndarray to
+    fast functions written in native C
+
+    Parameters
+    ---------
+    x : ndarray
+        Array of values to perform sliding window function
+    hw : int
+        Half window width
+    func_ptr: callable
+        Function to apply to each window 
+    """
     cdef int i
     cdef int n = x.shape[0]
 
@@ -45,52 +58,92 @@ cdef windowing_func(data_type_t [:] x, int hw, func_t func_ptr):
     return np.asarray(win)
 
 cpdef sum(data_type_t [:] x, int hw):
-    """Sum of values in sliding local window
-
-    :param x: (array) values
-    :param hw: (int) half window size to smooth
-    
-    :return: (float) summed values
+    """Sum of values in sliding local window across array
+
+    Paramters
+    ---------
+    x : ndarray
+        Values to perform windowed sum. Must be memory
+        contiguous (i.e., np.ascontiguousarray(x))
+    hw: int
+        Half window with to apply sum
+
+    Returns
+    -------
+    out : ndarray
+       Array of windowed summed values
     """
-
     return windowing_func(x, hw, fast_sum)
 
 cpdef product(data_type_t [:] x, int hw):
-    """Product of values in sliding local window
-
-    :param x: (array) values
-    :param hw: (int) half window size to smooth
-    
-    :return: (float) product of values
+    """Sum of values in sliding local window across array
+
+    Paramters
+    ---------
+    x : ndarray
+        Values to perform windowed product. Must be memory
+        contiguous (i.e., np.ascontiguousarray(x))
+    hw: int
+        Half window with to apply sum
+
+    Returns
+    -------
+    out : ndarray
+       Array of windowed product values
     """
-    
     return windowing_func(x, hw, fast_product)
 
 cpdef fishers_combined(data_type_t [:] x, int hw):
-    """Compute p-value for a window using Fisher's combined method
-
-    :param x: (array) p-values
-    :param hw: (int) half window size to smooth
-    
-    :return p: (float) combined p-values
+    """Compute p-values for a window using Fisher's combined method
+
+    Paramters
+    ---------
+    x : ndarray
+        P-values to perform Fisher's combined method. Array must 
+        be memory contiguous (i.e., np.ascontiguousarray(x))
+    hw: int
+        Half window with to combined p-values
+
+    Returns
+    -------
+    out : ndarray
+       Array of combined p-values
     """
-
-
     return windowing_func(x, hw, fast_fishers_combined)
 
 cpdef stouffers_z(data_type_t [:] x, int hw):
-    """Compute p-value for a window using Stouffer's method
-
-    :param x: (array) p-values
-    :param hw: (int) half window size to smooth
-    
-    :return p: (float) combined p-values
+    """Compute p-values for a window using Stouffer's Z-score method
+
+    Paramters
+    ---------
+    x : ndarray
+        P-values to perform Stouffer's Z-score method. Array must 
+        be memory contiguous (i.e., np.ascontiguousarray(x))
+    hw: int
+        Half window with to combined p-values
+
+    Returns
+    -------
+    out : ndarray
+       Array of combined p-values
     """
-
     return windowing_func(x, hw, fast_stouffers_z)
 
 cdef weighted_windowing_func(data_type_t [:] x, data_type_t [:] w, int hw, weighted_func_t func_ptr):
-    
+    """Weighted windowing function wrapper that passes ndarray to
+    fast functions written in native C
+
+    Parameters
+    ---------
+    x : ndarray
+        Array of values to perform sliding window function
+    w: ndarray
+        Weights for each array element
+    hw : int
+        Half window width
+    func_ptr: callable
+        Function to apply to each window 
+    """
     cdef int i
     cdef int n = x.shape[0]
 
@@ -107,11 +160,19 @@ cdef weighted_windowing_func(data_type_t [:] x, data_type_t [:] w, int hw, weigh
 cpdef weighted_stouffers_z(data_type_t [:] x, data_type_t [:] w, int hw):
     """Compute p-value for a window using weighted Stouffer's method
 
-    :param x: (array) p-values
-    :param w: (array) weights
-    :param hw: (int) half window size to smooth
-    
-    :return p: (float) combined p-values
+    Parameters
+    ----------
+    x : ndarray
+        P-values to perform Stouffer's Z-score method. Array must 
+        be memory contiguous (i.e., np.ascontiguousarray(x))
+    w: ndarray
+        Weights for each elements in `x`
+    hw: int
+        Half window with to combined p-values
+
+    Returns
+    -------
+    out : ndarray
+       Array of combined p-values
     """
-
     return weighted_windowing_func(x, w, hw, fast_weighted_stouffers_z)
diff --git a/requirements.txt b/requirements.txt
index 7fa3e6f..3a560e3 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -7,4 +7,8 @@ pwlf
 simplejson
 cython
 genome_tools>=1.0.2
-tqdm
\ No newline at end of file
+tqdm
+click
+click-option-group
+matplotlib
+yaspin
\ No newline at end of file
diff --git a/scripts/ftd-compute-deviation b/scripts/ftd-compute-deviation
deleted file mode 100644
index 2b33b2b..0000000
--- a/scripts/ftd-compute-deviation
+++ /dev/null
@@ -1,285 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2015-2021 Jeff Vierstra (jvierstra@altius.org)
-
-import sys
-import os
-import logging
-import random
-import math
-
-from argparse import ArgumentParser, Action, ArgumentError
-
-import multiprocessing as mp
-from tqdm import tqdm
-
-import numpy as np
-
-import pysam
-
-from genome_tools import bed, genomic_interval
-
-import footprint_tools
-from footprint_tools import cutcounts
-from footprint_tools.modeling import bias, predict, dispersion
-from footprint_tools.stats import fdr, windowing
-
-logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%m/%d/%Y %I:%M:%S %p', level=logging.INFO)
-logger = logging.getLogger(__name__)
-
-class bias_model_action(Action):
-    def __call__(self, parser, namespace, values, option_string = None):
-        try:
-            setattr(namespace, self.dest, bias.kmer_model(values[0]))
-        except IOError as e:
-             raise ArgumentError(self, str(e))
-
-class dispersion_model_action(Action):
-    def __call__(self, parser, namespace, values, option_string = None):
-        try:
-            setattr(namespace, self.dest, dispersion.read_dispersion_model(values[0]))
-        except IOError as e:
-             raise ArgumentError(self, str(e))
-
-def paired_ints(arg):
-    """
-    Parse a pair of integers delimited by a comma (used by argparser)
-    """
-
-    try:
-        fw, rev = list(map(int, arg.split(',')))
-        return (fw, rev)
-    except:
-        raise ArgumentTypeError("Offset argument must be a in the format of int,int")
-
-def parse_options(args):
-
-    parser = ArgumentParser(description = "Compute the per-nucleotide cleavage deviation statistics (version {})".format(footprint_tools.__version__),
-                epilog="Written by Jeff Vierstra (jvierstra@altius.org). See http://github.com/jvierstra/footprint-tools for extended documentation.")
-
-    parser.add_argument("bam_file", metavar = "bam_file", type = str,
-                        help = "Path to BAM-format tag sequence file")
-
-    parser.add_argument("fasta_file", metavar = "fasta_file", type = str, 
-                        help = "Path to genome FASTA file (requires associated"
-                        " FASTA index in same folder; see documentation on how"
-                        " to create an index)")
-
-    parser.add_argument("interval_file", metavar = "interval_file", type = str, 
-                        help = "File path to BED file")
-
-    grp_bm = parser.add_argument_group("Bias modeling options")
-
-    grp_bm.add_argument("--bm", metavar = "MODEL_FILE", dest = "bias_model", 
-                        nargs = 1, action = bias_model_action, default = bias.uniform_model(),
-                        help = "Use a k-mer model for local bias (supplied by file). If"
-                        " argument is not provided the model defaults to uniform sequence"
-                        " bias.")
-
-    grp_bm.add_argument("--half-win-width", metavar = "N", type = int, default = 5,
-                        help = "Half window width to apply bias model."
-                        " (default: %(default)s)")
-
-    grp_sm = parser.add_argument_group("Smoothing options")
-    
-    grp_sm.add_argument("--smooth-half-win-width", metavar = "N", type = int, default = 50,
-                        help = "Half window width to apply smoothing model. When set to"
-                        " zero no smoothing is applied. (default: %(default)s)")
-
-    grp_sm.add_argument("--smooth-clip", metavar = "N", type = float, default = 0.01,
-                        help = "Fraction of signal to clip when computing trimmed mean."
-                        " (default: %(default)s)")
-
-    grp_st = parser.add_argument_group("Statistics options")
-
-    grp_st.add_argument("--dm", nargs = 1, metavar = "MODEL_FILE", 
-                        dest = "dispersion_model", action = dispersion_model_action, default = None,
-                        help = "Dispersion model for negative binomial tests. If argument"
-                        " is not provided then no stastical output is provided. File is in"
-                        " JSON format and generated using the 'ftd-learn-dispersion-model'"
-                        " script included in the software package.")
-
-    # Not yet implemented
-    # grp_st.add_argument("--fdr_cutoffs", nargs = 1, metavar = "[N, ...]", dest = "fdr_cutoffs",
-    #                     action = list_parse_action, default = [],
-    #                     help = "FDR cutoff at which to report footprints.")
-
-    grp_st.add_argument("--fdr-shuffle-n", metavar = "N", type = int,
-                        dest = "fdr_shuffle_n", default = 50,
-                        help = "Number of times to shuffle data for FDR calculation."
-                        " (default: %(default)s)")
-
-    grp_rf = parser.add_argument_group("Read filtering options")
-
-    grp_rf.add_argument("--min-qual", metavar = "N", type = int,
-                        dest = "min_qual", default = 1,
-                        help = "Filter reads with mapping quality lower than this threshold."
-                        " (default: %(default)s)")
-
-    grp_rf.add_argument("--remove-dups", action = "store_true",
-                        dest = "remove_dups", help = "Remove duplicate reads from analysis"
-                        " (SAM flag -- 512)")
-    grp_rf.add_argument("--remove-qcfail", action = "store_true",
-                        dest = "remove_qcfail", help = "Remove QC-failed reads from analysis"
-                        " (SAM flag -- 1024)")
-
-    grp_ot = parser.add_argument_group("Other options")
-
-    grp_ot.add_argument("--bam-offset", metavar = "N", type = paired_ints,
-                        dest = "bam_offset", default = (0, -1),
-                        help = "BAM file offset (support for legacy BAM/SAM format)"
-                        " (default: %(default)s)")
-
-    grp_ot.add_argument("--seed", metavar = "N", type = int, dest = "seed",
-                        default = None, help = "Seed for random number generation"
-                        " (default: no seed)")
-
-    grp_ot.add_argument("--processors", metavar = "N", type = int,
-                        dest = "processors", default = max(1, mp.cpu_count()-1),
-                        help = "Number of processors to use."
-                        " (default: all available processors)")
-
-    return parser.parse_args(args)
-
-def read_func(bam_file, fasta_file, bm, dm, intervals, q, **kwargs):
-    """
-    Reads BAM file, computes expected cleavages and associated statistics and outputs to a multiprocessing pool queue
-    """
-
-    bam_kwargs = { k:kwargs.pop(k) for k in ["min_qual", "remove_dups", "remove_qcfail", "offset"] }
-    bam_reader = cutcounts.bamfile(bam_file, **bam_kwargs)
-    fasta_reader = pysam.FastaFile(fasta_file)
-
-    predict_kwargs = { k:kwargs.pop(k) for k in ["half_window_width", "smoothing_half_window_width", "smoothing_clip"] }
-    predictor = predict.prediction(bam_reader, fasta_reader, bm, **predict_kwargs)
-
-    # args used for FDR sampling procedure
-    seed = kwargs.pop("seed")
-    if seed:
-        random.seed(seed)
-        np.random.seed(seed)
-    fdr_shuffle_n = kwargs.pop("fdr_shuffle_n") 
-
-    win_pvals_func = lambda z: windowing.stouffers_z(np.ascontiguousarray(z), 3)
-    
-    # Read and process each region
-    for interval in intervals:
-
-        obs, exp, win = predictor.compute(interval)
-
-        obs = obs['+'][1:] + obs['-'][:-1]
-        exp = exp['+'][1:] + exp['-'][:-1]
-
-        n = len(obs)
-
-        if dm:
-
-            try:
-
-                pvals = dm.p_values(exp, obs)
-                counts_null, pvals_null = dm.resample_p_values(exp, fdr_shuffle_n)
-
-                win_pvals = win_pvals_func(pvals)
-                win_pvals_null = np.apply_along_axis(win_pvals_func, 0, pvals_null)
-                
-                efdr = fdr.emperical_fdr(win_pvals_null, win_pvals)
-
-            except Exception as e:
-
-                pvals = win_pvals = efdr = np.ones(n)
-
-            finally:
-
-                stats = np.column_stack((exp, obs, -np.log(pvals), -np.log(win_pvals), efdr))
-        else:
-
-            stats = np.column_stack((exp, obs))
-
-        q.put( (interval, stats) )
-
-        while q.qsize() > 100:
-            pass
-
-def write_func(q, total):
-    """
-    Function to write output from multiple threads to a single file
-    """    
-
-    f = sys.stdout
-
-    progress_desc="Regions processed"
-    with tqdm(total=total, desc=progress_desc, ncols=80) as progress_bar:
-    
-        while True:
-
-            data = q.get()
-            if data == None:
-                q.task_done()
-                break
-
-            interval, stats = data
-
-            for i in range(stats.shape[0]):
-                coords = "{}\t{:d}\t{:d}".format(interval.chrom, interval.start+i, interval.start+i+1)
-                val_string = "\t".join( ["{:0.4f}".format(val) for val in stats[i,:]])
-                print(coords + "\t" + val_string, file=f)
-
-            q.task_done()
-
-            progress_bar.update(1)
-
-def chunkify(l, nchunks):
-    """
-    Splits an iterable list into n chunks
-    """
-    
-    return [l[i::nchunks] for i in range(nchunks)]
-
-def main(argv = sys.argv[1:]):
-
-    args = parse_options(argv)
-
-    intervals = list(bed.bed3_iterator(open(args.interval_file)))
-    
-    logger.info("BED file contains {:,} regions".format(len(intervals)))
-
-    proc_kwargs = {
-        "min_qual": args.min_qual,
-        "remove_dups": args.remove_dups,
-        "remove_qcfail": args.remove_qcfail,
-        "offset": args.bam_offset,
-        "half_window_width": args.half_win_width,
-        "smoothing_half_window_width": args.smooth_half_win_width,
-        "smoothing_clip": args.smooth_clip,
-        "fdr_shuffle_n": args.fdr_shuffle_n,
-        "seed": args.seed
-        }
-
-    q = mp.JoinableQueue()
-
-    read_procs = []
-    for i, chunk in enumerate(chunkify(intervals, max(1, args.processors-1))):
-           p = mp.Process(target=read_func, args=(args.bam_file, args.fasta_file, args.bias_model, args.dispersion_model, chunk, q), kwargs=proc_kwargs)
-           read_procs.append(p)
-
-    write_proc = mp.Process(target=write_func, args=(q, len(intervals)))
-
-    logger.info("Using {} threads to compute footprint statistics".format(len(read_procs)))
-
-    [p.start() for p in read_procs]
-    write_proc.start()
-
-    try:
-        [p.join() for p in read_procs]
-        q.join() # block until queue is empty after processing is done
-        q.put(None) # send kill signal
-        write_proc.join() # wait for writer proc to return
-        logger.info("Finished computing and writing footprint statistics!")
-    except KeyboardInterrupt:
-        [p.terminate() for p in read_procs]
-        write_proc.terminate()
-        
-    return 0
-    
-if __name__ == "__main__":
-    sys.exit(main())
diff --git a/scripts/ftd-compute-posterior b/scripts/ftd-compute-posterior
deleted file mode 100644
index e953d9f..0000000
--- a/scripts/ftd-compute-posterior
+++ /dev/null
@@ -1,262 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2015-2021 Jeff Vierstra (jvierstra@altius.org)
-
-import sys, logging, os, os.path
-from argparse import ArgumentParser
-import multiprocessing as mp
-
-import numpy as np
-import pandas as pd
-
-import scipy.stats
-import pysam
-
-from genome_tools import bed, genomic_interval, genomic_interval_set
-
-import footprint_tools
-from footprint_tools.modeling import dispersion
-from footprint_tools.stats import bayesian, segment
-
-logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%m/%d/%Y %I:%M:%S %p', level=logging.INFO)
-logger = logging.getLogger(__name__)
-
-def parse_options(args):
-
-    parser = ArgumentParser(description = "Compute the posterior probability of cleavage data (version {})".format(footprint_tools.__version__),
-                epilog="Written by Jeff Vierstra (jvierstra@altius.org). See http://github.com/jvierstra/footprint-tools for extended documentation.")
-
-    parser.add_argument("metadata_file", metavar = "metadata_file", type = str,
-                        help = "Path to configuration file containing metadata for samples")
-
-    parser.add_argument("interval_file", metavar = "interval_file", type = str, 
-                        help = "File path to BED file")
-
-    grp_st = parser.add_argument_group("Statistical options")
-
-    grp_st.add_argument("--fdr-cutoff", metavar = "N", type = float,
-                        dest = "fdr_cutoff", default = 0.05,
-                        help = "Only consider nucleotides with a minimum FDR <= this value."
-                        " (default: %(default)s)")
-
-    grp_o = parser.add_argument_group("Output options")
-
-    grp_o.add_argument("--post-cutoff", metavar = "N", type = float,
-                        dest = "post_cutoff", default = 0.2,
-                        help = "Only output nucleotides with posterior probability <= this value in 1 or more samples"
-                        " (default: %(default)s)")
-
-    grp_ot = parser.add_argument_group("Other options")
-
-    grp_ot.add_argument("--processors", metavar = "N", type = int,
-                        dest = "processors", default = max(1, mp.cpu_count()-3),
-                        help = "Number of processors to use. Note that value excludes the"
-                        " minimum 3 threads that are dedicated to data I/O"
-                        " (default: all available processors)")
-
-    return parser.parse_args(args)
-
-def read_func(tabix_files, intervals, queue):
-    """
-    Reads TABIX files and outputs to a multiprocessing pool queue
-    """
-
-    tabix_handles = [pysam.TabixFile(f) for f in tabix_files]
-    n_datasets = len(tabix_handles)
-
-     # Write to input queue
-    for interval in intervals:
-
-        l = len(interval)
-
-        obs = np.zeros((n_datasets, l), dtype = np.float64)
-        exp = np.zeros((n_datasets, l), dtype = np.float64)
-        fdr = np.ones((n_datasets, l), dtype = np.float64)
-        w = np.zeros((n_datasets, l), dtype = np.float64)
-
-        i = 0
-        j = 0
-    
-        for tabix in tabix_handles:
-
-            try:
-                for row in tabix.fetch(interval.chrom, interval.start, interval.end, parser = pysam.asTuple()):
-                    j = int(row[1])-interval.start
-                    exp[i, j] = np.float64(row[3])
-                    obs[i, j] = np.float64(row[4])
-                    fdr[i, j] = np.float64(row[7])
-                    w[i, j] = 1
-            except:
-                pass
-
-            i += 1
-
-        queue.put( (interval, exp, obs, fdr, w) )
-
-        # Stop memory from getting out of control in the processing
-        while queue.qsize() > 100:
-            pass
-
-    [handle.close() for handle in tabix_handles]
-
-def process_func(disp_models, beta_priors, read_q, write_q, fdr_cutoff=0):
-    """
-    
-    """
-    
-    while True:
-
-        data = read_q.get()
-
-        if data == None:
-            read_q.task_done()
-            break
-
-        (interval, exp, obs, fdr, w) = data
-
-        prior = bayesian.compute_prior_weighted(fdr, w, cutoff = fdr_cutoff)   
-        scale = bayesian.compute_delta_prior(obs, exp, fdr, beta_priors, cutoff = fdr_cutoff)
-
-        ll_on = bayesian.log_likelihood(obs, exp, disp_models, delta = scale, w = 3) 
-        ll_off = bayesian.log_likelihood(obs, exp, disp_models, w = 3)
-
-        # Compute posterior
-        post = -bayesian.posterior(prior, ll_on, ll_off)
-        post[post <= 0] = 0.0
-
-        read_q.task_done()
-
-        write_q.put( (interval, post) )
-
-def write_func(q, total, log_post_cutoff=0):
-    """
-    
-    """
-    
-    handle = sys.stdout
-
-    progress_desc="Regions processed"
-    with tqdm(total=total, desc=progress_desc, ncols=80) as progress_bar:
-
-        while True:
-
-            data = q.get()
-            if data == None:
-                q.task_done()
-                break
-
-            interval, post = data
-            
-            # Write ouput; filter out positions by posterior cutoff
-            for i in np.where(np.nanmax(post, axis = 0) > log_post_cutoff)[0]:
-                print("{}\t{:d}\t{:d}\t".format(interval.chrom, interval.start+i, interval.start+i+1) + '\t'.join(map(str, post[:,i])), file = handle)
-
-            q.task_done()
-
-            progress_bar.update(1)
-
-def chunkify(l, nchunks):
-    """
-    Splits an iterable list into n chunks
-    """
-    
-    return [l[i::nchunks] for i in range(nchunks)]
-
-def main(argv = sys.argv[1:]):
-
-    args = parse_options(argv)
-
-    logger.info("Reading metadata file and verifying inputs")
-
-    metadata = pd.read_table(args.metadata_file, header=0, comment='#')
-
-    disp_models = []
-    beta_priors = np.ones((len(metadata), 2))
-
-    # Load and parse input files
-    for sample in metadata.itertuples():
-        if not os.path.exists(sample["tabix_file"]):
-            logger.critical("Fatal error: TABIX-file ({}) does not exists!".format(sample["tabix_file"]))
-            sys.exit(1)
-
-        if not os.path.exists(sample["dispersion_model"]):
-            logger.critical("Fatal error: Dispersion model ({}) does not exists!".format(sample["dispersion_model"]))
-            sys.exit(1)
-        else:
-            disp_models.append(dispersion.read_dispersion_model(sample["dispersion_model"]))
-       
-        if not os.path.exists(sample["beta_prior_file"]):
-            logger.critical("Fatal error: Beta-prior file ({}) does not exists!".format(sample["beta_prior_file"]))
-            sys.exit(1)
-        else:
-            with open(sample["beta_prior_file"], 'r') as f:
-                params = f.readline().strip().split('\t')
-                beta_priors[i,:] = np.array(params, dtype = np.float64)
-
-    # Load intervals file
-    intervals = list(bed.bed3_iterator(open(args.interval_file)))
-    
-    logger.info("BED file contains {:,} regions".format(len(intervals)))
-    
-    # Processing queues
-    read_q = mp.JoinableQueue()
-    write_q = mp.JoinableQueue()
-
-    #
-    read_procs = []
-    for i, chunk in enumerate(chunkify(intervals, 2))):
-        p = mp.Process(target=read_func, args=(metadata["tabix_file"].tolist(), chunk, read_q))
-        read_procs.append(p)
-
-    logger.info("Using {} threads to read footprint statistics".format(len(read_procs)))
-    [p.start() for p in read_procs]
-
-    #
-    process_procs = []
-    process_kwargs = {
-        "fdr_cutoff": args.fdr_cutoff,
-    }
-
-    for i in range(args.processors):
-        p = mp.Process(target=process_func, args=(disp_models, beta_priors, read_q, write_q), kwargs=process_kwargs)
-        process_procs.append(p)
-
-    logger.info("Using {} threads to compute posteriors".format(len(process_procs)))
-    [p.start() for p in process_procs]
-
-    #
-    write_kwargs = {
-        "log_post_cutoff": -np.log(args.post_cutoff)
-    }
-    write_proc = mp.Process(target=write_func, args=(write_q, len(intervals)), kwargs=write_kwargs)
-    write_proc.start()
-
-    try:
-        # Wait for readers to finish
-        [p.join() for p in read_procs]
-
-        # Block until all remaining regions are processed
-        read_q.join() # wait till read queue is empty
-
-        # Send a message to kill processing threads	and block until they exit
-        [read_q.put(None) for i in range(len(process_procs))]
-        [p.join() for p in process_procs]
-        
-        # Block until all remaining regions are written to output and thread exits
-        write_q.join() # wait till write queue is empty
-        write_q.put(None) # sends kill signal
-        write_proc.join() # block until thread exits
-
-        logger.info("Finished computing and writing footprint posteriors!")
-
-    except KeyboardInterrupt:
-        [p.terminate() for p in read_procs]
-        [p.terminate() for p in process_procs]
-        write_proc.terminate()
-
-    return 0
-
-if __name__ == "__main__":
-
-    sys.exit(main())
-
diff --git a/scripts/ftd-learn-beta-prior b/scripts/ftd-learn-beta-prior
deleted file mode 100644
index c34659f..0000000
--- a/scripts/ftd-learn-beta-prior
+++ /dev/null
@@ -1,78 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2015-2021 Jeff Vierstra (jvierstra@altius.org)
-
-import sys, os, logging
-from argparse import ArgumentParser
-import multiprocessing as mp
-from functools import partial
-
-import numpy as np
-import scipy.stats
-
-import pysam
-
-import footprint_tools 
-
-logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%m/%d/%Y %I:%M:%S %p', level=logging.INFO)
-logger = logging.getLogger(__name__)
-
-def parse_options(args):
-
-    parser = ArgumentParser(description = "Learn the genome-wide beta model prior for nucleotide protection (version {})".format(footprint_tools.__version__),
-                epilog="Written by Jeff Vierstra (jvierstra@altius.org). See http://github.com/jvierstra/footprint-tools for extended documentation.")
-
-    parser.add_argument("tabix_file", metavar = "tabix_file", type = str,
-                        help = "Path to TABIX-format file (e.g., ouput"
-                        " from 'ftd-compute-deviation')")
-
-    grp_m = parser.add_argument_group("Model options")
-
-    grp_m.add_argument("--fdr-cutoff", metavar = "N", type = float,
-                        dest = "fdr_cutoff", default = 0.05,
-                        help = "Only consider nucleotides with FDR <= this value."
-                        " (default: %(default)s)")
-
-    grp_m.add_argument("--exp-cutoff", metavar = "N", type = int,
-                        dest = "exp_cutoff", default = 10,
-                        help = "Only consider nucleotides with expected cleavages >= this value."
-                        " (default: %(default)s)")
-
-    return parser.parse_args(args)
-
-##
-
-def compute_beta_prior(filename, fdr_cutoff, exp_cutoff):
-
-	handle = pysam.TabixFile(filename)
-
-	ratios = []
-
-	for row in handle.fetch(parser = pysam.asTuple()):
-
-		exp = np.float64(row[3])
-		obs = np.float64(row[4])
-		fdr = np.float64(row[7])
-
-		if fdr <= fdr_cutoff and exp >= exp_cutoff:
-			ratios.append( (obs+1)/(exp+1) )
-
-	handle.close()
-
-	ratios = np.array(ratios)
-	(a, b) = scipy.stats.beta.fit(ratios[np.isfinite(ratios) & (ratios>0) & (ratios<1)], floc = 0, fscale = 1)[0:2]
-
-	return [a, b]
-
-def main(argv = sys.argv[1:]):
-
-	args = parse_options(argv)
-
-	prior = compute_beta_prior(args.tabix_file, fdr_cutoff = args.fdr_cutoff, exp_cutoff = args.exp_cutoff)
-
-	print("%0.4f\t%0.4f" % (prior[0], prior[1]), file = sys.stdout)
-
-	return 0
-
-if __name__ == "__main__":
-    sys.exit(main())
diff --git a/scripts/ftd-learn-dispersion-model b/scripts/ftd-learn-dispersion-model
deleted file mode 100644
index 4dbe632..0000000
--- a/scripts/ftd-learn-dispersion-model
+++ /dev/null
@@ -1,203 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2015-2021 Jeff Vierstra (jvierstra@altius.org)
-
-import sys
-import os
-import logging
-import random
-import math
-
-from argparse import ArgumentParser, Action, ArgumentError, ArgumentTypeError
-
-import multiprocessing as mp
-from tqdm import tqdm
-
-import numpy as np
-
-import pysam
-
-from genome_tools import bed, genomic_interval
-
-import footprint_tools
-from footprint_tools import cutcounts
-from footprint_tools.modeling import bias, predict, dispersion
-
-logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%m/%d/%Y %I:%M:%S %p', level=logging.INFO)
-logger = logging.getLogger(__name__)
-
-class bias_model_action(Action):
-    def __call__(self, parser, namespace, values, option_string = None):
-        try:
-            setattr(namespace, self.dest, bias.kmer_model(values[0]))
-        except IOError as e:
-             raise ArgumentError(self, str(e))
-
-def paired_ints(arg):
-	try:
-		fw, rev = list(map(int, arg.split(',')))
-		return (fw, rev)
-	except:
-		raise ArgumentTypeError("Offset argument must be a in the format of int,int")
-
-def parse_options(args):
-
-    parser = ArgumentParser(description = "Learn a negative binomial dispersion model from data corrected for intrinsic sequence preference (version {})".format(footprint_tools.__version__),
-                epilog="Written by Jeff Vierstra (jvierstra@altius.org). See http://github.com/jvierstra/footprint-tools for extended documentation.")
-
-    parser.add_argument("bam_file", metavar = "bam_file", type = str,
-                        help = "Path to BAM-format tag sequence file")
-
-    parser.add_argument("fasta_file", metavar = "fasta_file", type = str, 
-                        help = "Path to genome FASTA file (requires associated"
-                        " FASTA index in same folder; see documentation on how"
-                        " to create an index)")
-
-    parser.add_argument("interval_file", metavar = "interval_file", type = str, 
-                        help = "File path to BED file")
-
-    grp_bm = parser.add_argument_group("Bias modeling options")
-
-    grp_bm.add_argument("--bm", metavar = "MODEL_FILE", dest = "bias_model", 
-                        nargs = 1, action = bias_model_action, default = bias.uniform_model(),
-                        help = "Use a k-mer model for local bias (supplied by file). If"
-                        " argument is not provided the model defaults to uniform sequence"
-                        " bias.")
-
-    grp_bm.add_argument("--half-win-width", metavar = "N", type = int, default = 5,
-                        help = "Half window width to apply bias model."
-                        " (default: %(default)s)")
-
-    grp_rf = parser.add_argument_group("Read filtering options")
-
-    grp_rf.add_argument("--min-qual", metavar = "N", type = int,
-                        dest = "min_qual", default = 1,
-                        help = "Ignore reads with mapping quality lower than this threshold."
-                        " (default: %(default)s)")
-
-    grp_rf.add_argument("--remove-dups", action = "store_true",
-                        dest = "remove_dups", help = "Remove duplicate reads from analysis"
-                        " (SAM flag -- 512)")
-
-    grp_rf.add_argument("--remove-qcfail", action = "store_true",
-                        dest = "remove_qcfail", help = "Remove QC-failed reads from analysis"
-                        " (SAM flag -- 1024)")
-
-    grp_ot = parser.add_argument_group("Other options")
-
-    grp_ot.add_argument("--bam-offset", metavar = "N", type = paired_ints,
-                        dest = "bam_offset", default = (0, -1),
-                        help = "BAM file offset (enables support for other datatypes -- e.g. Tn5/ATAC)"
-                        " (default: %(default)s)")
-
-    grp_ot.add_argument("--processors", metavar = "N", type = int,
-                        dest = "processors", default = mp.cpu_count(),
-                        help = "Number of processors to use."
-                        " (default: all available processors)")
-
-    return parser.parse_args(args)
-
-
-
-class process_callback(object):
-	"""
-	Class that aggregates results from processor functions
-	In this case, it just sums up the 2d histogram matricies
-	"""
-
-	def __init__(self, hist_size):
-		self.x = np.zeros(hist_size)
-
-	def __call__(self, other):
-		self.x += other
-
-def process_func(bam_file, fasta_file, bm, intervals, hist_size, proc_id, **kwargs):
-	"""
-	Function that can be used to parallize reading 
-	the BAM files and computing expected counts
-	"""
-
-	bam_kwargs = { k:kwargs.pop(k) for k in ["min_qual", "remove_dups", "remove_qcfail", "offset"] }
-	bam_reader = cutcounts.bamfile(bam_file, **bam_kwargs)
-	fasta_reader = pysam.FastaFile(fasta_file)
-
-	predictor = predict.prediction(bam_reader, fasta_reader, bm, **kwargs)
-
-	res_hist = np.zeros(hist_size)
-
-	progress_desc = "Chunk {}".format(proc_id)
-
-	with tqdm(total=len(intervals), desc=progress_desc, position=proc_id+1, ncols=80) as progress_bar:
-
-		for interval in intervals:
-
-			obs, exp, win = predictor.compute(interval)
-			
-			obs = obs['+'][1:] + obs['-'][:-1]
-			exp = exp['+'][1:] + exp['-'][:-1]
-		
-			for o, e in zip(obs, exp):
-				try:
-					res_hist[int(e), int(o)] += 1.0
-				except IndexError:
-					pass
-
-			progress_bar.update(1)
-
-	return res_hist
-
-def chunkify(l, nchunks):
-	"""Splits an iterable list into n chunks"""
-	return [l[i::nchunks] for i in range(nchunks)]
-
-def main(argv = sys.argv[1:]):
-
-	args = parse_options(argv)
-
-	hist_size = (200, 1000) # hard coded histogram size -- for now...
-	hist_agg = process_callback(hist_size)
-
-	intervals = list(bed.bed3_iterator(open(args.interval_file)))
-
-	logger.info("BED file contains {:,} regions".format(len(intervals)))
-
-	proc_kwargs = {
-		"min_qual": args.min_qual,
-		"remove_dups": args.remove_dups,
-		"remove_qcfail": args.remove_qcfail,
-		"offset": args.bam_offset,
-		"half_window_width": args.half_win_width,
-		}
-
-	pool = mp.Pool(args.processors)
-
-	logger.info("Using {} threads to computed expected cleavage counts".format(args.processors))
-
-	for i, chunk in enumerate(chunkify(intervals, args.processors)):
-		pool.apply_async(process_func, args=(args.bam_file, args.fasta_file, args.bias_model, chunk, hist_size, i), kwds=proc_kwargs, callback=hist_agg)
-
-	pool.close()
-	pool.join()
-
-	logger.info("Finished computing expected cleavage counts!")
-
-	# Learn model from histogram
-
-	logger.info("Learning dispersion model")
-
-	model = dispersion.learn_dispersion_model(hist_agg.x)
-
-	# Write model
-
-	model_file = os.path.abspath(os.path.join(os.getcwd(), "dm.json"))
-
-	logger.info("Writing dispersion model to {}".format(model_file))
-
-	with open(model_file, "w") as f:
-		print(dispersion.write_dispersion_model(model), file = f)
-
-	# Success!
-	return 0
-
-if __name__ == "__main__":
-	sys.exit(main())
diff --git a/setup.py b/setup.py
index fee8b3d..8ca0cd8 100644
--- a/setup.py
+++ b/setup.py
@@ -13,8 +13,6 @@
 
 import numpy as np
 
-__version__ = "1.2.1"
-
 if sys.version_info[0] != 3 or sys.version_info[1] < 6:
     print("Package requires Python version 3.6+")
     sys.exit(1)
@@ -26,49 +24,64 @@
 # see MANIFEST.in -- a slight hack to include all of the header files in the source distrubution
 
 modules = [
-	dict(name="footprint_tools.modeling.predict", sources=["footprint_tools/modeling/predict.pyx"]),
-	dict(name="footprint_tools.modeling.predict", sources=["footprint_tools/modeling/predict.pyx"]),
-	dict(name="footprint_tools.modeling.dispersion", sources=["footprint_tools/modeling/dispersion.pyx"]),
-	dict(name="footprint_tools.stats.windowing", sources = ["footprint_tools/stats/windowing.pyx"]),
-	dict(name="footprint_tools.stats.distributions.nbinom", sources=["footprint_tools/stats/distributions/nbinom.pyx"]),
-	dict(name="footprint_tools.stats.segment", sources=["footprint_tools/stats/segment.pyx"]),
-	dict(name="footprint_tools.stats.fdr.bisect", sources=["footprint_tools/stats/fdr/bisect.pyx"]),
-	dict(name="footprint_tools.stats.differential", sources=["footprint_tools/stats/differential.pyx"])
+    dict(name="footprint_tools.modeling.predict", sources=["footprint_tools/modeling/predict.pyx"]),
+    dict(name="footprint_tools.modeling.dispersion", sources=["footprint_tools/modeling/dispersion.pyx"]),
+    dict(name="footprint_tools.stats.utils", sources=["footprint_tools/stats/utils.pyx"]),
+    dict(name="footprint_tools.stats.windowing", sources = ["footprint_tools/stats/windowing.pyx"]),
+    dict(name="footprint_tools.stats.distributions.nbinom", sources=["footprint_tools/stats/distributions/nbinom.pyx"]),
+    dict(name="footprint_tools.stats.distributions.invchi2", sources=["footprint_tools/stats/distributions/invchi2.pyx"]),
+    dict(name="footprint_tools.stats.distributions.t", sources=["footprint_tools/stats/distributions/t.pyx"]),
+    dict(name="footprint_tools.stats.differential", sources=["footprint_tools/stats/differential.pyx"]),
 ]
 
-scripts = ["scripts/ftd-learn-dispersion-model",
-	"scripts/ftd-compute-deviation",
-	"scripts/ftd-learn-beta-prior",
-	"scripts/ftd-compute-posterior"]
-
-install_requires = ["cython", "numpy>=1.10", "scipy>=0.17", "pandas", "pysam>=0.15", "statsmodels", "genome_tools>=1.0.2", "pwlf", "simplejson", "tqdm"]
+install_requires = [
+    "cython",
+    "numpy>=1.10",
+    "scipy>=0.17",
+    "pandas",
+    "pysam>=0.15",
+    "statsmodels",
+    "genome_tools>=1.0.2",
+    "pwlf",
+    "simplejson",
+    "tqdm",
+    "click",
+    "click-option-group",
+    "matplotlib",
+    "yaspin"
+]
 
 setup(
-	name = "footprint_tools",
-	version = __version__,
-	license = "GPL-3.0-or-later",
-	description = "Genomic footprint detection",
-	long_description=(pathlib.Path(__file__).parent / "README.md").read_text(),
-	long_description_content_type="text/markdown",
-	author = "Jeff Vierstra",
-	author_email = "jvierstra@altius.org",
-	url = "https://github.com/jvierstra/footprint-tools",
-	download_url = "https://github.com/jvierstra/footprint-tools/archive/v{}.tar.gz".format(__version__),
-	keywords = ["genomic footprints", "bioinformatics"],
-	zip_safe = False,
-	packages =  find_packages(),
-	libraries = [cehpes_lib],
+    name = "footprint_tools",
+    version = "1.3.0",
+    license = "GPL-3.0-or-later",
+    description = "Genomic footprint detection",
+    long_description=(pathlib.Path(__file__).parent / "README.md").read_text(),
+    long_description_content_type="text/markdown",
+    author = "Jeff Vierstra",
+    author_email = "jvierstra@altius.org",
+    url = "https://github.com/jvierstra/footprint-tools",
+    download_url = "https://github.com/jvierstra/footprint-tools/archive/v{}.tar.gz".format(__version__),
+    keywords = ["genomic footprints", "bioinformatics", "chromatin", "transcription factors", "dnase"],
+    zip_safe = False,
+    packages =  find_packages(),
+    libraries = [cehpes_lib],
     ext_modules = [Extension(**opts) for opts in modules],
     include_dirs=[np.get_include(), cephes_include],
-	package_data={k:["*.pxd"] for k in find_packages()},
+    package_data={'footprint_tools':['*.pxd', 'logging.conf']},
     cmdclass = {'build_clib': build_clib, 'build_ext': build_ext},
     install_requires = install_requires,
-    scripts = scripts,
+    include_package_data = True,
+    #scripts = scripts,
+    entry_points = {"console_scripts": ["ftd = footprint_tools.__main__:main"]},
     classifiers=[
-	    'Development Status :: 5 - Production/Stable', 
-	    'Intended Audience :: Science/Research', 
-	    'Topic :: Scientific/Engineering :: Bio-Informatics',
-	    'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)',
-	    'Programming Language :: Python :: 3.6',
-],
+        'Development Status :: 5 - Production/Stable', 
+        'Intended Audience :: Science/Research', 
+        'Topic :: Scientific/Engineering :: Bio-Informatics',
+        'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)',
+        'Programming Language :: Python :: 3.6',
+        'Programming Language :: Python :: 3.7',
+        'Programming Language :: Python :: 3.8',
+        'Programming Language :: Python :: 3.9',
+    ],
 )


Test patch:
diff --git a/test/test_predict.py b/test/test_predict.py
deleted file mode 100644
index 2468573..0000000
--- a/test/test_predict.py
+++ /dev/null
@@ -1,34 +0,0 @@
-#!/usr/bin/which python3
-
-
-from genome_tools import genomic_interval, genomic_interval_set
-
-from footprint_tools import cutcounts
-from footprint_tools.modeling import predict, bias
-
-from pysam import FastaFile
-
-
-base_url = "https://resources.altius.org/~jvierstra/projects/footprinting.2020/per.dataset/CD20+-DS18208/"
-bam_file_url = base_url + "reads.bam"
-
-bm_file_path = "data/vierstra_et_al.6mer-model.txt"
-fasta_file_path = "/home/jvierstra/data/genomes/hg38/hg.all.fa"
-
-
-interval = genomic_interval('chr19', 48363826, 48364602)
-
-fasta_fn = FastaFile(fasta_file_path)
-reads_fn = cutcounts.bamfile(bam_file_url)
-
-bm = bias.kmer_model(bm_file_path)
-
-
-print(reads_fn[interval])
-print(fasta_fn.fetch(interval.chrom, interval.start, interval.end))
-
-pred = predict.prediction(reads_fn, fasta_fn, bm, half_window_width = 5, smoothing_half_window_width = 50, smoothing_clip = 0.01)
-
-obs, exp, win = pred.compute(interval)
-
-print(obs['+'].shape, exp['+'].shape, win['+'].shape)