diff --git a/libhts/__init__.py b/libhts/__init__.py
index c342524162c4ee06ac8f6bd45c2c1a0bf773ade5..5a45ee3d20d74e0d790407fe56946f55080bea73 100644
--- a/libhts/__init__.py
+++ b/libhts/__init__.py
@@ -1,5 +1,7 @@
 from .libhts import (
-    do_deseq2, gtf_2_genes_exon_lengths, median_ratio_to_pseudo_ref_size_factors,
+    do_deseq2, gtf_2_genes_exon_lengths,
+    make_empty_bigwig,
+    median_ratio_to_pseudo_ref_size_factors,
     plot_boxplots, plot_counts_distribution, plot_histo,
     plot_lfc_distribution, plot_MA,
     plot_norm_correlations, plot_paired_scatters, plot_scatter,
diff --git a/libhts/libhts.py b/libhts/libhts.py
index 2597cf331387ceb6b91f2715732dad301425921e..b3775766fea5c565be3b82cc8a08cbb27e24f24f 100644
--- a/libhts/libhts.py
+++ b/libhts/libhts.py
@@ -39,6 +39,7 @@ from rpy2.rinterface import RRuntimeError
 from rpy2.robjects.packages import importr
 deseq2 = importr("DESeq2")
 from pybedtools import BedTool
+import pyBigWig
 import networkx as nx
 
 
@@ -175,6 +176,14 @@ def spikein_gtf_2_lengths(spikein_gtf):
         name=("union_exon_len")).rename_axis("gene"))
 
 
+def make_empty_bigwig(filename, chrom_sizes):
+    bw_out = pyBigWig.open(filename, "w")
+    bw_out.addHeader(list(chrom_sizes.items()))
+    for (chrom, chrom_len) in bw_out.chroms().items():
+        bw_out.addEntries(chrom, 0, values=np.nan_to_num(np.zeros(chrom_len)[0::10]), span=10, step=10)
+    bw_out.close()
+
+
 def do_deseq2(cond_names, conditions, counts_data,
               formula=None, contrast=None, deseq2_args=None):
     """Runs a DESeq2 differential expression analysis."""