Interventions per Gbp and percent assembly in chromosomes stats

Author: jgrg

src/tola/assembly/assembly_stats.py

@@ -1,30 +1,40 @@
 import csv
 import io
 import logging
+from typing import TypeAlias
 
 from tola.assembly.assembly import Assembly, AssemblyDict
 
 log = logging.getLogger(__name__)
 
 
+RankedNameLengths: TypeAlias = dict[int, dict[str, int]]
+
+
 class AssemblyStatsError(Exception):
     """Error from AssemblyStats"""
 
 
 class AssemblyStats:
-    def __init__(self, autosome_prefix: str = "SUPER_") -> None:
+    def __init__(self, autosome_prefix: str = "SUPER_"):
         self.autosome_prefix = autosome_prefix
         self.input_assembly: Assembly | None = None
         self.cuts = 0
-        self.breaks = 0
-        self.joins = 0
+        self.breaks = None
+        self.joins = None
+        self.interventions_per_gbp = None
+        self.percent_assembly_in_chromosomes = None
         self.per_assembly_stats = {}
         self.assembly_scaffold_lengths = {}
 
-    def make_stats(self, output_assemblies: AssemblyDict) -> None:
+    def make_stats(self, output_assemblies: AssemblyDict):
         if not self.input_assembly:
             msg = "Missing input_assembly attribute"
             raise AssemblyStatsError(msg)
+        self.__build_junction_stats(output_assemblies)
+        self.__build_length_stats(output_assemblies)
+
+    def __build_junction_stats(self, output_assemblies: AssemblyDict):
 
         # These stats are going to be wrong if re-curating assemblies with
         # Fragment names beginning with "SUPER_".
@@ -63,6 +73,25 @@ def make_stats(self, output_assemblies: AssemblyDict) -> None:
                     "manual_joins": len((junc_set - input_asm_set) & total_joins),
                 }
 
+    def __build_length_stats(self, output_assemblies: AssemblyDict):
+        input_asm_length = self.input_assembly.fragments_length
+
+        # Calculate the number of breaks and joins made per Gbp of the input assembly
+        self.interventions_per_gbp = round(
+            (self.breaks + self.joins) / (input_asm_length / 1e9), 3
+        )
+
+        # Calculate the percentage of the assembly that was placed in chromosomes
+        tier_1_and_2_length = 0
+        for hap, asm in output_assemblies.items():
+            ranked_scffld_lengths = self.get_assembly_scaffold_lengths(hap, asm)
+            for rank in (1, 2):
+                if scffld_lengths := ranked_scffld_lengths.get(rank):
+                    tier_1_and_2_length += sum(scffld_lengths.values())
+        self.percent_assembly_in_chromosomes = round(
+            100 * (tier_1_and_2_length / input_asm_length), 1
+        )
+
     def log_curation_stats(self):
         cut_plural = "cut in a contig" if self.cuts == 1 else "cuts in contigs"
         break_plural = "break at a gap" if self.breaks == 1 else "breaks at gaps"
@@ -71,6 +100,10 @@ def log_curation_stats(self):
             f"Curation made {self.cuts} {cut_plural}, {self.breaks}"
             f" {break_plural} and {self.joins} {join_plural}"
         )
+        log.info(
+            f"Assembly placed in chromosomes = {self.percent_assembly_in_chromosomes}%"
+        )
+        log.info(f"Interventions per Gbp = {self.interventions_per_gbp}")
 
     def ranked_scaffolds(self, asm: Assembly):
         ranked_scaffolds = {}
@@ -79,7 +112,7 @@ def ranked_scaffolds(self, asm: Assembly):
             ranked_scaffolds.setdefault(rank, []).append(scffld)
         return ranked_scaffolds
 
-    def build_assembly_scaffold_lengths(self, asm: Assembly):
+    def build_assembly_scaffold_lengths(self, asm: Assembly) -> RankedNameLengths:
         ranked_scaffolds = self.ranked_scaffolds(asm)
 
         ranked_names_lengths = {}
@@ -105,7 +138,11 @@ def build_assembly_scaffold_lengths(self, asm: Assembly):
 
         return ranked_names_lengths
 
