Version 3.4 (#59)

- Fix bug where MM/ML tags are off in supplementary alignments 
- Fix bug where the last base of a read is sometimes used in phasing 
- Fix bug in variant calling where in some reads a wrong base is chosen between primary and supplementary alignments 
- Fix bug where some alleles may include redundant haplotype names 
- Fix bug where it's not possible to lower the minimum variant frequency for variant calling for individual target regions. This change affects opn1lw and ikbkg. 
- Improve large deletion calling in reads. This change may affect all targets, particularly ikbkg and pms2. 
- Rename haplotypes to label gene1 and gene2 for regions with fusion calling(GBA, CYP2D6, CYP11B1 and CFH/CFHR3) 
- For smn1, consider more scenarios for adjusting SMN2 copy number 
- For ncf1, adjust copy number for the scenario where three haplotypes are found and all are present at two copies.
- Update two copy haplotypes for CFHclust

Author: xiao-chen-xc

paraphase/__init__.py

@@ -1 +1 @@
-__version__ = "3.3.4"
+__version__ = "3.4.0"

paraphase/genes/cfhclust.py

@@ -52,6 +52,8 @@ def call(self):
                     len(self.cfh["final_haplotypes"]),
                     len(self.cfhr3["final_haplotypes"]),
                 )
+                if total_cn < self.cfh["total_cn"] or total_cn < self.cfhr3["total_cn"]:
+                    two_cp_haps = []
 
         return self.GeneCall(
             total_cn,

paraphase/genes/ncf1_phaser.py

@@ -119,7 +119,23 @@ def call(self):
         haplotypes = tmp
 
         two_cp_haps = []
-        if counter_gene == 1:
+        # scenario where only three haplotypes are found, possibly each at CN2
+        if total_cn == 3:
+            if self.mdepth is not None:
+                probs = self.depth_prob(self.region_avg_depth.median, self.mdepth)
+                if probs is not None and probs[2] + probs[3] > 0.95:
+                    two_cp_haps = list(ass_haps.values())
+                    for hap in two_cp_haps:
+                        total_cn += 1
+                        if "pseudo" in hap:
+                            counter_pseudo += 1
+                        else:
+                            counter_gene += 1
+            if counter_gene == 1 and counter_pseudo == 2:
+                counter_gene = None
+                total_cn = None
+
+        elif counter_gene == 1:
             two_cp_hap_candidate = self.compare_depth(haplotypes, ass_haps)
             if two_cp_hap_candidate == []:
                 # check if one haplotype has more reads than others
@@ -131,18 +147,15 @@ def call(self):
                 counter_gene += 1
                 total_cn += 1
 
-        if self.mdepth is not None:
+        if self.mdepth is not None and counter_gene is not None:
             prob = self.depth_prob(gene_reads, self.mdepth / 2)
-            if prob[0] < 0.9 and counter_gene == 1:
-                counter_gene = None
-                total_cn = None
-            if prob[0] > 0.95 and counter_gene > 1 and two_cp_haps != []:
-                counter_gene = None
-                total_cn = None
-        # scenario where only three haplotypes are found, possibly each at CN2
-        if counter_gene == 1 and counter_pseudo == 2 and total_cn == 3:
-            counter_gene = None
-            total_cn = None
+            if prob is not None:
+                if prob[0] < 0.9 and counter_gene == 1:
+                    counter_gene = None
+                    total_cn = None
+                if prob[0] > 0.95 and counter_gene > 1 and two_cp_haps != []:
+                    counter_gene = None
+                    total_cn = None
 
         # homozygous case
         if total_cn == 0:

paraphase/genes/pms2_phaser.py

@@ -28,7 +28,10 @@ def call(self):
                 phase_region=f"{self.genome_build}:{self.nchr}:{self.left_boundary}-{self.right_boundary}",
             )
         self.get_homopolymer()
