From fa9e08bc2fca11b0bc83bb420405d880d01842f3 Mon Sep 17 00:00:00 2001 From: Victor Alexandru Date: Mon, 20 Dec 2021 17:01:11 -0500 Subject: [PATCH 1/2] feat: capture all reads and gaps on coverage --- src/epivizFileParser/BamFile.py | 82 +++++++++++++++++++++++++++------ 1 file changed, 68 insertions(+), 14 deletions(-) diff --git a/src/epivizFileParser/BamFile.py b/src/epivizFileParser/BamFile.py index 45016a4..f500c0f 100644 --- a/src/epivizFileParser/BamFile.py +++ b/src/epivizFileParser/BamFile.py @@ -76,29 +76,83 @@ def getRange(self, chr, start, end, bins=2000, zoomlvl=-1, metric="AVG", respTyp result = [] try: - iter = self.file.pileup(chr, start, end) + iter = self.file.pileup( + chr, + start, + end, + truncate=True, + stepper='nofilter', + min_base_quality=0, + ) + + chrTemp = valueTemp = None + startTemp = endTemp = start - chrTemp = startTemp = endTemp = valueTemp = None for x in iter: + # get data + _current_reference_name = x.reference_name + _current_reference_pos = x.reference_pos + _current_valueTemp = x.get_num_aligned() # or x.nsegments - if not x.reference_pos >= start or not x.reference_pos <= end: + # print(x.pileups[0]) + + for pileupread in x.pileups: + # or look at x.get_query_sequences() + _current_valueTemp -= pileupread.is_del + + print( + f"{_current_reference_pos} : {_current_valueTemp} : {x.get_query_sequences()}") + + # is first? + if valueTemp is None: + chrTemp = _current_reference_name + startTemp = _current_reference_pos + valueTemp = _current_valueTemp + + endTemp = _current_reference_pos + 1 continue - - if x.reference_pos == end and valueTemp is not None: + + # is it at the end? did it finish before end? + if x.reference_pos == end-1 and valueTemp is not None: + # last min change + if valueTemp is not _current_valueTemp: + result.append((chrTemp, startTemp, endTemp, valueTemp)) + chrTemp = _current_reference_name + startTemp = _current_reference_pos + valueTemp = _current_valueTemp + result.append((chrTemp, startTemp, end, valueTemp)) + continue - if valueTemp is None: - chrTemp = x.reference_name - startTemp = x.reference_pos - valueTemp = x.get_num_aligned() - elif valueTemp is not x.get_num_aligned(): + # gap detection + if _current_reference_pos != endTemp: + print('caught gap') + # print(f"{_current_reference_pos} : {endTemp}") + result.append((chrTemp, startTemp, endTemp, valueTemp)) + chrTemp = _current_reference_name + startTemp = _current_reference_pos + valueTemp = _current_valueTemp + + endTemp = _current_reference_pos + 1 + continue + + # change of coverage value + if valueTemp is not _current_valueTemp: result.append((chrTemp, startTemp, endTemp, valueTemp)) - chrTemp = x.reference_name - startTemp = x.reference_pos - valueTemp = x.get_num_aligned() + chrTemp = _current_reference_name + startTemp = _current_reference_pos + valueTemp = _current_valueTemp + + endTemp = _current_reference_pos + 1 + continue + + # no change + endTemp = _current_reference_pos + 1 - endTemp = x.reference_pos+1 + # last insertion check + if len(result) == 0 or startTemp != result[len(result) - 1][1]: + result.append((chrTemp, startTemp, endTemp, valueTemp)) self.get_col_names() From 5a3ecec7ff0cff77e5a3e2cc76c0129c29deb5c7 Mon Sep 17 00:00:00 2001 From: Victor Alexandru Date: Wed, 22 Dec 2021 12:10:27 -0500 Subject: [PATCH 2/2] fix: junction calculation pileup and tests --- src/epivizFileParser/BamFile.py | 11 +++++------ src/epivizFileParser/SplicingBamFile.py | 25 +++++++++++++++++++++++-- tests/test_BamFile.py | 14 ++++++++++++-- tests/test_SplicingBamFile.py | 20 ++++++-------------- 4 files changed, 46 insertions(+), 24 deletions(-) diff --git a/src/epivizFileParser/BamFile.py b/src/epivizFileParser/BamFile.py index f500c0f..8942325 100644 --- a/src/epivizFileParser/BamFile.py +++ b/src/epivizFileParser/BamFile.py @@ -81,8 +81,6 @@ def getRange(self, chr, start, end, bins=2000, zoomlvl=-1, metric="AVG", respTyp start, end, truncate=True, - stepper='nofilter', - min_base_quality=0, ) chrTemp = valueTemp = None @@ -100,8 +98,8 @@ def getRange(self, chr, start, end, bins=2000, zoomlvl=-1, metric="AVG", respTyp # or look at x.get_query_sequences() _current_valueTemp -= pileupread.is_del - print( - f"{_current_reference_pos} : {_current_valueTemp} : {x.get_query_sequences()}") + # print( + # f"{_current_reference_pos} : {_current_valueTemp} : {x.get_query_sequences()}") # is first? if valueTemp is None: @@ -127,7 +125,7 @@ def getRange(self, chr, start, end, bins=2000, zoomlvl=-1, metric="AVG", respTyp # gap detection if _current_reference_pos != endTemp: - print('caught gap') + # print('caught gap') # print(f"{_current_reference_pos} : {endTemp}") result.append((chrTemp, startTemp, endTemp, valueTemp)) chrTemp = _current_reference_name @@ -152,7 +150,8 @@ def getRange(self, chr, start, end, bins=2000, zoomlvl=-1, metric="AVG", respTyp # last insertion check if len(result) == 0 or startTemp != result[len(result) - 1][1]: - result.append((chrTemp, startTemp, endTemp, valueTemp)) + if valueTemp is not None: + result.append((chrTemp, startTemp, endTemp, valueTemp)) self.get_col_names() diff --git a/src/epivizFileParser/SplicingBamFile.py b/src/epivizFileParser/SplicingBamFile.py index f2e7b49..850a863 100644 --- a/src/epivizFileParser/SplicingBamFile.py +++ b/src/epivizFileParser/SplicingBamFile.py @@ -36,11 +36,32 @@ def get_col_names(self): return self.columns def _get_junctions(self, chr, start, end): - iter = self.file.fetch(chr, start, end) + # iter = self.file.fetch(chr, start, end) + iter = self.file.pileup( + chr, + start, + end, + truncate=True, + ) + + reads_set = set() + reads_refs = [] + + for pileupcolumn in iter: + for pileupread in pileupcolumn.pileups: + _name = pileupread.alignment.query_name + + if _name in reads_set: + pass + else: + reads_set.add(_name) + reads_refs.append(pileupread.alignment) junctions = {} - for read in iter: + for read in reads_refs: + # print(read.query_name) + if read.cigartuples is None: # Skipping read with no CIGAR string continue diff --git a/tests/test_BamFile.py b/tests/test_BamFile.py index 5b0e24b..58f10de 100644 --- a/tests/test_BamFile.py +++ b/tests/test_BamFile.py @@ -3,6 +3,8 @@ import pytest import os +import pandas as pd + from epivizFileParser import BamFile __author__ = "jkanche, alexanv8" @@ -18,11 +20,19 @@ def test_getRange(): - res, error = bb.getRange("chr10", 1, 100000000) + coverage, error = bb.getRange("chr10", 1, 100000000) assert(error == None) + assert(len(bb.columns) == 4) assert(bb.columns == ["chr", "start", "end", "value"]) - assert(len(res) == 8) + assert(len(coverage) == 25) + + assert (coverage['value'].tolist() == [0, 1, 0, 1, 0, 1, 0, + 0, 1, 0, 0, 1, 2, 1, 2, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]) + assert (coverage['start'].tolist() == [10265460, 10265473, 10266155, 10266200, 10266296, 10266311, 10266332, 10266594, 10266646, 10266743, 10266857, + 10267044, 10267063, 10267104, 10267116, 10267117, 10267118, 10267119, 10267121, 10267122, 10267142, 10267224, 10267690, 10267691, 10268169]) + assert (coverage['end'].tolist() == [10265473, 10266155, 10266200, 10266296, 10266311, 10266332, 10266374, 10266646, 10266743, 10266808, 10267044, + 10267063, 10267104, 10267116, 10267117, 10267118, 10267119, 10267121, 10267122, 10267142, 10267224, 10267690, 10267691, 10268169, 10268186]) def test_input_error(): diff --git a/tests/test_SplicingBamFile.py b/tests/test_SplicingBamFile.py index 17e7b4a..78eedf5 100644 --- a/tests/test_SplicingBamFile.py +++ b/tests/test_SplicingBamFile.py @@ -24,21 +24,13 @@ def test_getRange(): assert(len(bb.columns) == 6) assert(bb.columns == ["chr", "region1_start", "region2_start", "region1_end", "region2_end", "value"]) - assert(len(coverage) == 8) - assert(len(junctions) == 14) + assert(len(junctions) == 8) - assert (coverage['value'].tolist() == [0, 1, 0, 1, 0, 1, 2, 1]) - assert (coverage['start'].tolist() == [10265460, 10265473, - 10266332, 10266646, 10266743, 10267044, 10267063, 10267142]) - assert (coverage['end'].tolist() == [10265473, 10266332, - 10266646, 10266743, 10267044, 10267063, 10267142, 10268169]) - - assert (junctions['value'].tolist() == [ - 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) - assert (junctions['region2_start'].tolist() == [10266201, 10266312, 10266657, 10266723, 10266972, - 10266978, 10267047, 10266897, 10267023, 10267117, 10267120, 10267123, 10267225, 10267692]) - assert (junctions['region1_end'].tolist() == [10266155, 10266296, 10266641, 10266707, 10266961, - 10266972, 10267043, 10266881, 10267007, 10267104, 10267117, 10267121, 10267118, 10267690]) + assert (junctions['value'].tolist() == [1, 1, 1, 1, 1, 1, 1, 1]) + assert (junctions['region2_start'].tolist() == [ + 10266201, 10266312, 10267023, 10267117, 10267120, 10267123, 10267225, 10267692]) + assert (junctions['region1_end'].tolist() == [ + 10266155, 10266296, 10267007, 10267104, 10267117, 10267121, 10267118, 10267690]) # same test as bam file (copy paste)