refactor: compare alignments between mapped reads only

zavolanlab · Dec 13, 2023 · 04fcab2 · 04fcab2
1 parent 930b741
commit 04fcab2
Showing 1 changed file with 36 additions and 25 deletions.
diff --git 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()