diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000000000000000000000000000000000000..291bf67cb5b3e3cde43061e58d429d61fa0e7e46
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1 @@
+*.bw filter=lfs diff=lfs merge=lfs -text
diff --git a/libhts/__init__.py b/libhts/__init__.py
index 1e4256757ef0d4cbe7b16385035a109fa8193f8a..83eb4951548b23eff68e1f823e2a2efe14695e7e 100644
--- a/libhts/__init__.py
+++ b/libhts/__init__.py
@@ -1,6 +1,6 @@
__copyright__ = "Copyright (C) 2020 Blaise Li"
__licence__ = "GNU GPLv3"
-__version__ = 0.2
+__version__ = 0.3
from .libhts import (
aligner2min_mapq,
gtf_2_genes_exon_lengths,
@@ -8,6 +8,7 @@ from .libhts import (
make_empty_bigwig,
make_seeding_function,
median_ratio_to_pseudo_ref_size_factors,
+ paste_bigwig_region,
plot_boxplots, plot_counts_distribution, plot_histo,
plot_lfc_distribution, plot_MA,
plot_norm_correlations, plot_paired_scatters, plot_scatter,
diff --git a/libhts/data/klp7co_RPF.bw b/libhts/data/klp7co_RPF.bw
new file mode 100644
index 0000000000000000000000000000000000000000..27fd7f77252b39b78b1e4f276e9a996c5bac83ed
--- /dev/null
+++ b/libhts/data/klp7co_RPF.bw
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fcc0e76307e387181256b9e2c3294d93e88f628270649126a55ef93b392f554a
+size 12941493
diff --git a/libhts/libhts.py b/libhts/libhts.py
index 8dc859bb78a49c5ab8cd8d8124cbae85247a6b8c..1cb75107986a2d008e4132c653f3000e9f6262cd 100644
--- a/libhts/libhts.py
+++ b/libhts/libhts.py
@@ -18,6 +18,8 @@ from functools import reduce
import warnings
import numpy as np
import pandas as pd
+# To compute bins in bigwig data
+from scipy.stats import binned_statistic
# To compute correlation coefficient, and compute linear regression
from scipy.stats.stats import pearsonr, linregress
# To compute geometric mean
@@ -249,6 +251,69 @@ def make_empty_bigwig(filename, chrom_sizes):
bw_out.close()
+# Possible improvements: use an itearble of (from_region, to_region) pairs,
+# or two zippable iterables
+def paste_bigwig_region(
+ from_fname, to_fname, from_region, to_region,
+ dest_fname, nanmean=False):
+ """
+
+ Take values from a region *from_region* in a bigwig file *from_fname* and
+ paste them into a region *to_region* in another bigwig file *to_fname*, and
+ write the results in a third bigwig file *dest_fname*.
+
+ A region should be specified as a (chrom, start, stop) triplet, where start
+ is zero-based and stop is 1-based (like BED coordinates or Python slices).
+
+ *dest_fname* will have the same chromosomes as *to_fname*.
+ The values in *to_region* are substituted.
+
+ Current limitation: The regions should have the same length.
+ Option *nanmean* is likely not working correctly.
+ """
+ (from_chrom, from_start, from_stop) = from_region
+ from_len = from_stop - from_start
+ (to_chrom, to_start, to_stop) = to_region
+ assert to_stop - to_start == from_len, (
+ "Regions should have the same lengths.")
+ from_bw = pyBigWig.open(from_fname)
+ to_bw = pyBigWig.open(to_fname)
+ chrom_sizes = list(to_bw.chroms().items())
+ dest_bw = pyBigWig.open(dest_fname, "w")
+ dest_bw.addHeader(chrom_sizes)
+ for (chrom, chrom_len) in chrom_sizes:
+ nb_bins = ceil(chrom_len / 10)
+ # values = to_bw.values(chrom, 0, chrom_len)
+ # Original values, plus some nans for binned_statistics to properly set the bin boundaries
+ values = np.pad(
+ to_bw.values(chrom, 0, chrom_len),
+ # pad zero on the left, and what is needed to complete the bin on the right
+ (0, 10 * nb_bins - chrom_len),
+ constant_values=np.nan)
+ if chrom == to_chrom:
+ # Replace the values in the desired region
+ values[to_start:to_stop] = from_bw.values(
+ from_chrom, from_start, from_stop)
+ if nanmean:
+ bin_means = binned_statistic(
+ range(0, 10 * nb_bins), values,
+ statistic=np.nanmean, bins=nb_bins).statistic
+ else:
+ bin_means = binned_statistic(
+ # range(0, chrom_len), np.nan_to_num(values),
+ range(0, 10 * nb_bins), values,
+ statistic="mean", bins=nb_bins).statistic
+ dest_bw.addEntries(
+ chrom, 0,
+ # Mean for each bin of size 10
+ values=bin_means,
+ # values=np.nan_to_num(np.zeros(chrom_len)[0::10]),
+ span=10, step=10)
+ dest_bw.close()
+ to_bw.close()
+ from_bw.close()
+
+
#################
# Bowtie2 stuff #
#################
diff --git a/libhts/test_libhts.py b/libhts/test_libhts.py
new file mode 100644
index 0000000000000000000000000000000000000000..5cacf63a927e34508784d7dffc5b8974d91a37a2
--- /dev/null
+++ b/libhts/test_libhts.py
@@ -0,0 +1,97 @@
+import os
+import shutil
+import tempfile
+import unittest
+import libhts
+import pyBigWig
+import numpy as np
+
+
+class TestPasteBigWig(unittest.TestCase):
+ def setUp(self):
+ self.bw_in_path = "data/klp7co_RPF.bw"
+ self.from_region = ("klp7co", 0, 2459)
+ self.to_region = ("III", (10808346-20), (10808346-20)+2459)
+ self.tmpdir = tempfile.mkdtemp()
+ self.bw_out_path = os.path.join(
+ self.tmpdir, "pasted.bw")
+ libhts.paste_bigwig_region(
+ self.bw_in_path, self.bw_in_path,
+ self.from_region, self.to_region,
+ self.bw_out_path)
+
+ def tearDown(self):
+ shutil.rmtree(self.tmpdir)
+
+ def test_preserved_region_num_bins(self):
+ bw_in = pyBigWig.open(self.bw_in_path)
+ bw_out = pyBigWig.open(self.bw_out_path)
+ bins_in = {
+ (start, stop): val
+ for (start, stop, val)
+ in bw_in.intervals(*self.from_region)
+ if (not np.isnan(val))}
+ # if (not np.isnan(val)) and val > 0}
+ # contains extra (start, stop, nan),
+ # that's where there was no bin
+ bins_out = {
+ (start, stop): val
+ for (start, stop, val)
+ in bw_out.intervals(*self.from_region)
+ if (not np.isnan(val))}
+ self.assertTrue(len(bins_in) == len(bins_out))
+
+ def test_preserved_region_bin_contents(self):
+ bw_in = pyBigWig.open(self.bw_in_path)
+ bw_out = pyBigWig.open(self.bw_out_path)
+ bins_in = [
+ (start, stop, val)
+ for (start, stop, val)
+ in bw_in.intervals(*self.from_region)
+ if (not np.isnan(val))]
+ bins_out = [
+ (start, stop, val)
+ for (start, stop, val)
+ in bw_out.intervals(*self.from_region)
+ if (not np.isnan(val))]
+ for (i, (bin_in, bin_out)) in enumerate(zip(bins_in, bins_out)):
+ self.assertTrue(bin_in == bin_out)
+
+ # There are differences, maybe due to phase differences with respect to bins
+ # def test_pasted_region_num_bins(self):
+ # bw_in = pyBigWig.open(self.bw_in_path)
+ # bw_out = pyBigWig.open(self.bw_out_path)
+ # bins_in = {
+ # (start, stop): val
+ # for (start, stop, val)
+ # in bw_in.intervals(*self.from_region)
+ # if (not np.isnan(val))}
+ # # if (not np.isnan(val)) and val > 0}
+ # # contains extra (start, stop, nan),
+ # # that's where there was no bin
+ # bins_out = {
+ # (start, stop): val
+ # for (start, stop, val)
+ # in bw_out.intervals(*self.to_region)
+ # if (not np.isnan(val))}
+ # self.assertTrue(len(bins_in) == len(bins_out))
+
+ # def test_pasted_region_bin_contents(self):
+ # bw_in = pyBigWig.open(self.bw_in_path)
+ # bw_out = pyBigWig.open(self.bw_out_path)
+ # bins_in = [
+ # (start, stop, val)
+ # for (start, stop, val)
+ # in bw_in.intervals(*self.from_region)
+ # if (not np.isnan(val))]
+ # bins_out = [
+ # (start, stop, val)
+ # for (start, stop, val)
+ # in bw_out.intervals(*self.to_region)
+ # if (not np.isnan(val))]
+ # for (i, (bin_in, bin_out)) in enumerate(zip(bins_in, bins_out)):
+ # self.assertTrue(bin_in == bin_out)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/scripts/copypaste_bigwig_regions.py b/scripts/copypaste_bigwig_regions.py
new file mode 100755
index 0000000000000000000000000000000000000000..a3d06ca2f6c01206e4c2f9c4e5fa03b1be16db08
--- /dev/null
+++ b/scripts/copypaste_bigwig_regions.py
@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+# Copyright (C) 2020 Blaise Li
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+"""Copy bigwig data from a region in a bigwig file into another region.
+
+The other region can be in the same or in another bigwig file.
+The result will be written in another bigwig file.
+"""
+
+from argparse import (
+ ArgumentParser,
+ ArgumentDefaultsHelpFormatter, RawTextHelpFormatter)
+import sys
+from operator import itemgetter
+from yaml import safe_load as yload
+from libhts import paste_bigwig_region
+
+CONFIG_EXAMPLE = """
+# bigwig file containing the region to copy
+from_bw: "data/klp7co_RPF.bw"
+# bigwig file containing the region to replace
+to_bw: "data/klp7co_RPF.bw"
+# bigwig file in which to write the result
+dest_bw: "/tmp/pasted.bw"
+# region from which to copy (BED-like coordinates)
+from_region:
+ chrom: "klp7co"
+ # zero-based coordinate
+ start: 0
+ # 1-based coordinate
+ end: 2459
+# region in which to copy (BED-like coordinates)
+to_region:
+ chrom: "III"
+ # zero-based coordinate
+ start: 10808326
+ # 1-based coordinate
+ end: 10810785
+"""
+
+
+class Formatter(ArgumentDefaultsHelpFormatter, RawTextHelpFormatter):
+ """
+ Hybrid help formatter.
+ See comments of <https://stackoverflow.com/a/3853776/1878788>
+ """
+
+
+def main():
+ """Run the command-line script."""
+ # parser = ArgumentParser(
+ # description=__doc__,
+ # formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+ parser = ArgumentParser(
+ description=f"{__doc__}\n\nExample configuration:\n{CONFIG_EXAMPLE}",
+ formatter_class=Formatter)
+ parser.add_argument(
+ "-c", "--configfile",
+ help="""YAML-formatted file containing information about:
+ - input bigwig files
+ - output bigwig file
+ - region from which to copy
+ - region where to paste
+""")
+ parser.add_argument(
+ "-e", "--example_config",
+ action="store_true",
+ help="display an example configuration file")
+ args = parser.parse_args()
+ if args.example_config:
+ print(CONFIG_EXAMPLE)
+ return 0
+ if not args.configfile:
+ print("You must provide a configuration file using option -c. "
+ "See option -e for an example.", file=sys.stderr)
+ return 1
+ with open(args.configfile) as config_fh:
+ config = yload(config_fh)
+ region_extractor = itemgetter("chrom", "start", "end")
+ from_region = region_extractor(config["from_region"])
+ to_region = region_extractor(config["to_region"])
+ paste_bigwig_region(
+ config["from_bw"], config["to_bw"],
+ from_region, to_region,
+ config["dest_bw"])
+
+ return 0
+
+
+if __name__ == "__main__":
+ sys.exit(main())
diff --git a/setup.py b/setup.py
index 20425384fe8e3f7966749d4f46fc0e0dcedb059d..fe572f53fb83373fb8957e436859448e4c8a079f 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,9 @@ setup(
author_email="blaise.li@normalesup.org",
license="GNU GPLv3",
packages=find_packages(),
- scripts=["scripts/bam2bigwig.sh", "scripts/sam2indexedbam.sh"],
+ scripts=[
+ "scripts/bam2bigwig.sh", "scripts/sam2indexedbam.sh",
+ "scripts/copypaste_bigwig_regions.py"],
install_requires=[
#"libworkflows @ git+https://gitlab+deploy-token-31:isEzpsgbNf2sJMdUDy2g@gitlab.pasteur.fr/bli/libworkflows.git@744dd79b579577cb6e131653260d7990946be3ad#egg=libworkflows-0.1",
#"libworkflows @ git+https://gitlab+deploy-token-31:isEzpsgbNf2sJMdUDy2g@gitlab.pasteur.fr/bli/libworkflows.git#egg=libworkflows-0.1",