-        self.find_big_deletion(min_size=2900)
+        self.find_big_deletion(min_size=1000)
+        pms2cl_clip_site = self.clip_3p_positions[0]
+        self.find_clip_site()
+        self.clip_3p_positions = [pms2cl_clip_site]
 
         if self.deletion1_size is not None:
             self.del1_reads, self.del1_reads_partial = self.get_long_del_reads(
@@ -52,6 +55,7 @@ def call(self):
         # for distinguishing pms2 from pms2cl
         raw_read_haps = self.get_haplotypes_from_reads(
             clip_buffer=50,
+            min_clip_len=100,
             check_clip=True,
             kept_sites=homo_sites_to_add,
             add_sites=self.add_sites,
@@ -93,13 +97,22 @@ def call(self):
         for hap in ass_haps:
             clip_position = self.get_3pclip_from_hap(hap)
             if clip_position is None:
-                counter_unknown += 1
-                hap_name = f"{self.gene}_unknownhap{counter_unknown}"
-            elif clip_position in self.clip_3p_positions:
+                num_pms2_bases = 0
+                for i, base in enumerate(hap):
+                    pos = int(self.het_sites[i].split("_")[0])
+                    if pos > pms2cl_clip_site + 50 and base not in ["0", "x"]:
+                        num_pms2_bases += 1
+                if num_pms2_bases >= 3:
+                    counter_gene += 1
+                    hap_name = f"{self.gene}_pms2hap{counter_gene}"
+                else:
+                    counter_unknown += 1
+                    hap_name = f"{self.gene}_unknownhap{counter_unknown}"
+            elif clip_position == pms2cl_clip_site:
                 counter_pseudo += 1
                 hap_name = f"{self.gene}_pms2clhap{counter_pseudo}"
             else:
-                assert clip_position == 0
+                # assert clip_position == 0
                 counter_gene += 1
                 hap_name = f"{self.gene}_pms2hap{counter_gene}"
             tmp.setdefault(hap, hap_name)

paraphase/genes/smn1_phaser.py

@@ -311,21 +311,22 @@ def adjust_smn1_cn(self, smn1_cn, smn2_cn, hcn, ass_haps, read_counts, smn1_haps
             haploid_depth = genome_depth / 2
             depth1 = self.smn1_reads_splice
             copy_number_probs = self.depth_prob(depth1, haploid_depth)
-            if smn1_cn == 1:
-                # when there is only one haplotype, check if depth is
-                # consistent with haploid depth
-                copy_one_prob = copy_number_probs[0]
-                if copy_one_prob < 0.05:
-                    two_cp_haps = list(smn1_haps.values())
-                    return 2, two_cp_haps
-                if copy_one_prob < 0.25:
-                    new_smn1_cn = None
-            elif smn1_cn == 2:
-                # scenario where two two-copy alleles are identical
-                copy_four_prob = copy_number_probs[3]
-                if copy_four_prob > 0.75:
-                    two_cp_haps = list(smn1_haps.values())
-                    return 4, two_cp_haps
+            if copy_number_probs is not None:
+                if smn1_cn == 1:
+                    # when there is only one haplotype, check if depth is
+                    # consistent with haploid depth
+                    copy_one_prob = copy_number_probs[0]
+                    if copy_one_prob < 0.05:
+                        two_cp_haps = list(smn1_haps.values())
+                        return 2, two_cp_haps
+                    if copy_one_prob < 0.25:
+                        new_smn1_cn = None
+                elif smn1_cn == 2:
+                    # scenario where two two-copy alleles are identical
+                    copy_four_prob = copy_number_probs[3]
+                    if copy_four_prob > 0.75:
+                        two_cp_haps = list(smn1_haps.values())
+                        return 4, two_cp_haps
         if smn1_cn in [2, 3] and read_counts is not None and self.targeted is False:
             # check if one smn1 haplotype has more reads than others
             haps = list(read_counts.keys())
@@ -335,7 +336,7 @@ def adjust_smn1_cn(self, smn1_cn, smn2_cn, hcn, ass_haps, read_counts, smn1_haps
             if cp2_hap in smn1_haps:
                 others_max = sorted(counts, reverse=True)[1]
                 probs = self.depth_prob(max_count, others_max)
-                if probs[0] < 0.0005 and others_max >= 10:
+                if probs is not None and probs[0] < 0.0005 and others_max >= 10:
                     two_cp_haps.append(smn1_haps[cp2_hap])
                     return smn1_cn + 1, two_cp_haps
 
@@ -346,12 +347,13 @@ def adjust_smn1_cn(self, smn1_cn, smn2_cn, hcn, ass_haps, read_counts, smn1_haps
                 copy_number_probs = self.depth_prob(depth1, haploid_depth)
                 # when there is only one haplotype, check if depth is
                 # consistent with haploid depth
-                copy_one_prob = copy_number_probs[0]
-                if copy_one_prob < 0.25:
-                    new_smn1_cn = None
-                    if smn2_cn > 1 and copy_one_prob < 0.05:
-                        two_cp_haps = list(smn1_haps.values())
-                        return 2, two_cp_haps
+                if copy_number_probs is not None:
+                    copy_one_prob = copy_number_probs[0]
+                    if copy_one_prob < 0.25:
+                        new_smn1_cn = None
+                        if smn2_cn > 1 and copy_one_prob < 0.05:
+                            two_cp_haps = list(smn1_haps.values())
+                            return 2, two_cp_haps
         # if we see more haplotypes in other regions of the gene
         if smn1_cn == 1 and hcn > len(ass_haps):
             if self.has_smn2 is False:
@@ -361,7 +363,7 @@ def adjust_smn1_cn(self, smn1_cn, smn2_cn, hcn, ass_haps, read_counts, smn1_haps
                 new_smn1_cn = None
         return new_smn1_cn, two_cp_haps
 
-    def adjust_smn2_cn(self, smn1_cn, smn2_cn, smn2_haps):
+    def adjust_smn2_cn(self, smn1_cn, smn2_cn, read_counts, smn2_haps):
         """
         Adjust SMN2 CN if there is only one haplotype found.
         """
@@ -377,19 +379,43 @@ def adjust_smn2_cn(self, smn1_cn, smn2_cn, smn2_haps):
             depth1 = self.smn2_reads_splice
             copy_number_probs = self.depth_prob(depth1, haploid_depth)
             # if smn1 cn is 3, then no need to adjust smn2 cn
-            if (
-                smn2_cn == 1
-                and len(self.smn2_del_reads_partial) <= 1
-                and (smn1_cn is None or smn1_cn <= 2)
-            ):
-                # when there is only one haplotype, check if depth is
-                # consistent with haploid depth
-                copy_one_prob = copy_number_probs[0]
-                if copy_one_prob < 0.25:
-                    new_smn2_cn = None
-                    if copy_one_prob < 0.05:
+            if copy_number_probs is not None:
+                if (
+                    smn2_cn == 1
+                    and len(self.smn2_del_reads_partial) <= 1
+                    and (smn1_cn is None or smn1_cn <= 2)
+                ):
+                    # when there is only one haplotype, check if depth is
+                    # consistent with haploid depth
+                    copy_one_prob = copy_number_probs[0]
+                    if copy_one_prob < 0.25:
+                        new_smn2_cn = None
+                        if copy_one_prob < 0.05:
+                            two_cp_haps = list(smn2_haps.values())
+                            return 2, two_cp_haps
+                if smn1_cn == 0 and smn2_cn == 2:
+                    # scenario where two two-copy alleles are identical
+                    copy_four_prob = copy_number_probs[3]
+                    if copy_four_prob > 0.75:
                         two_cp_haps = list(smn2_haps.values())
-                        return 2, two_cp_haps
+                        return 4, two_cp_haps
+        if (
+            smn1_cn in [0, 1]
+            and smn2_cn in [2, 3]
+            and read_counts is not None
+            and self.targeted is False
+        ):
+            # check if one smn2 haplotype has more reads than others
+            haps = list(read_counts.keys())
+            counts = list(read_counts.values())
+            max_count = max(counts)
+            cp2_hap = haps[counts.index(max_count)]
+            if cp2_hap in smn2_haps:
+                others_max = sorted(counts, reverse=True)[1]
+                probs = self.depth_prob(max_count, others_max)
+                if probs is not None and probs[0] < 0.0005 and others_max >= 10:
+                    two_cp_haps.append(smn2_haps[cp2_hap])
+                    return smn2_cn + 1, two_cp_haps
         if smn1_cn is not None:
             if smn2_cn == 1 and smn1_cn == 2 and len(self.smn2_del_reads_partial) <= 1:
                 haploid_depth = self.smn1_reads_splice / smn1_cn
@@ -398,12 +424,13 @@ def adjust_smn2_cn(self, smn1_cn, smn2_cn, smn2_haps):
                     copy_number_probs = self.depth_prob(depth1, haploid_depth)
                     # when there is only one haplotype, check if depth is
                     # consistent with haploid depth
-                    copy_one_prob = copy_number_probs[0]
-                    if copy_one_prob < 0.25:
-                        new_smn2_cn = None
-                        if copy_one_prob < 0.05:
-                            two_cp_haps = list(smn2_haps.values())
-                            return 2, two_cp_haps
+                    if copy_number_probs is not None:
+                        copy_one_prob = copy_number_probs[0]
+                        if copy_one_prob < 0.25:
+                            new_smn2_cn = None
+                            if copy_one_prob < 0.05:
+                                two_cp_haps = list(smn2_haps.values())
+                                return 2, two_cp_haps
         return new_smn2_cn, two_cp_haps
 
     def get_best_match(
@@ -629,7 +656,7 @@ def call(self):
                     cp2_hap = haps[counts.index(max_count)]
                     others_max = sorted(counts, reverse=True)[1]
                     probs = self.depth_prob(max_count, others_max)
-                    if probs[0] < 0.05 and others_max >= 10:
+                    if probs is not None and probs[0] < 0.05 and others_max >= 10:
                         two_cp_haps.append(haps_to_compare[cp2_hap])
         for hap in two_cp_haps:
             if "smn1hap" in hap:
@@ -645,7 +672,9 @@ def call(self):
         for hap in two_cp_haps_smn1:
             if hap not in two_cp_haps:
                 two_cp_haps.append(hap)
-        smn2_cn, two_cp_haps_smn2 = self.adjust_smn2_cn(smn1_cn_old, smn2_cn, smn2_haps)
+        smn2_cn, two_cp_haps_smn2 = self.adjust_smn2_cn(
+            smn1_cn_old, smn2_cn, read_counts, smn2_haps
+        )
         for hap in two_cp_haps_smn2:
             if hap not in two_cp_haps:
                 two_cp_haps.append(hap)