feat: fingerprint --extendReads uses fragsize from respective bamfile (order independence) (#1483)

Author: WardDeb

pydeeptools/deeptools/countReadsPerBin.py

@@ -1,4 +1,5 @@
 import contextlib
+import functools
 import multiprocessing
 import os
 import shutil
@@ -198,35 +199,17 @@ def __init__(self, bamFilesList, binLength=50, numberOfSamples=None, numberOfPro
         self.bed_and_bin = bed_and_bin
         self.genomeChunkSize = genomeChunkSize
 
-        if extendReads and len(bamFilesList):
-            from deeptools.getFragmentAndReadSize import get_read_and_fragment_length
-            frag_len_dict, read_len_dict = get_read_and_fragment_length(bamFilesList[0],
-                                                                        return_lengths=False,
-                                                                        blackListFileName=blackListFileName,
-                                                                        numberOfProcessors=numberOfProcessors,
-                                                                        verbose=verbose)
-            if extendReads is True:
-                # try to guess fragment length if the bam file contains paired end reads
-                if frag_len_dict:
-                    self.defaultFragmentLength = int(frag_len_dict['median'])
-                else:
-                    sys.exit("*ERROR*: library is not paired-end. Please provide an extension length.")
-                if verbose:
-                    print("Fragment length based on paired en data "
-                          "estimated to be {}".format(frag_len_dict['median']))
-
-            elif extendReads < read_len_dict['median']:
-                sys.stderr.write("*WARNING*: read extension is smaller than read length (read length = {}). "
-                                 "Reads will not be extended.\n".format(int(read_len_dict['median'])))
-                self.defaultFragmentLength = 'read length'
-
-            elif extendReads > 2000:
-                sys.exit(f"*ERROR*: read extension must be smaller that 2000. Value give: {extendReads} ")
-            else:
-                self.defaultFragmentLength = int(extendReads)
-
+        if extendReads:
+            self.defaultFragmentLengths = [self.get_default_fragment_length(fname, extendReads,
+                                                                            blackListFileName=blackListFileName,
+                                                                            numberOfProcessors=numberOfProcessors,
+                                                                            verbose=verbose)
+                                           for fname in bamFilesList]
         else:
-            self.defaultFragmentLength = 'read length'
+            self.defaultFragmentLengths = ['read length'] * len(bamFilesList)
+
+        # direct call fallback
+        self.defaultFragmentLength = self.defaultFragmentLengths[0] if self.defaultFragmentLengths else 'read length'
 
         self.numberOfProcessors = numberOfProcessors
         self.verbose = verbose
@@ -255,12 +238,8 @@ def __init__(self, bamFilesList, binLength=50, numberOfSamples=None, numberOfPro
         if numberOfSamples is None and stepSize is None and bedFile is None:
             raise ValueError("either stepSize, numberOfSamples or bedFile have to be set")
 
-        if self.defaultFragmentLength != 'read length':
-            self.maxPairedFragmentLength = 4 * self.defaultFragmentLength
-        else:
-            self.maxPairedFragmentLength = 1000
-        if self.maxFragmentLength > 0:
-            self.maxPairedFragmentLength = self.maxFragmentLength
+        self.maxPairedFragmentLengths = [self.get_max_paired_fragment_length(x) for x in self.defaultFragmentLengths]
+        self.maxPairedFragmentLength = self.get_max_paired_fragment_length(self.defaultFragmentLength)
 
         if len(self.mappedList) == 0:
             try:
@@ -273,6 +252,50 @@ def __init__(self, bamFilesList, binLength=50, numberOfSamples=None, numberOfPro
                 self.mappedList = []
                 self.statsList = []
 
+    @staticmethod
+    def get_default_fragment_length(bamFile, extendReads, blackListFileName=None,
+                                    numberOfProcessors=1, verbose=False):
+        """
+        Returns the length to which reads of the given file are extended,
+        or 'read length' if reads are not to be extended.
+        If extendReads is True, the median fragment length of the file is used.
+        """
+        from deeptools.getFragmentAndReadSize import get_read_and_fragment_length
+        frag_len_dict, read_len_dict = get_read_and_fragment_length(bamFile,
+                                                                    return_lengths=False,
+                                                                    blackListFileName=blackListFileName,
+                                                                    numberOfProcessors=numberOfProcessors,
+                                                                    verbose=verbose)
+        if extendReads is True:
+            # try to guess fragment length if the bam file contains paired end reads
+            if not frag_len_dict:
+                sys.exit(f"*ERROR*: library {bamFile} is not paired-end. Please provide an extension length.")
+            if verbose:
+                print(f"Fragment length of {bamFile} based on paired en data "
+                      f"estimated to be {frag_len_dict['median']}")
+            return int(frag_len_dict['median'])
+
+        if extendReads < read_len_dict['median']:
+            sys.stderr.write(f"*WARNING*: read extension is smaller than read length (read length = {int(read_len_dict['median'])}) "
+                             f"for {bamFile}. Reads will not be extended.\n")
+            return 'read length'
+
+        if extendReads > 2000:
+            sys.exit(f"*ERROR*: read extension must be smaller that 2000. Value give: {extendReads} ")
+
+        return int(extendReads)
+
+    def get_max_paired_fragment_length(self, defaultFragmentLength):
+        """
+        Returns the maximum fragment length for which read pairs are
+        considered proper pairs, given the default fragment length of a file.
+        """
+        if self.maxFragmentLength > 0:
+            return self.maxFragmentLength
+        if defaultFragmentLength != 'read length':
+            return 4 * defaultFragmentLength
+        return 1000
+
     def get_chunk_length(self, bamFilesHandles, genomeSize, chromSizes, chrLengths):
         # Try to determine an optimal fraction of the genome (chunkSize) that is sent to
         # workers for analysis. If too short, too much time is spent loading the files
@@ -499,9 +522,11 @@ def count_reads_in_region(self, chrom, start, end, bed_regions_list=None):
             else:
                 _file_name = ''
 
-            for bam in bam_handles:
+            for bam_idx, bam in enumerate(bam_handles):
                 for trans in transcriptsToConsider:
-                    tcov = self.get_coverage_of_region(bam, chrom, trans)
+                    tcov = self.get_coverage_of_region(bam, chrom, trans,
+                                                       defaultFragmentLength=self.defaultFragmentLengths[bam_idx],
+                                                       maxPairedFragmentLength=self.maxPairedFragmentLengths[bam_idx])
                     if bed_regions_list is not None and not self.bed_and_bin:
                         subnum_reads_per_bin.append(np.sum(tcov))
                     else:
@@ -537,11 +562,17 @@ def count_reads_in_region(self, chrom, start, end, bed_regions_list=None):
         return subnum_reads_per_bin, _file_name
 
     def get_coverage_of_region(self, bamHandle, chrom, regions,
-                               fragmentFromRead_func=None):
+                               fragmentFromRead_func=None,
+                               defaultFragmentLength=None,
+                               maxPairedFragmentLength=None):
         """
         Returns a numpy array that corresponds to the number of reads
         that overlap with each tile.
 
+        defaultFragmentLength and maxPairedFragmentLength are the values for
+        the file of bamHandle. If not given, self.defaultFragmentLength and
+        self.maxPairedFragmentLength are used.
+
         >>> test = Tester()
         >>> import pysam
         >>> c = CountReadsPerBin([], stepSize=1, extendReads=300)
@@ -572,8 +603,14 @@ def get_coverage_of_region(self, bamHandle, chrom, regions,
 
 
         """
+        if defaultFragmentLength is None:
+            defaultFragmentLength = self.defaultFragmentLength
+        if maxPairedFragmentLength is None:
+            maxPairedFragmentLength = self.maxPairedFragmentLength
         if not fragmentFromRead_func:
-            fragmentFromRead_func = self.get_fragment_from_read
+            fragmentFromRead_func = functools.partial(self.get_fragment_from_read,
+                                                      defaultFragmentLength=defaultFragmentLength,
+                                                      maxPairedFragmentLength=maxPairedFragmentLength)
         nbins = len(regions)
         if len(regions[0]) == 3:
             nbins = 0
@@ -583,10 +620,10 @@ def get_coverage_of_region(self, bamHandle, chrom, regions,
                     nbins += 1
         coverages = np.zeros(nbins, dtype='float64')
 
-        if self.defaultFragmentLength == 'read length':
+        if defaultFragmentLength == 'read length':
             extension = 0
         else:
-            extension = self.maxPairedFragmentLength
+            extension = maxPairedFragmentLength
 
         blackList = None
         if self.blackListFileName is not None:
@@ -779,7 +816,7 @@ def is_proper_pair(read, maxPairedFragmentLength):
                 return True
         return False
 
-    def get_fragment_from_read(self, read):
+    def get_fragment_from_read(self, read, defaultFragmentLength=None, maxPairedFragmentLength=None):
         """Get read start and end position of a read.
         If given, the reads are extended as follows:
         If reads are paired end, each read mate is extended to match
@@ -823,6 +860,14 @@ def get_fragment_from_read(self, read):
         ----------
         read : pysam read object
 
+        defaultFragmentLength : int or 'read length'
+            The fragment length used to extend reads of the file of `read`.
+            If None, self.defaultFragmentLength is used.
+
+        maxPairedFragmentLength : int
+            Pairs with a larger fragment length are extended as single reads.
+            If None, self.maxPairedFragmentLength is used.
+
 
         Returns
         -------
@@ -865,11 +910,16 @@ def get_fragment_from_read(self, read):
         # E.g for a cigar of 40M260N22M
         # get blocks return two elements for the first 40 matches
         # and the for the last 22 matches.
-        if self.defaultFragmentLength == 'read length':
+        if defaultFragmentLength is None:
+            defaultFragmentLength = self.defaultFragmentLength
+        if maxPairedFragmentLength is None:
+            maxPairedFragmentLength = self.maxPairedFragmentLength
+
+        if defaultFragmentLength == 'read length':
             return read.get_blocks()
 
         else:
-            if self.is_proper_pair(read, self.maxPairedFragmentLength):
+            if self.is_proper_pair(read, maxPairedFragmentLength):
                 if read.is_reverse:
                     fragmentStart = read.next_reference_start
                     fragmentEnd = read.reference_end
@@ -882,11 +932,11 @@ def get_fragment_from_read(self, read):
             # Extend using the default fragment length
             else:
                 if read.is_reverse:
-                    fragmentStart = read.reference_end - self.defaultFragmentLength
+                    fragmentStart = read.reference_end - defaultFragmentLength
                     fragmentEnd = read.reference_end
                 else:
                     fragmentStart = read.reference_start
-                    fragmentEnd = read.reference_start + self.defaultFragmentLength
+                    fragmentEnd = read.reference_start + defaultFragmentLength
 
         if self.center_read:
             fragmentCenter = fragmentEnd - (fragmentEnd - fragmentStart) / 2

pydeeptools/deeptools/sumCoveragePerBin.py

@@ -1,3 +1,4 @@
+import functools
 import multiprocessing
 import time
 
@@ -13,11 +14,17 @@ class SumCoveragePerBin(countReadsPerBin.CountReadsPerBin):
     There, we need to sum the per-base coverage.
     """
     def get_coverage_of_region(self, bamHandle, chrom, regions,
-                               fragmentFromRead_func=None):
+                               fragmentFromRead_func=None,
+                               defaultFragmentLength=None,
+                               maxPairedFragmentLength=None):
         """
         Returns a numpy array that corresponds to the number of reads
         that overlap with each tile.
 
+        defaultFragmentLength and maxPairedFragmentLength are the values for
+        the file of bamHandle. If not given, self.defaultFragmentLength and
+        self.maxPairedFragmentLength are used.
+
         >>> test = Tester()
         >>> import pysam
         >>> c = SumCoveragePerBin([], stepSize=1, extendReads=300)
@@ -38,19 +45,25 @@ def get_coverage_of_region(self, bamHandle, chrom, regions,
 
 
         """
+        if defaultFragmentLength is None:
+            defaultFragmentLength = self.defaultFragmentLength
+        if maxPairedFragmentLength is None:
+            maxPairedFragmentLength = self.maxPairedFragmentLength
         if not fragmentFromRead_func:
-            fragmentFromRead_func = self.get_fragment_from_read
+            fragmentFromRead_func = functools.partial(self.get_fragment_from_read,
+                                                      defaultFragmentLength=defaultFragmentLength,
+                                                      maxPairedFragmentLength=maxPairedFragmentLength)
         nbins = len(regions)
         if len(regions[0]) == 3:
             nbins = 0
             for reg in regions:
                 nbins += (reg[1] - reg[0]) // reg[2]
         coverages = np.zeros(nbins, dtype='float64')
 
-        if self.defaultFragmentLength == 'read length':
+        if defaultFragmentLength == 'read length':
             extension = 0
         else:
-            extension = self.maxPairedFragmentLength
+            extension = maxPairedFragmentLength
 
         blackList = None
         if self.blackListFileName is not None:

pydeeptools/deeptools/test/test_plotFingerprint.py

@@ -1,13 +1,13 @@
 import os
 from tempfile import NamedTemporaryFile
 
+import numpy as np
 from matplotlib.testing.compare import compare_images
 
 import deeptools.plotFingerprint
 
 TEST_DATA = os.path.dirname(os.path.abspath(__file__)) + "/test_data/"
 ROOT = os.path.dirname(os.path.abspath(__file__)) + "/test_plotFingerprint/"
-
 tolerance = 13
 
 
@@ -135,4 +135,46 @@ def test_plotFingerprint_quality_metrics_and_JSD():
             assert abs(float(rows["test2"][jsd]) - 0.078613413909822) < 1e-4
 
         finally:
-            cleanup(plotfile.name, qcfile.name)
+            cleanup(plotfile.name, qcfile.name)
+
+
+def run_plotFingerprint_extendReads(bamfiles, labels, rawfile, qcfile):
+    args = (
+        ["-b"] + [TEST_DATA + x for x in bamfiles]
+        + ["-l"] + labels
+        + ["--extendReads",
+           "--region", "chr2:4999000:5003000", "--binSize", "10", "--numberOfSamples", "400", "-p", "1",
+           "--outRawCounts", rawfile, "--outQualityMetrics", qcfile]
+    )
+    deeptools.plotFingerprint.main(args)
+
+    counts = np.loadtxt(rawfile, skiprows=2, ndmin=2)
+    raw = {label: counts[:, idx] for idx, label in enumerate(labels)}
+    with open(qcfile) as _foo:
+        lines = [line.rstrip("\n").split("\t") for line in _foo]
+    qc = {row[0]: np.array(row[1:], dtype=float) for row in lines[1:]}
+
+    return raw, qc
+
+
+def test_plotFingerprint_extendReads_single_bam(tmp_path):
+    raw, qc = run_plotFingerprint_extendReads(["test_paired2.bam"], ["test_paired2"],
+                                              str(tmp_path / "raw.txt"), str(tmp_path / "qc.txt"))
+    expected_counts = np.loadtxt(ROOT + "test_plotFingerprint_extendReads_raw.txt", skiprows=2)
+    np.testing.assert_array_equal(raw["test_paired2"], expected_counts)
+    expected_qc = [0.10334787097692963, 0.4761528725688839, 0.595, 2.0396151481800195e-14,
+                   0.745, 0.5373923764988261, 0.324362640621008]
+    np.testing.assert_allclose(qc["test_paired2"], expected_qc, rtol=1e-6, atol=1e-12)
+
+
+def test_plotFingerprint_extendReads_order_independent(tmp_path):
+    bamfiles = ["test_paired.bam", "test_paired2.bam", "test_paired2.cram"]
+    labels = ["paired", "paired2", "paired2_cram"]
+    raw1, qc1 = run_plotFingerprint_extendReads(bamfiles, labels,
+                                                str(tmp_path / "raw1.txt"), str(tmp_path / "qc1.txt"))
+    raw2, qc2 = run_plotFingerprint_extendReads(bamfiles[::-1], labels[::-1],
+                                                str(tmp_path / "raw2.txt"), str(tmp_path / "qc2.txt"))
+
+    for label in labels:
+        np.testing.assert_array_equal(raw1[label], raw2[label])
+        np.testing.assert_array_equal(qc1[label], qc2[label])