From 0f548e631cdc2a539e9e0a862a0bdf9ea573c9e5 Mon Sep 17 00:00:00 2001 From: IanSudbery Date: Sat, 16 May 2020 22:20:55 +0100 Subject: [PATCH 1/3] Fixes to BamTools.geneprofile.RangerCounterBigWig to work with pyBigWig --- cgat/BamTools/geneprofile.pyx | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/cgat/BamTools/geneprofile.pyx b/cgat/BamTools/geneprofile.pyx index 670362e45..92ce2320a 100644 --- a/cgat/BamTools/geneprofile.pyx +++ b/cgat/BamTools/geneprofile.pyx @@ -470,17 +470,24 @@ class RangeCounterBigWig(RangeCounter): for start, end in ranges: length = end - start + if contig not in wigfile.chroms(): + continue + values = wigfile.values(contig, start, end) if values == None: continue - for tstart, tend, value in values: - rstart = tstart - start + current_offset - rend = tend - start + current_offset - for i from rstart <= i < rend: counts[i] = value + for tstart, value in enumerate(values): + rstart = tstart + current_offset + + # check for nan values + if value != value: + value = 0.0 + counts[rstart] = value current_offset += length - + + class IntervalsCounter: '''aggregate reads/intervals within a certain arrangement of intervals (sections). Sections are defined by a transcript-model From 445b45354aa2415cae2779fdebc25bb235450536 Mon Sep 17 00:00:00 2001 From: IanSudbery Date: Sat, 16 May 2020 22:22:48 +0100 Subject: [PATCH 2/3] Tidying up extraneous cdefs in BamTools.geneprofile.RangeCounterBigWid --- cgat/BamTools/geneprofile.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cgat/BamTools/geneprofile.pyx b/cgat/BamTools/geneprofile.pyx index 92ce2320a..b52e3882b 100644 --- a/cgat/BamTools/geneprofile.pyx +++ b/cgat/BamTools/geneprofile.pyx @@ -455,8 +455,7 @@ class RangeCounterBigWig(RangeCounter): def count(self, counts, files, contig, ranges ): # collect pileup profile in region bounded by start and end. - cdef int i - cdef int rstart, rend, start, end, tstart, tend + cdef int rstart, start, end, tstart if len(ranges) == 0: return @@ -483,6 +482,7 @@ class RangeCounterBigWig(RangeCounter): # check for nan values if value != value: value = 0.0 + counts[rstart] = value current_offset += length From fe620f2c9a2e4406078c444f952d03a675789d82 Mon Sep 17 00:00:00 2001 From: IanSudbery Date: Sun, 17 May 2020 14:41:31 +0100 Subject: [PATCH 3/3] bam2geneprofile BigWig Counter - deal with zero length intervals, and intervals that go off the end of the contig. Currently silently skipped. --- cgat/BamTools/geneprofile.pyx | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/cgat/BamTools/geneprofile.pyx b/cgat/BamTools/geneprofile.pyx index b52e3882b..41432ec02 100644 --- a/cgat/BamTools/geneprofile.pyx +++ b/cgat/BamTools/geneprofile.pyx @@ -453,7 +453,7 @@ class RangeCounterBigWig(RangeCounter): "getTotal not implemented for RangeCounterBigwig") def count(self, counts, files, contig, ranges ): - + # collect pileup profile in region bounded by start and end. cdef int rstart, start, end, tstart @@ -465,14 +465,28 @@ class RangeCounterBigWig(RangeCounter): for wigfile in files: + if contig not in wigfile.chroms(): + continue + + contig_length = wigfile.chroms(contig) + current_offset = 0 for start, end in ranges: + + start = min(start, contig_length) + end = min(end, contig_length) length = end - start - if contig not in wigfile.chroms(): + if length <= 0: continue - - values = wigfile.values(contig, start, end) + + try: + values = wigfile.values(contig, start, end) + except RuntimeError: + E.error("Invalid interval found: %s\t%i\t%i contig length: %i" % + (contig, start, end, wigfile.chroms(contig))) + raise + if values == None: continue