Fix #94

Author: inumanag

.gitignore

@@ -40,3 +40,4 @@ aldy_old/*
 scripts
 aldy/indelpost/*.c
 build/*
+/_*

aldy/__main__.py

@@ -31,11 +31,15 @@
 def get_version():
     return "{} {}".format(
         platform.system() if platform.system() != "Darwin" else "macOS",
-        platform.platform()[6:]
-        if platform.system() == "Linux"
-        else platform.mac_ver()[0]
-        if platform.system() == "Darwin"
-        else platform.platform(),
+        (
+            platform.platform()[6:]
+            if platform.system() == "Linux"
+            else (
+                platform.mac_ver()[0]
+                if platform.system() == "Darwin"
+                else platform.platform()
+            )
+        ),
     )
 
 

aldy/common.py

@@ -5,7 +5,7 @@
 
 
 from typing import Iterable, Any, List
-import pkg_resources
+import importlib.resources
 import re
 import time
 import pprint
@@ -115,7 +115,7 @@ def samtools(self, pad_left=500, pad_right=1, prefix="") -> str:
         :param prefix: Chromosome prefix."""
 
         return "{}:{}-{}".format(
-            prefix + self.chr, self.start - pad_left, self.end + pad_right
+            prefix + self.chr, max(1, self.start - pad_left), self.end + pad_right
         )
 
     def __str__(self):
@@ -202,7 +202,7 @@ def script_path(key: str) -> str:
     components = key.split("/")
     if len(components) < 2:
         raise AldyException(f'"{key}"" is not valid resource name')
-    return pkg_resources.resource_filename(components[0], "/".join(components[1:]))
+    return str(importlib.resources.files(components[0]) / "/".join(components[1:]))
 
 
 def colorize(text: str, color: str = "green") -> str:

aldy/diplotype.py

@@ -176,15 +176,19 @@ def write_vcf(
                     "GT": "|".join(str(all_mutations[m][mi][i]) for i in range(nall)),
                     "DP": str(coverage[m]),
                     "MA": ",".join(
-                        f"*{minor.solution[i].major}"
-                        if all_mutations[m][mi][i] > 0
-                        else "-"
+                        (
+                            f"*{minor.solution[i].major}"
+                            if all_mutations[m][mi][i] > 0
+                            else "-"
+                        )
                         for i in range(nall)
                     ),
                     "MI": ",".join(
-                        f"*{minor.solution[i].minor}"
-                        if all_mutations[m][mi][i] > 0
-                        else "-"
+                        (
+                            f"*{minor.solution[i].minor}"
+                            if all_mutations[m][mi][i] > 0
+                            else "-"
+                        )
                         for i in range(nall)
                     ),
                 }
@@ -307,6 +311,7 @@ def key(x):
     solution.set_diplotype(diplotype)
     return diplotype
 
+
 def estimate_cpic(gene: Gene, solution: MinorSolution) -> str:
     """Calculate the CPIC functionality for a minor solution.
     :returns: CPIC functionality.
@@ -335,6 +340,3 @@ def estimate_cpic(gene: Gene, solution: MinorSolution) -> str:
             if len(i) == 1:
                 return (0, i[0])
     return (0, "indeterminate")
-
-
-

aldy/gene.py

@@ -280,9 +280,11 @@ def get_functional(self, mut, infer=True) -> Optional[str]:
             if self.seq[pos] != op[0]:
                 log.warn(f"Bad mutation: {op[0]} != {self.seq[pos]}")
             seq = "".join(
-                self.seq[s:pos] + op[2] + self.seq[pos + 1 : e]
-                if s <= pos < e
-                else self.seq[s:e]
+                (
+                    self.seq[s:pos] + op[2] + self.seq[pos + 1 : e]
+                    if s <= pos < e
+                    else self.seq[s:e]
+                )
                 for s, e in self.exons
             )
             amino = seq_to_amino(seq)
@@ -409,6 +411,8 @@ def _init_basic(self, yml) -> None:
                 seq[pos - 1] = nuc
             self.seq = "".join(seq)
 
+        if self.genome not in yml["reference"]["mappings"]:
+            raise AldyException("Gene {self.name} not compatible with {self.genome}")
         self.chr, start, end, strand, cigar = yml["reference"]["mappings"][self.genome]
         self.chr_to_ref = {}
         self.ref_to_chr = {}
@@ -433,12 +437,14 @@ def _init_basic(self, yml) -> None:
         self._lookup_range = (start - 1, end - 1)
         self._lookup_seq = "".join(
             (
-                rev_comp(self.seq[self.chr_to_ref[i]])
-                if self.strand < 0
-                else self.seq[self.chr_to_ref[i]]
+                (
+                    rev_comp(self.seq[self.chr_to_ref[i]])
+                    if self.strand < 0
+                    else self.seq[self.chr_to_ref[i]]
+                )
+                if i in self.chr_to_ref
+                else "N"
             )
-            if i in self.chr_to_ref
-            else "N"
             for i in range(*self._lookup_range)
         )
 
@@ -474,7 +480,7 @@ def _init_regions(self, yml) -> None:
             for name, coord in yml["structure"]["regions"][self.genome].items():
                 if name[0] == "e" and name[1:].isdigit():  # this is exon
                     num_exons += 1
-                if not coord[i * 2] <= coord[i * 2 + 1]:
+                if i * 2 + 1 < len(coord) and not coord[i * 2] <= coord[i * 2 + 1]:
                     raise AldyException(
                         f"Malformed YML file {self.name} (structure:regions:{name})"
                     )

aldy/genotype.py

@@ -7,7 +7,7 @@
 from typing import List, Optional, Any, Set, Dict
 import os
 import sys
-import pkg_resources
+import importlib.resources
 import datetime
 import time
 
@@ -98,17 +98,17 @@ def genotype(
 
     avail_genes = []
     if gene_db == "all":
-        avail_genes = pkg_resources.resource_listdir("aldy.resources", "genes")
+        avail_genes = importlib.resources.files("aldy.resources.genes").iterdir()
         avail_genes = [
             i[:-4]
             for i in avail_genes
             if len(i) > 4 and i.endswith(".yml") and not i.startswith("pharma-")
         ]
         avail_genes = sorted(avail_genes)
     elif gene_db == "pharmacoscan":
-        avail_genes = pkg_resources.resource_listdir(
-            "aldy.resources.genes", "pharmacoscan"
-        )
+        avail_genes = importlib.resources.files(
+            "aldy.resources.genes.pharmacoscan"
+        ).iterdir()
         avail_genes = [
             f"pharmacoscan/{i[:-4]}" for i in avail_genes if i.endswith(".yml")
         ]
@@ -362,7 +362,10 @@ def genotype(
                 cpic += f"; cpic={cpic_fun}"
                 if gene.cpic_scores and cpic_fun != "indeterminate":
                     cpic += f"; cpic_score={cpic_score}"
-            print(f"#Solution {i + 1}: {minor_sol._solution_nice()}{cpic}", file=output_file)
+            print(
+                f"#Solution {i + 1}: {minor_sol._solution_nice()}{cpic}",
+                file=output_file,
+            )
             diplotype.write_decomposition(
                 sample.name, gene, sample.coverage, i + 1, minor_sol, output_file
             )

aldy/indelpost/alleles.py

@@ -32,18 +32,22 @@ def phase_nearby_variants(
     # no phasable variants
     variants_to_phase = contig.mismatches + contig.non_target_indels
     if not variants_to_phase:
-        return  make_target_obj_from_contig(target, indexed_contig)
+        return make_target_obj_from_contig(target, indexed_contig)
 
     # phase all phasables within the target exon (hard phasing)
     if hard:
         variants_list = []
-        cleaned, variant_list = precleaning(indexed_contig, variants_list, target_pos_on_contig, pileup, target)
+        cleaned, variant_list = precleaning(
+            indexed_contig, variants_list, target_pos_on_contig, pileup, target
+        )
         return greedy_phasing(target, cleaned)
     else:
-        indexed_contig, variants_to_phase = precleaning(indexed_contig, variants_to_phase, target_pos_on_contig, pileup)
+        indexed_contig, variants_to_phase = precleaning(
+            indexed_contig, variants_to_phase, target_pos_on_contig, pileup
+        )
 
     if not variants_to_phase:
-        return  make_target_obj_from_contig(target, indexed_contig)
+        return make_target_obj_from_contig(target, indexed_contig)
     else:
         variants_in_non_targets, mut_frac = variants_in_non_target_pileup(
             pileup, target, basequalthresh, to_complex
@@ -69,13 +73,25 @@ def phase_nearby_variants(
 
     remove_deletables(indexed_contig, lt_end, target_pos_on_contig, rt_end)
 
-    mismatches_to_phase = [var for var in variants_to_phase if not var.is_indel and indexed_contig.get(var.pos, False)]
-    non_target_indels_to_phase = [var for var in variants_to_phase if var.is_indel and indexed_contig.get(var.pos, False) and var != target]
+    mismatches_to_phase = [
+        var
+        for var in variants_to_phase
+        if not var.is_indel and indexed_contig.get(var.pos, False)
+    ]
+    non_target_indels_to_phase = [
+        var
+        for var in variants_to_phase
+        if var.is_indel and indexed_contig.get(var.pos, False) and var != target
+    ]
 
     if variants_to_phase:
         if not non_target_indels_to_phase:
             peak_locs = locate_mismatch_cluster_peaks(
-                indexed_contig, mismatches_to_phase, target, snv_neighborhood, to_complex
+                indexed_contig,
+                mismatches_to_phase,
+                target,
+                snv_neighborhood,
+                to_complex,
             )
 
             if peak_locs:
@@ -93,10 +109,20 @@ def phase_nearby_variants(
             if max(target_len, non_target_max_len) < 4:
                 indel_neighborhood = int(indel_neighborhood / 2) + 1
 
-            remove_common_substrings(indexed_contig, target_pos_on_contig, indel_neighborhood)
+            remove_common_substrings(
+                indexed_contig, target_pos_on_contig, indel_neighborhood
+            )
 
-            lt_end = end_point(indexed_contig, mismatches_to_phase, target, snv_neighborhood, left=True)
-            rt_end = end_point(indexed_contig, mismatches_to_phase, target, snv_neighborhood, left=False)
+            lt_end = end_point(
+                indexed_contig, mismatches_to_phase, target, snv_neighborhood, left=True
+            )
+            rt_end = end_point(
+                indexed_contig,
+                mismatches_to_phase,
+                target,
+                snv_neighborhood,
+                left=False,
+            )
 
             remove_deletables(indexed_contig, lt_end, target_pos_on_contig, rt_end)
 
@@ -111,10 +137,13 @@ def phase_nearby_variants(
 def make_target_obj_from_contig(target, indexed_contig):
     try:
         data = indexed_contig[target.pos]
-        return Variant(target.chrom, target.pos, data[0], data[1], target.reference).normalize()
+        return Variant(
+            target.chrom, target.pos, data[0], data[1], target.reference
+        ).normalize()
     except:
         return target.normalize()
 
+
 def greedy_phasing(target, indexed_contig):
 
     cpos = 0
@@ -142,7 +171,9 @@ def seq_complexity(contig, snv_neighborhood, indel_neighorhood):
     )
 
 
-def precleaning(genome_indexed_contig, variants_list, target_pos, pileup, limit_to_target_exon=True):
+def precleaning(
+    genome_indexed_contig, variants_list, target_pos, pileup, limit_to_target_exon=True
+):
     lt_loci, rt_loci = [], []
 
     # filter low qual loci
@@ -234,9 +265,10 @@ def locate_mismatch_cluster_peaks(
     else:
         return None
 
-
     lt_peak_pos = target.pos if lt_peak_pos == -np.inf else lt_peak_pos
-    rt_peak_pos = target.pos + len(target.ref) - 1 if rt_peak_pos == np.inf else rt_peak_pos
+    rt_peak_pos = (
+        target.pos + len(target.ref) - 1 if rt_peak_pos == np.inf else rt_peak_pos
+    )
 
     return (lt_peak_pos - 1, rt_peak_pos + 1)
 
@@ -248,7 +280,7 @@ def calc_peak(indexed_contig, mismatches, target, snv_neighborhood, left):
         loci = [k for k, v in indexed_contig.items() if k <= target_pos][::-1]
         snv_loci = [var.pos for var in mismatches if var.pos < target_pos]
     else:
-        del_adjust = len(target.ref) -1
+        del_adjust = len(target.ref) - 1
         loci = [k for k, v in indexed_contig.items() if k > target_pos + del_adjust]
         snv_loci = [var.pos for var in mismatches if var.pos > target_pos]
 
@@ -319,7 +351,7 @@ def variants_in_non_target_pileup(pileup, target, basequalthresh, to_complex):
     nontarget_pileup = [
         findall_mismatches(read, end_trim=10)
         for read in pileup
-        if not read["is_target" ] and read["is_covering"] and not read["is_dirty"]
+        if not read["is_target"] and read["is_covering"] and not read["is_dirty"]
     ]
 
     if not nontarget_pileup:
@@ -340,27 +372,26 @@ def variants_in_non_target_pileup(pileup, target, basequalthresh, to_complex):
     indels = [
         indel
         for indel, cnt in Counter(indels).items()
-        if (cnt > 2 and cnt/ len(nontarget_pileup) > 0.15)
-        or cnt > 5
+        if (cnt > 2 and cnt / len(nontarget_pileup) > 0.15) or cnt > 5
     ]
 
     mismatches = [
         Variant(target.chrom, v[0], v[1], v[2], target.reference)
         for read in nontarget_pileup
-        for v in read["mismatches"] if v[3] > basequalthresh
+        for v in read["mismatches"]
+        if v[3] > basequalthresh
     ]
 
-    nontarget_pileup_vol = sum(
-        max(0, len(read["ref_seq"]) - 20) for read in nontarget_pileup
-    ) + 1
+    nontarget_pileup_vol = (
+        sum(max(0, len(read["ref_seq"]) - 20) for read in nontarget_pileup) + 1
+    )
 
     mutation_frac = (len(mismatches) + len(indels)) / nontarget_pileup_vol
 
     mismatches = [
         var
         for var, cnt in Counter(mismatches).items()
-        if (cnt > 2 and cnt / len(nontarget_pileup) > 0.15)
-        or cnt > 5
+        if (cnt > 2 and cnt / len(nontarget_pileup) > 0.15) or cnt > 5
     ]
 
     return set(indels + mismatches), mutation_frac
@@ -396,9 +427,9 @@ def get_freq(freqinfo):
 def remove_deletables(indexed_contig, lt_end, target_pos, rt_end):
     tmp = indexed_contig.copy()
 
-    #if lt_end == -np.inf:
+    # if lt_end == -np.inf:
     #    lt_end = target_pos - 1
-    #if rt_end == np.inf:
+    # if rt_end == np.inf:
     #    rt_end = target_pos + 1
 
     for k, v in tmp.items():
@@ -448,7 +479,6 @@ def trim_common(indexed_contig, commons, max_common_str_len, left):
         else:
             start = search_nearest_lt_locus(indexed_contig, sub_str[0], left)
 
-
         end = sub_str[-1]
         start_event = indexed_contig[start]
         end_event = indexed_contig[end]
@@ -579,9 +609,7 @@ def end_point(indexed_contig, mismatches, target, snv_neighborhood, left):
             return peak_pos + 1
 
 
-
 def get_end_most_indel(indexed_contig, target):
     for k, v in indexed_contig.items():
         if len(v[0]) != len(v[1]):
             return Variant(target.chrom, k, v[0], v[1], target.reference)
-