Version 3.2.1 (#34)

- Fix problem working with CRAM
Other minor changes:
- Add RN tag to output bam. RN stands for region name, indicating reads used to analyze one region (paralog group)
- Sort some output fields in the JSON, making it consistent from run to run
- Clean up reported alleles, filtering out cases where all haplotypes are linked into one allele, or there are missing haplotypes when two alleles are reported
- Fix getting sample ID from bam header when there are blank spaces
- Fix logic for writing homozygous haplotypes in VCF (no haplotypes phased -> no haplotypes phased and no heterozygous variant sites)
- Update two_cp_haplotypes when adjusting total_cn from 2 to 4 based on depth
- Use all reads instead of unique reads for variant calling at edges of clipped haplotypes

Author: xiao-chen-xc

paraphase/__init__.py

@@ -1 +1 @@
-__version__ = "3.2.0"
+__version__ = "3.2.1"

paraphase/genes/cfc1_phaser.py

@@ -102,5 +102,5 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
-            None,
+            self.init_het_sites,
         )

paraphase/genes/cfhclust.py

@@ -72,5 +72,7 @@ def call(self):
             None,
             None,
             None,
+            None,
+            None,
             fusions,
         )

paraphase/genes/f8_phaser.py

@@ -5,13 +5,12 @@
 from collections import namedtuple
 import pysam
 from ..phaser import Phaser
+from paraphase.prepare_bam_and_vcf import pysam_handle
 
 
 class F8Phaser(Phaser):
     new_fields = copy.deepcopy(Phaser.fields)
     new_fields.remove("gene_cn")
-    new_fields.remove("alleles_final")
-    new_fields.remove("hap_links")
     new_fields.insert(3, "sv_called")
     new_fields.insert(3, "flanking_summary")
     new_fields.insert(3, "exon1_to_exon22_depth")
@@ -33,6 +32,9 @@ def __init__(
         Phaser.__init__(
             self, sample_id, outdir, args, genome_depth, genome_bam, sample_sex
         )
+        self.reference_fasta = None
+        if args is not None:
+            self.reference_fasta = args.reference
 
     def set_parameter(self, config):
         super().set_parameter(config)
@@ -41,11 +43,11 @@ def set_parameter(self, config):
             "extract_regions"
         ].split()
 
-    def get_read_positions(self, min_extension=5000):
+    def get_read_positions(self, genome_bamh, min_extension=5000):
         """Get mapped region of the part of reads not overlapping repeat"""
         dpos5 = {}
         dpos3 = {}
-        genome_bamh = pysam.AlignmentFile(self.genome_bam, "rb")
+
         for i, extract_region in enumerate(
             [self.extract_region1, self.extract_region2, self.extract_region3]
         ):
@@ -79,18 +81,17 @@ def get_read_positions(self, min_extension=5000):
                             dpos3.setdefault(read_name, []).append(pos_name)
                         else:
                             dpos5.setdefault(read_name, []).append(pos_name)
-        genome_bamh.close()
         return dpos5, dpos3
 
     def call(self):
         if self.check_coverage_before_analysis() is False:
             return self.GeneCall()
 
-        genome_bamh = pysam.AlignmentFile(self.genome_bam, "rb")
+        genome_bamh = pysam_handle(self.genome_bam, self.reference_fasta)
         e1_e22_depth = self.get_regional_depth(genome_bamh, self.depth_region)[0].median
-        genome_bamh.close()
 
-        dpos5, dpos3 = self.get_read_positions()
+        dpos5, dpos3 = self.get_read_positions(genome_bamh)
+        genome_bamh.close()
         self.get_homopolymer()
         self.get_candidate_pos()
 
@@ -114,6 +115,7 @@ def call(self):
 
         self.het_sites = sorted(list(self.candidate_pos))
         self.remove_noisy_sites()
+        self.init_het_sites = [a for a in self.het_sites]
 
         raw_read_haps = self.get_haplotypes_from_reads(
             check_clip=True,
@@ -215,6 +217,8 @@ def call(self):
             e1_e22_depth,
             flanking_sum,
             sv_hap,
+            None,
+            None,
             hcn,
             original_haps,
             self.het_sites,
@@ -227,4 +231,5 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
+            self.init_het_sites,
         )