-    def get_assembly_scaffold_lengths(self, asm_key: str | None, asm: Assembly):
+    def get_assembly_scaffold_lengths(
+        self,
+        asm_key: str | None,
+        asm: Assembly,
+    ) -> RankedNameLengths:
         scaff_lengths = self.assembly_scaffold_lengths.get(asm_key)
         if not scaff_lengths:
             scaff_lengths = self.assembly_scaffold_lengths[asm_key] = (
@@ -250,7 +287,7 @@ def log_sanity_checks(self, hap_asm: dict[str | None, Assembly]) -> None:
                     log.warning(msg)
 
     def check_consistent_autosome_count(
-        self, hap_asm: dict[str | None, Assembly]
+        self, hap_asm: AssemblyDict
     ) -> list[str] | None:
         chr_counts = {}
         for hap, asm in hap_asm.items():
@@ -267,9 +304,7 @@ def check_consistent_autosome_count(
                 ]
         return None
 
-    def check_for_large_haplotigs(
-        self, hap_asm: dict[str | None, Assembly]
-    ) -> list[str] | None:
+    def check_for_large_haplotigs(self, hap_asm: AssemblyDict) -> list[str] | None:
         htigs = hap_asm.get("Haplotig")
         if not htigs:
             return None

src/tola/assembly/scripts/pretext_to_asm.py

@@ -481,7 +481,7 @@ def write_chr_csv_files(
                 csv_fh.write(chr_names)
 
 
-def write_info_yaml(output_file, stats, out_assemblies: AssemblyDict, clobber):
+def write_info_yaml(output_file, stats: AssemblyStats, out_assemblies: AssemblyDict, clobber):
     asm_stats = stats.per_assembly_stats
     info = {"assemblies": asm_stats}
     if len(asm_stats) > 1:
@@ -492,6 +492,8 @@ def write_info_yaml(output_file, stats, out_assemblies: AssemblyDict, clobber):
     if h_asm := out_assemblies.get("Haplotig"):
         haplotig_count = len(h_asm.scaffolds)
     info["manual_haplotig_removals"] = haplotig_count
+    info["percent_assembly_in_chromosomes"] = stats.percent_assembly_in_chromosomes
+    info["interventions_per_gbp"] = stats.interventions_per_gbp
 
     yaml_file = output_file.with_name(output_file.stem + ".info.yaml")
     with get_output_filehandle(yaml_file, clobber) as yaml_fh:

tests/data/bChlMac1_3/bChlMac1-pretext-to-tpf.3.info.yaml

@@ -8,3 +8,5 @@ assemblies:
 manual_breaks: 15
 manual_joins: 29
 manual_haplotig_removals: 2
+percent_assembly_in_chromosomes: 95.6
+interventions_per_gbp: 18.446

tests/data/bChlMac1_3/bChlMac1-pretext-to-tpf.3.log

@@ -50,5 +50,7 @@ bChlMac1-pretext-to-tpf.3.haplotigs
           4,748,144  bp total
 
 Curation made 1 cut in a contig, 15 breaks at gaps and 29 joins
+Assembly placed in chromosomes = 95.6%
+Interventions per Gbp = 18.446
 Mismatch in autosome count between Hap1 = 36 and Hap2 = 35
 Haplotig H_1 (Scaffold_75) is 4,560,100 bp which is longer than the shortest chromosome (292,585 bp)

tests/data/csSphGirg1/csSphGirg1-pretext-to-tpf.info.yaml

@@ -3,3 +3,5 @@ assemblies:
     manual_breaks: 3
     manual_joins: 4
 manual_haplotig_removals: 0
+percent_assembly_in_chromosomes: 96.6
+interventions_per_gbp: 15.651

tests/data/csSphGirg1/csSphGirg1-pretext-to-tpf.log

@@ -19,3 +19,5 @@ csSphGirg1-pretext-to-tpf.1.contaminants
          15,413,452  bp total
 
 Curation made 0 cuts in contigs, 3 breaks at gaps and 4 joins
+Assembly placed in chromosomes = 96.6%
+Interventions per Gbp = 15.651

tests/data/eaAstIrre1/eaAstIrre1-pretext-to-tpf.info.yaml

@@ -3,3 +3,5 @@ assemblies:
     manual_breaks: 71
     manual_joins: 97
 manual_haplotig_removals: 44
+percent_assembly_in_chromosomes: 93.4
+interventions_per_gbp: 336.493

tests/data/eaAstIrre1/eaAstIrre1-pretext-to-tpf.log

@@ -148,3 +148,5 @@ eaAstIrre1-pretext-to-tpf.1.contaminants
           2,742,977  bp total
 
 Curation made 20 cuts in contigs, 71 breaks at gaps and 97 joins
+Assembly placed in chromosomes = 93.4%
+Interventions per Gbp = 336.493

tests/data/idDilFebr1/idDilFebr1-pretext-to-tpf.info.yaml

@@ -3,3 +3,5 @@ assemblies:
     manual_breaks: 20
     manual_joins: 31
 manual_haplotig_removals: 7
+percent_assembly_in_chromosomes: 93.8
+interventions_per_gbp: 163.495

tests/data/idDilFebr1/idDilFebr1-pretext-to-tpf.log

@@ -28,3 +28,5 @@ idDilFebr1-pretext-to-tpf.1.additional_haplotigs
           3,119,814  bp total
 
 Curation made 0 cuts in contigs, 20 breaks at gaps and 31 joins
+Assembly placed in chromosomes = 93.8%
+Interventions per Gbp = 163.495