Add indel storing and fix cleanup

Author: ThijsMaas

Diff for: iss/app.py

@@ -131,7 +131,7 @@ def generate_reads(args):
                 args.output + ".vcf",
                 "##fileformat=VCFv4.1\n" + "\t".join(["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"]),
             )
-        full_tmp_list = temp_R1 + temp_R2
+        full_tmp_list = temp_R1 + temp_R2 + temp_mut
         full_tmp_list.append(genome_file)
         if os.path.exists("%s.memmap" % args.output):
             full_tmp_list.append("%s.memmap" % args.output)

Diff for: iss/error_models/__init__.py

@@ -171,7 +171,7 @@ def introduce_indels(self, record, orientation, full_seq, bounds):
             bounds (tuple): the position of the read in the full_sequence
 
         Returns:
-            Seq: a sequence with (eventually) indels
+            SeqRecord: a sequence record with indel errors
         """
 
         # get the right indel arrays
@@ -194,11 +194,35 @@ def introduce_indels(self, record, orientation, full_seq, bounds):
                     if random.random() < prob:
                         # we want to insert after the base read
                         mutable_seq.insert(position + 1, str(nucl_to_insert))
+                        if self.store_mutations:
+                            record.annotations["mutations"].append(
+                                {
+                                    "id": record.id,
+                                    "position": position,
+                                    "ref": mutable_seq[position],
+                                    "alt": mutable_seq[position] + nucl_to_insert,
+                                    "quality": ".",
+                                    "type": "ins",
+                                }
+                            )
+
                 if random.random() < deletions[position][mutable_seq[nucl].upper()]:
                     mutable_seq.pop(position)
+                    if self.store_mutations:
+                        record.annotations["mutations"].append(
+                            {
+                                "id": record.id,
+                                "position": position,
+                                "ref": mutable_seq[position],
+                                "alt": ".",
+                                "quality": ".",
+                                "type": "del",
+                            }
+                        )
                 position += 1
             except IndexError:
                 continue
 
         seq = self.adjust_seq_length(mutable_seq, orientation, full_seq, bounds)
-        return seq
+        record.seq = seq
+        return record

Diff for: iss/generator.py

@@ -153,7 +153,7 @@ def simulate_read(record, error_model, i, cpu_number, sequence_type):
     forward.annotations["original"] = str(forward.seq)
 
     # add the indels, the qual scores and modify the record accordingly
-    forward.seq = error_model.introduce_indels(forward, "forward", sequence, bounds)
+    forward = error_model.introduce_indels(forward, "forward", sequence, bounds)
     forward = error_model.introduce_error_scores(forward, "forward")
     forward = error_model.mut_sequence(forward, "forward")
 
@@ -185,7 +185,7 @@ def simulate_read(record, error_model, i, cpu_number, sequence_type):
     reverse.annotations["original"] = str(reverse.seq)
 
     # add the indels, the qual scores and modify the record accordingly
-    reverse.seq = error_model.introduce_indels(reverse, "reverse", sequence, bounds)
+    reverse = error_model.introduce_indels(reverse, "reverse", sequence, bounds)
     reverse = error_model.introduce_error_scores(reverse, "reverse")
     reverse = error_model.mut_sequence(reverse, "reverse")
 

Diff for: iss/test/test_error_model.py

@@ -66,7 +66,7 @@ def test_introduce_indels():
     bounds = (5, 130)
     read = SeqRecord(Seq(str("ATATA" * 25)), id="read_1", description="test read")
     ref_genome = SeqRecord(Seq(str("ATATA" * 100)), id="ref_genome", description="test reference")
-    read.seq = err_mod.introduce_indels(read, "forward", ref_genome, bounds)
+    read = err_mod.introduce_indels(read, "forward", ref_genome, bounds)
     assert len(read.seq) == 125
     assert read.seq[:10] == "ATGATAATAT"
 
@@ -81,7 +81,7 @@ def test_adjust_seq_length_extend():
     bounds = (480, 500)
     read = SeqRecord(Seq(str("ATTTA" * 4)), id="read_1", description="test read")
     ref_genome = SeqRecord(Seq(str("ATTTA" * 100)), id="ref_genome", description="test reference")
-    read.seq = err_mod.introduce_indels(read, "forward", ref_genome, bounds)
+    read = err_mod.introduce_indels(read, "forward", ref_genome, bounds)
     assert len(read.seq) == 20
     assert read.seq[:10] == "TTAATTTAAT"
     assert read.seq[10:] == "TTAATTTAAA"
@@ -100,7 +100,7 @@ def test_introduce_indels_rev():
 
     ref_genome = SeqRecord(Seq("GG" + str("GTACC" * 100) + "GG"), id="ref_genome", description="test reference")
     read = SeqRecord(Seq(rev_comp(str(ref_genome.seq[484:504]))), id="read_1", description="test read")
-    read.seq = err_mod.introduce_indels(read, "reverse", ref_genome, bounds)
+    read = err_mod.introduce_indels(read, "reverse", ref_genome, bounds)
     assert len(read.seq) == 20
     assert read.seq == "CGTACGGTACGGTACGGTAC"