paraphase/genes/hba_phaser.py

@@ -5,6 +5,7 @@
 from collections import namedtuple
 import pysam
 from ..phaser import Phaser
+from paraphase.prepare_bam_and_vcf import pysam_handle
 
 
 class HbaPhaser(Phaser):
@@ -31,6 +32,9 @@ def __init__(
         Phaser.__init__(
             self, sample_id, outdir, args, genome_depth, genome_bam, sample_sex
         )
+        self.reference_fasta = None
+        if args is not None:
+            self.reference_fasta = args.reference
 
     def set_parameter(self, config):
         super().set_parameter(config)
@@ -39,7 +43,7 @@ def set_parameter(self, config):
     def call(self):
         if self.check_coverage_before_analysis() is False:
             return self.GeneCall()
-        genome_bamh = pysam.AlignmentFile(self.genome_bam, "rb")
+        genome_bamh = pysam_handle(self.genome_bam, self.reference_fasta)
         surrounding_region_depth = self.get_regional_depth(
             genome_bamh, self.depth_region
         )[0].median
@@ -227,13 +231,15 @@ def call(self):
 
         # phase
         alleles = []
+        linked_haps = []
         hap_links = {}
         if self.to_phase is True:
             (
                 alleles,
                 hap_links,
                 _,
                 _,
+                linked_haps,
             ) = self.phase_alleles(
                 uniquely_supporting_reads,
                 nonuniquely_supporting_reads,
@@ -245,6 +251,7 @@ def call(self):
             if len(alleles) == 1:
                 if len(alleles[0]) == len(ass_haps):
                     alleles = []
+                    linked_haps = []
         self.close_handle()
 
         return self.GeneCall(
@@ -268,4 +275,6 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
+            self.init_het_sites,
+            linked_haps,
         )

paraphase/genes/ikbkg_phaser.py

@@ -71,6 +71,7 @@ def call(self):
         kept_sites = [self.add_sites[1]]
         self.het_sites = sorted(self.het_sites)
         self.remove_noisy_sites()
+        self.init_het_sites = [a for a in self.het_sites]
         homo_sites_to_add = self.add_homo_sites(min_no_var_region_size=6000)
         # print(homo_sites_to_add)
         kept_sites += homo_sites_to_add
@@ -189,13 +190,23 @@ def call(self):
             hap_links,
             _,
             _,
+            linked_haps,
         ) = self.phase_alleles(
             uniquely_supporting_reads,
             nonuniquely_supporting_reads,
             raw_read_haps,
             ass_haps,
             reverse=self.is_reverse,
         )
+        # all haplotypes are phased into one allele
+        if (
+            len(alleles) == 1
+            and sorted(alleles[0]) == sorted(ass_haps.values())
+            and self.sample_sex is not None
+            and self.sample_sex == "female"
+        ):
+            alleles = []
+            linked_haps = []
 
         if gene_counter == 0 or pseudo_counter == 0:
             total_cn = None
@@ -228,4 +239,6 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
+            self.init_het_sites,
+            linked_haps,
         )

paraphase/genes/ncf1_phaser.py

@@ -52,6 +52,7 @@ def call(self):
         self.get_candidate_pos()
         self.het_sites = sorted(list(self.candidate_pos))
         self.remove_noisy_sites()
+        self.init_het_sites = [a for a in self.het_sites]
 
         raw_read_haps = self.get_haplotypes_from_reads(add_sites=self.add_sites)
 
@@ -155,4 +156,5 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
+            self.init_het_sites,
         )

paraphase/genes/neb_phaser.py

@@ -40,6 +40,7 @@ def call(self):
         self.get_candidate_pos()
         self.het_sites = sorted(list(self.candidate_pos))
         self.remove_noisy_sites()
