diff --git a/pydeeptools/deeptools/countReadsPerBin.py b/pydeeptools/deeptools/countReadsPerBin.py index 29ba785cd5..85c2a42c89 100644 --- a/pydeeptools/deeptools/countReadsPerBin.py +++ b/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 diff --git a/pydeeptools/deeptools/sumCoveragePerBin.py b/pydeeptools/deeptools/sumCoveragePerBin.py index 75350c0851..aea4b29bd3 100644 --- a/pydeeptools/deeptools/sumCoveragePerBin.py +++ b/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,8 +45,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 @@ -47,10 +60,10 @@ def get_coverage_of_region(self, bamHandle, chrom, 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: diff --git a/pydeeptools/deeptools/test/test_plotFingerprint.py b/pydeeptools/deeptools/test/test_plotFingerprint.py index 1292d1a6c8..664fbcb280 100644 --- a/pydeeptools/deeptools/test/test_plotFingerprint.py +++ b/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) \ No newline at end of file + 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]) diff --git a/pydeeptools/deeptools/test/test_plotFingerprint/test_plotFingerprint_extendReads_raw.txt b/pydeeptools/deeptools/test/test_plotFingerprint/test_plotFingerprint_extendReads_raw.txt new file mode 100644 index 0000000000..a0fd648ebd --- /dev/null +++ b/pydeeptools/deeptools/test/test_plotFingerprint/test_plotFingerprint_extendReads_raw.txt @@ -0,0 +1,402 @@ +#plotFingerprint --outRawCounts +'pydeeptools/deeptools/test/test_data/test_paired2.bam' +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +20 +20 +20 +20 +20 +20 +20 +20 +20 +20 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +14 +20 +20 +20 +20 +20 +20 +20 +44 +60 +60 +60 +60 +60 +60 +60 +60 +60 +60 +60 +60 +60 +60 +60 +20 +22 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +20 +28 +40 +40 +40 +40 +40 +40 +40 +42 +114 +132 +140 +150 +166 +180 +180 +180 +180 +180 +200 +226 +240 +220 +220 +220 +220 +220 +220 +240 +250 +220 +220 +240 +240 +240 +240 +240 +240 +272 +300 +300 +300 +300 +240 +240 +200 +200 +200 +200 +184 +170 +170 +150 +150 +110 +90 +90 +90 +90 +90 +90 +70 +70 +70 +70 +70 +71 +80 +80 +60 +60 +20 +20 +20 +20 +20 +29 +30 +30 +30 +30 +30 +30 +30 +30 +30 +30 +30 +20 +20 +20 +20 +20 +20 +20 +20 +20 +20 +20 +20 +20 +10 +10 +10 +10 +10 +10 +10 +10 +10 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/pyproject.toml b/pyproject.toml index ded513f661..56d35d64d7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,7 +33,7 @@ authors = [ {name="Andreas S Richter"}, {name="Devon P Ryan"}, {name="Pavankumar Videm"}, - {email="bioinfo-core@ie-freiburg.mpg.de"} + {email="w@rddeboutte.com"} ] requires-python = ">= 3.12" dependencies = [