Convert coverage file to read_count dict

iss/abundance.py

@@ -315,3 +315,20 @@ def expand_draft_abundance(abundance_dic, draft, mode="abundance"):
                     elif mode == "coverage":
                         draft_dic[record.id] = abundance
     return draft_dic
+
+
+def to_readcount(coverage_dic, read_length, genome_file):
+    logger = logging.getLogger(__name__)
+    readcount_dic = {}
+    f = open(genome_file, "r")  # re-opens the file
+    with f:
+        fasta_file = SeqIO.parse(f, "fasta")
+        for record in fasta_file:
+            try:
+                record_coverage = coverage_dic[record.id]
+            except KeyError as e:
+                logger.warning("Fasta record not found in coverage file: %s" % e)
+                continue
+            n_reads = int((record_coverage * len(record.seq)) / read_length)
+            readcount_dic[record.id] = n_reads
+    return readcount_dic

iss/generator.py

@@ -555,9 +555,11 @@ def load_readcount_or_abundance(
             coverage_dic = abundance.parse_abundance_file(coverage_file)
             complete_genomes_dic = {k: v for k, v in coverage_dic.items() if k not in draft}
             draft_dic = abundance.expand_draft_abundance(coverage_dic, draft, mode="coverage")
-            abundance_dic = {**complete_genomes_dic, **draft_dic}
+            coverage_dic = {**complete_genomes_dic, **draft_dic}
         else:
-            abundance_dic = abundance.parse_abundance_file(coverage_file)
+            coverage_dic = abundance.parse_abundance_file(coverage_file)
+        # generate new n_reads from coverage
+        readcount_dic = abundance.to_readcount(coverage_dic, error_model.read_length, genome_file)
     elif coverage in abundance_dispatch:
         # todo coverage distribution with --draft
         logger.info("Using %s coverage distribution" % coverage)