+        self.init_het_sites = [a for a in self.het_sites]
         homo_sites_to_add = self.add_homo_sites()
         raw_read_haps = self.get_haplotypes_from_reads(
             kept_sites=homo_sites_to_add,
@@ -120,13 +121,15 @@ def call(self):
             tri3.append(tri3[0])
 
         alleles = []
+        linked_haps = []
         hap_links = {}
         if two_cp_haps == []:
             (
                 alleles,
                 hap_links,
                 _,
                 _,
+                linked_haps,
             ) = self.phase_alleles(
                 uniquely_supporting_reads,
                 nonuniquely_supporting_reads,
@@ -139,6 +142,7 @@ def call(self):
         # incorrect phasing suggests haplotypes with cn > 1
         if len(alleles) == 1 and sorted(alleles[0]) == sorted(ass_haps.values()):
             alleles = []
+            linked_haps = []
             total_cn = None
         elif len(tri1) > 2 or len(tri3) > 2:
             total_cn = None
@@ -155,6 +159,7 @@ def call(self):
                     total_cn = 6
                     two_cp_haps = [tri2[0]]
                     alleles = []
+                    linked_haps = []
 
         # homozygous case
         if total_cn == 0:
@@ -181,4 +186,6 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
+            self.init_het_sites,
+            linked_haps,
         )

paraphase/genes/opn1lw_phaser.py

@@ -1,13 +1,13 @@
 # paraphase
 # Author: Xiao Chen <xchen@pacificbiosciences.com>
 
-
+import copy
 from collections import namedtuple
 from ..phaser import Phaser
 
 
 class Opn1lwPhaser(Phaser):
-    fields = (
+    new_fields = (
         "total_cn",
         "opn1lw_cn",
         "opn1mw_cn",
@@ -25,23 +25,13 @@ class Opn1lwPhaser(Phaser):
         "directional_links",
         "links_loose",
         "alleles_raw",
-        "highest_total_cn",
-        "assembled_haplotypes",
-        "sites_for_phasing",
-        "unique_supporting_reads",
-        "het_sites_not_used_in_phasing",
-        "homozygous_sites",
-        "haplotype_details",
-        "nonunique_supporting_reads",
-        "read_details",
-        "genome_depth",
-        "region_depth",
-        "sample_sex",
     )
+    # new_fields = copy.deepcopy(Phaser.fields)
+    new_fields += tuple(Phaser.fields[6:])
     GeneCall = namedtuple(
         "GeneCall",
-        fields,
-        defaults=(None,) * len(fields),
+        new_fields,
+        defaults=(None,) * len(new_fields),
     )
     pathogenic_haps = [
         "LIAVA",
@@ -189,6 +179,7 @@ def call(self):
         self.get_candidate_pos(min_vaf=0.08)
         self.het_sites = sorted(list(self.candidate_pos))
         self.remove_noisy_sites()
+        self.init_het_sites = [a for a in self.het_sites]
         homo_sites_to_add = self.add_homo_sites()
         raw_read_haps = self.get_haplotypes_from_reads(
             check_clip=True,
@@ -284,13 +275,15 @@ def call(self):
                 hap_links,
                 directional_links,
                 directional_links_loose,
+                linked_haps,
             ) = self.phase_alleles(
                 uniquely_supporting_reads,
                 nonuniquely_supporting_reads,
                 raw_read_haps,
                 ass_haps,
                 reverse=self.is_reverse,
             )
+            alleles = linked_haps
 
             # incorrect phasing suggests haplotypes with cn > 1
             if (
@@ -455,4 +448,6 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
+            self.init_het_sites,
+            alleles_1st_2nd,
         )

paraphase/genes/pms2_phaser.py

@@ -42,6 +42,7 @@ def call(self):
         self.get_candidate_pos()
         self.het_sites = sorted(list(self.candidate_pos))
         self.remove_noisy_sites()
+        self.init_het_sites = [a for a in self.het_sites]
         homo_sites_to_add = self.add_homo_sites(min_no_var_region_size=6000)
         # for distinguishing pms2 from pms2cl
         raw_read_haps = self.get_haplotypes_from_reads(
@@ -149,5 +150,5 @@ def call(self):
             self.mdepth,
             self.region_avg_depth._asdict(),
             self.sample_sex,
-            None,
+            self.init_het_sites,
         )