diff --git a/htsinfer/get_library_type.py b/htsinfer/get_library_type.py index 4d4257c..545c4d2 100644 --- a/htsinfer/get_library_type.py +++ b/htsinfer/get_library_type.py @@ -144,6 +144,8 @@ def _align_mates(self): samfile1 = pysam.AlignmentFile(str(alignment_1), 'r') samfile2 = pysam.AlignmentFile(str(alignment_2), 'r') + seq_id1 = None + seq_id2 = None previous_seq_id1 = None previous_seq_id2 = None @@ -154,43 +156,52 @@ def _align_mates(self): concordant = 0 for read1 in samfile1: - seq_id1 = read1.query_name - if seq_id1 != previous_seq_id1 \ - and previous_seq_id1 is not None: - mate1.append(reads1.copy()) - reads1.clear() - if read1.reference_end: - reads1.append(read1) + if not read1.flag & (1 << 2) and isinstance(read1.query_name, str): + seq_id1 = read1.query_name + if ( + seq_id1 != previous_seq_id1 and + previous_seq_id1 is not None + ): + mate1.append(reads1.copy()) + reads1.clear() + if read1.reference_end: + reads1.append(read1) previous_seq_id1 = seq_id1 mate1.append(reads1.copy()) read_counter = 0 for read2 in samfile2: - seq_id2 = read2.query_name - if seq_id2 != previous_seq_id2 \ - and previous_seq_id2 is not None: - if self._compare_alignments(mate1[read_counter], reads2): - concordant += 1 - reads2.clear() - read_counter += 1 - if read2.reference_end: - reads2.append(read2) + if not read2.flag & (1 << 2) and isinstance(read2.query_name, str): + seq_id2 = read2.query_name + if seq_id2 != previous_seq_id2 \ + and previous_seq_id2 is not None: + if self._compare_alignments(mate1[read_counter], reads2): + concordant += 1 + reads2.clear() + read_counter += 1 + if read2.reference_end: + reads2.append(read2) previous_seq_id2 = seq_id2 if self._compare_alignments(mate1[read_counter], reads2): concordant += 1 - if (concordant / read_counter) >= self.cutoff: - self.results.relationship = ( - StatesTypeRelationship.split_mates - ) - self.mapping.library_type.relationship \ - = StatesTypeRelationship.split_mates - self.mapping.mapped = False - self.mapping.star_dirs = [] + if read_counter > 0: + if (concordant / read_counter) >= self.cutoff: + self.results.relationship = ( + StatesTypeRelationship.split_mates + ) + self.mapping.library_type.relationship \ + = StatesTypeRelationship.split_mates + self.mapping.mapped = False + self.mapping.star_dirs = [] + else: + self.results.relationship = ( + StatesTypeRelationship.not_mates + ) else: self.results.relationship = ( - StatesTypeRelationship.not_mates + StatesTypeRelationship.not_available ) samfile1.close()