From cb62a73653ba13a98e8facf1dcc97b753a8db5dd Mon Sep 17 00:00:00 2001 From: msbentsen Date: Wed, 7 Aug 2019 09:59:33 +0200 Subject: [PATCH] Bug fix for regions close to chromosome borders --- CHANGES | 3 +++ tobias/__init__.py | 2 +- tobias/footprinting/atacorrect.py | 7 +++++-- tobias/footprinting/atacorrect_functions.py | 2 +- tobias/footprinting/scorebigwig.py | 13 +++++++++---- tobias/utils/regions.py | 3 +++ tobias/utils/utilities.py | 1 + 7 files changed, 23 insertions(+), 8 deletions(-) diff --git a/CHANGES b/CHANGES index e9adb91..04b1599 100644 --- a/CHANGES +++ b/CHANGES @@ -1,3 +1,6 @@ +## 0.6.4 (2019-08-07) +- Bug fix for error with regions very close to chromosome borders for ATACorrect and ScoreBigwig. + ## 0.6.3 (2019-07-16) - Increased size of texts in BINDetect volcano plot and moved label into plot diff --git a/tobias/__init__.py b/tobias/__init__.py index 63af887..364e7ba 100644 --- a/tobias/__init__.py +++ b/tobias/__init__.py @@ -1 +1 @@ -__version__ = "0.6.3" +__version__ = "0.6.4" diff --git a/tobias/footprinting/atacorrect.py b/tobias/footprinting/atacorrect.py index e9059ed..289e023 100644 --- a/tobias/footprinting/atacorrect.py +++ b/tobias/footprinting/atacorrect.py @@ -235,10 +235,13 @@ def run_atacorrect(args): output_regions = RegionList().from_bed(args.regions_out) else: output_regions = deepcopy(peak_regions) - - output_regions.apply_method(OneRegion.extend_reg, args.extend) + + #Extend regions to make sure extend + flanking for window/flank are within boundaries + flank_extend = args.k_flank + int(args.window/2.0) + output_regions.apply_method(OneRegion.extend_reg, args.extend + flank_extend) output_regions.merge() output_regions.apply_method(OneRegion.check_boundary, bam_chrom_info, "cut") + output_regions.apply_method(OneRegion.extend_reg, -flank_extend) #Cut to needed size knowing that the region will be extended in function #Remove blacklisted regions and chromosomes not in common blacklist_regions = RegionList().from_bed(args.blacklist) if args.blacklist != None else RegionList([]) #fill in with regions from args.blacklist diff --git a/tobias/footprinting/atacorrect_functions.py b/tobias/footprinting/atacorrect_functions.py index e9c0d11..908abb2 100644 --- a/tobias/footprinting/atacorrect_functions.py +++ b/tobias/footprinting/atacorrect_functions.py @@ -194,7 +194,7 @@ def bias_correction(regions_list, params, bias_obj): for region_obj in regions_list: region_obj.extend_reg(f_extend) - region_obj.check_boundary(chrom_lengths, "cut") + #region_obj.check_boundary(chrom_lengths, "cut") #moved to outside of function reg_len = region_obj.get_length() #length including flanking reg_key = (region_obj.chrom, region_obj.start+f_extend, region_obj.end-f_extend) #output region out_signals[reg_key] = {"uncorrected":{}, "bias":{}, "expected":{}, "corrected":{}} diff --git a/tobias/footprinting/scorebigwig.py b/tobias/footprinting/scorebigwig.py index f77e443..7d31cab 100644 --- a/tobias/footprinting/scorebigwig.py +++ b/tobias/footprinting/scorebigwig.py @@ -74,6 +74,8 @@ def add_scorebigwig_arguments(parser): #------------------------------------------------------------------# def calculate_scores(regions, args): + logger = TobiasLogger("", args.verbosity, args.log_q) + pybw_signal = pyBigWig.open(args.signal) #cutsites signal pybw_header = pybw_signal.chroms() chrom_lengths = {chrom:int(pybw_header[chrom]) for chrom in pybw_header} @@ -84,6 +86,8 @@ def calculate_scores(regions, args): #Go through each region for i, region in enumerate(regions): + logger.debug("Calculating scores for region: {0}".format(region)) + #Extend region with necessary flank region.extend_reg(flank) reg_key = (region.chrom, region.start+flank, region.end-flank) #output region @@ -168,6 +172,7 @@ def run_scorebigwig(args): pybw_signal = pyBigWig.open(args.signal) pybw_header = pybw_signal.chroms() chrom_info = {chrom:int(pybw_header[chrom]) for chrom in pybw_header} + logger.debug("Chromosome lengths from input bigwig: {0}".format(chrom_info)) #Decide regions logger.info("- Getting output regions ready") @@ -175,7 +180,8 @@ def run_scorebigwig(args): regions = RegionList().from_bed(args.regions) regions.apply_method(OneRegion.extend_reg, args.extend) regions.merge() - regions.apply_method(OneRegion.check_boundary, chrom_info) + regions.apply_method(OneRegion.check_boundary, chrom_info, "cut") + else: regions = RegionList().from_list([OneRegion([chrom, 0, chrom_info[chrom]]) for chrom in chrom_info]) @@ -190,9 +196,8 @@ def run_scorebigwig(args): #Go through each region for i, region in enumerate(regions): region.extend_reg(args.region_flank) - region.check_boundary(chrom_info, "cut") - region.start = region.start + args.region_flank - region.end = region.end - args.region_flank + region = region.check_boundary(chrom_info, "cut") + region.extend_reg(-args.region_flank) #Information for output bigwig reference_chroms = sorted(list(chrom_info.keys())) diff --git a/tobias/utils/regions.py b/tobias/utils/regions.py index 492a069..bdb3053 100644 --- a/tobias/utils/regions.py +++ b/tobias/utils/regions.py @@ -139,6 +139,9 @@ def check_boundary(self, boundaries_dict, action="cut"): else: exit("Error in regions.check_boundary: unknown action") + self[1] = self.start + self[2] = self.end + return(self) def get_signal(self, pybw, numpy_bool = True): diff --git a/tobias/utils/utilities.py b/tobias/utils/utilities.py index 796d4d9..0432665 100644 --- a/tobias/utils/utilities.py +++ b/tobias/utils/utilities.py @@ -168,6 +168,7 @@ def bigwig_writer(q, key_file_dict, header, regions, args): if i_to_write[akey] != no_regions - 1: logger.error("Wrote {0} regions but there are {1} in total".format(i_to_write[akey], len(region_tups))) logger.error("Ready_to_write[{0}]: {1}".format(akey, len(ready_to_write[akey]))) + sys.exit() break #Save key:region:signal to ready_to_write