New pluviometer

bin/FeatureOutputWriter.py

@@ -0,0 +1,193 @@
+from MultiCounter import MultiCounter
+from Bio.SeqFeature import SeqFeature
+from Bio.SeqFeature import ExactPosition
+from collections import defaultdict
+from typing import TextIO
+from utils import BASE_TYPES, MATCH_MISMATCH_TYPES
+
+FEATURE_OUTPUT_FIELDS = [
+    "SeqID",
+    "Parents",
+    "FeatureID",
+    "Type",
+    "Start",
+    "End",
+    "Strand",
+    "CoveredSites",
+    f"GenomeBases[{','.join(BASE_TYPES)}]",
+    f"SiteBasePairs[{','.join(MATCH_MISMATCH_TYPES)}]",
+    f"ReadBasePairs[{','.join(MATCH_MISMATCH_TYPES)}]",
+]
+
+FEATURE_METADATA_OUTPUT_FIELDS = [
+    "SeqID",
+    "ParentsIDs",
+    "FeatureID",
+    "Type",
+    "Start",
+    "End",
+    "Strand",
+]
+
+FEATURE_DATA_OUTPUT_FIELDS = [
+    "CoveredSites",
+    f"GenomeBases[{','.join(BASE_TYPES)}]",
+    f"SiteBasePairs[{','.join(MATCH_MISMATCH_TYPES)}]",
+    f"ReadBasePairs[{','.join(MATCH_MISMATCH_TYPES)}]",
+]
+
+AGGREGATE_METADATA_OUTPUT_FIELDS = [
+    "SeqID",
+    "ParentsIDs",
+    "FeatureID",
+    "ParentType",
+    "AggregateType",
+]
+
+AGGREGATE_DATA_OUTPUT_FIELDS = [
+    "CoveredSites",
+    f"GenomeBases[{','.join(BASE_TYPES)}]",
+    f"SiteBasePairs[{','.join(MATCH_MISMATCH_TYPES)}]",
+    f"ReadBasePairs[{','.join(MATCH_MISMATCH_TYPES)}]",
+]
+
+STR_ZERO_BASE_FREQS = ",".join('0' for _ in range(len(BASE_TYPES)))
+STR_ZERO_EDIT_FREQS = ",".join('0' for _ in range(len(MATCH_MISMATCH_TYPES)))
+
+
+def make_parent_path(parent_list: list[str]) -> str:
+    """
+    Create a path string from an ordered list of parent IDs.
+    The separator is a comma, chosen because it is one of the few invalid characters in tag=value entries of the attributes field in the GFF3 format.
+
+    Consult the GFF3 specification for details: https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md
+    """
+    return ','.join(parent_list)
+
+
+class RainFileWriter:
+    def __init__(
+        self, handle: TextIO, metadata_fields: list[str], data_fields: list[str]
+    ):
+        self.handle = handle
+        self.metadata_fields: list[str] = metadata_fields
+        self.n_metadata: int = len(self.metadata_fields)
+        self.data_fields: list[str] = data_fields
+        self.n_data: int = len(self.data_fields)
+
+        return None
+
+    def write_header(self) -> int:
+        b: int = self.handle.write("\t".join(self.metadata_fields))
+        b += self.handle.write("\t")
+        b += self.handle.write("\t".join(self.data_fields))
+        b += self.handle.write("\n")
+
+        return b
+
+    def write_comment(self, comment: str) -> int:
+        b = self.handle.write("# ")
+        b += self.handle.write(comment)
+        b += self.handle.write("\n")
+
+        return b
+
+    def write_metadata(self, *metadata_values) -> int:
+        b: int = 0
+        for val in metadata_values:
+            b += self.handle.write(val)
+            b += self.handle.write("\t")
+
+        return b
+
+    def write_data(self, *data_values) -> int:
+        b: int = 0
+        for val in data_values[:-1]:
+            b += self.handle.write(val)
+            b += self.handle.write("\t")
+
+        b += self.handle.write(data_values[-1])
+        b += self.handle.write("\n")
+
+        return b
+
+
+class FeatureFileWriter(RainFileWriter):
+    def __init__(self, handle: TextIO):
+        super().__init__(
+            handle, FEATURE_METADATA_OUTPUT_FIELDS, FEATURE_DATA_OUTPUT_FIELDS
+        )
+
+        return None
+
+    def write_metadata(self, record_id: str, feature: SeqFeature) -> int:
+        return super().write_metadata(
+            record_id,
+            make_parent_path(feature.parent_list),
+            feature.id,
+            feature.type,
+            str(feature.location.parts[0].start + ExactPosition(1)),
+            str(feature.location.parts[-1].end),
+            str(feature.location.strand),
+        )
+
+    def write_row_with_data(
+        self, record_id: str, feature: SeqFeature, counter: MultiCounter
+    ) -> int:
+        return self.write_metadata(record_id, feature) + self.write_data(
+            str(counter.genome_base_freqs.sum()),
+            ",".join(map(str, counter.genome_base_freqs.flat)),
+            ",".join(map(str, counter.edit_site_freqs.flat)),
+            ",".join(map(str, counter.edit_read_freqs.flat)),
+        )
+
+    def write_row_without_data(self, record_id: str, feature: SeqFeature) -> int:
+        return self.write_metadata(record_id, feature) + self.write_data(
+            '0', STR_ZERO_BASE_FREQS, STR_ZERO_EDIT_FREQS, STR_ZERO_EDIT_FREQS
+        )
+
+class AggregateFileWriter(RainFileWriter):
+    def __init__(self, handle: TextIO):
+        super().__init__(
+            handle, AGGREGATE_METADATA_OUTPUT_FIELDS, AGGREGATE_DATA_OUTPUT_FIELDS
+        )
+
+        return None
+
+    def write_metadata(self, seq_id: str, feature: SeqFeature, aggregate_type: str) -> int:
+        return super().write_metadata(seq_id, make_parent_path(feature.parent_list), feature.id, feature.type, aggregate_type)
+
+    def write_rows_with_feature_and_data(self, record_id: str, feature: SeqFeature, counter_dict: defaultdict[str,MultiCounter]) -> int:
+        b: int = 0
+
+        for aggregate_type, aggregate_counter in counter_dict.items():
+            b += self.write_metadata(record_id, feature, aggregate_type)
+            b += self.write_data(
+                str(aggregate_counter.genome_base_freqs.sum()),
+                ",".join(map(str, aggregate_counter.genome_base_freqs.flat)),
+                ",".join(map(str, aggregate_counter.edit_site_freqs.flat)),
+                ",".join(map(str, aggregate_counter.edit_read_freqs.flat)),
+            )
+
+        return b
+
+    def write_rows_with_data(
+            self,
+            record_id: str,
+            parent_list: list[str],
+            feature_id: str,
+            feature_type: str,
+            counter_dict: defaultdict[str,MultiCounter]
+            ) -> int:
+        b: int = 0
+
+        for aggregate_type, aggregate_counter in counter_dict.items():
+            b += super().write_metadata(record_id, make_parent_path(parent_list), feature_id, feature_type, aggregate_type)
+            b += self.write_data(
+                str(aggregate_counter.genome_base_freqs.sum()),
+                ",".join(map(str, aggregate_counter.genome_base_freqs.flat)),
+                ",".join(map(str, aggregate_counter.edit_site_freqs.flat)),
+                ",".join(map(str, aggregate_counter.edit_read_freqs.flat)),
+            )
+
+        return b

bin/MultiCounter.py

@@ -0,0 +1,67 @@
+from utils import SiteVariantData
+import numpy as np
+from numpy.typing import NDArray
+from SiteFilter import SiteFilter
+from typing import TextIO
+
+class MultiCounter:
+    """Holds the counter data and logic for a feature, feature aggregate, or record"""
+
+    def __init__(self, site_filter: SiteFilter) -> None:
+        """
+        Tallies of the numbers of reads per edit type
+        This is a numpy matrix where the rows represent the reference base and the columns the edited base
+        Rows and column indices correspond to bases in alphabetic order (ACGT)
+        Row-columns corresponding to the same base (e.g. (0,0) -> (A,A)) do not represent edits, and should remain 0
+        """
+        self.edit_read_freqs: NDArray[np.int64] = np.zeros((5, 5), dtype=np.int64)
+        self.edit_site_freqs: NDArray[np.int64] = np.zeros((5, 5), dtype=np.int64)
+
+        self.genome_base_freqs: NDArray[np.int64] = np.zeros(5, dtype=np.int64)
+
+        self.filter = site_filter
+
+        return None
+
+    def update(self, variant_data: SiteVariantData) -> None:
+        """Increment the counters from the data in a SiteVariantData object."""
+        i: int = variant_data.reference
+
+        self.edit_read_freqs[i, :] += variant_data.frequencies
+
+        self.filter.apply(variant_data)
+        self.edit_site_freqs[i, :] += self.filter.frequencies
+
+        self.genome_base_freqs[i] += 1
+
+        return None
+
+    def merge(self, other_counter: "MultiCounter") -> None:
+        """
+        Add to this counter the values of another.
+        """
+        self.edit_read_freqs[:] += other_counter.edit_read_freqs
+        self.edit_site_freqs[:] += other_counter.edit_site_freqs
+        self.genome_base_freqs[:] += other_counter.genome_base_freqs
+
+        return None
+
+    def report(self, output_handle: TextIO) -> int:
+        b = 0
+
+        # Write the number of covered sites
+        b += output_handle.write(str(self.genome_base_freqs.sum()))
+        b += output_handle.write("\t")
+
+        # Write the base frequencies in the genome
+        b += write_base_array(output_handle, self.genome_base_freqs)
+        b += output_handle.write("\t")
+
+        # Write edited sites
+        b += write_edit_array(output_handle, self.edit_site_freqs)
+        b += output_handle.write("\t")
+
+        # Write edit frequencies
+        b += write_edit_array(output_handle, self.edit_read_freqs)
+
+        return b

bin/SeqFeature_extensions.py

@@ -0,0 +1,104 @@
+# Custom methods and attributes for Bio.SeqFeature
+from Bio.SeqFeature import SeqFeature, SimpleLocation, CompoundLocation
+from utils import location_union
+from typing import Optional
+import logging
+
+logger = logging.getLogger(__name__)
+
+setattr(SeqFeature, "level", 0)
+
+
+def get_transcript_like(self: SeqFeature) -> list[tuple[str, str, int]]:
+    """
+    Return a list with information about sub-features that are transcript-like (i.e. their contain children of type "exon" or "CDS").
+
+    List items are tuples that contain the ID of the transcript-like feature, the type of the transcript-like feature, and the total exon or CDS of the transcript-like feature.
+    """
+    transcript_like_list: list[tuple[str, str, int]] = []
+    for transcript_candidate in self.sub_features:
+        total_exon_length: int = 0
+        total_cds_length: int = 0
+        for child in transcript_candidate.sub_features:
+            if child.type == "exon":
+                total_exon_length += len(child)
+            elif child.type == "CDS":
+                total_cds_length += len(child)
+
+        if total_cds_length > 0:
+            transcript_like_list.append(
+                (transcript_candidate.id, "CDS", total_cds_length)
+            )
+        elif total_exon_length > 0:
+            transcript_like_list.append(
+                (transcript_candidate.id, "exon", total_exon_length)
+            )
+
+    return transcript_like_list
+
+
+setattr(SeqFeature, "get_transcript_like", get_transcript_like)
+
+setattr(SeqFeature, "parent_list", [""])
+
+def make_chimaera(self: SeqFeature) -> None:
+    """
+    If the feature contains
+    """
+    if hasattr(self, "sub_features"):
+        if len(self.sub_features) == 0:
+            return None
+    else:
+        return None
+
+    transcript_like_list: list[SeqFeature] = list(
+        filter(
+            lambda transcript: any(map(lambda part: part.type == "CDS", transcript.sub_features)),
+            self.sub_features,
+        )
+    )
+
+    if len(transcript_like_list) == 0:
+        chimaeric_type: str = "exon"
+        transcript_like_list: list[SeqFeature] = list(
+            filter(
+                lambda transcript: any(map(lambda part: part.type == "exon", transcript.sub_features)),
+                self.sub_features,
+            )
+        )
+    else:
+        chimaeric_type: str = "CDS"
+
+    if len(transcript_like_list) == 0:
+        return None
+
+
+    target_locations: list[SimpleLocation | CompoundLocation] = []
+    for transcript in transcript_like_list:
+        target_locations.extend(
+            list(map(
+                lambda part: part.location,
+                filter(lambda part: part.type == chimaeric_type, transcript.sub_features),
+            ))
+        )
+
+    chimaeric_location: SimpleLocation | CompoundLocation = location_union(
+        target_locations
+    )
+    logging.info(f"Created {chimaeric_type} chimaera of feature {self.id}: {len(transcript_like_list)} transcripts were merged into one transcript of {len(chimaeric_location.parts)} elements")
+
+    chimaeric_feature: SeqFeature = SeqFeature(
+        location=chimaeric_location,
+        type=chimaeric_type + "-chimaera",
+        id=self.id + "-chimaera",
+        qualifiers={"Parent": self.id},
+    )
+
+    chimaeric_feature.sub_features = []
+
+    self.sub_features.append(chimaeric_feature)
+
+    return None
+
+
+setattr(SeqFeature, "make_chimaera", make_chimaera)

bin/SiteFilter.py

@@ -0,0 +1,21 @@
+from utils import SiteVariantData
+import numpy as np
+from numpy.typing import NDArray
+
+class SiteFilter:
+    def __init__(self, cov_threshold: int, edit_threshold: int) -> None:
+        self.cov_threshold: int = cov_threshold
+        self.edit_threshold: int = edit_threshold
+        self.frequencies: NDArray[np.int32] = np.zeros(5, np.int32)
+
+    def apply(self, variant_data: SiteVariantData) -> None:
+        if variant_data.coverage >= self.cov_threshold:
+            np.copyto(
+                self.frequencies,
+                variant_data.frequencies * variant_data.frequencies
+                >= self.edit_threshold,
+            )
+        else:
+            self.frequencies.fill(0)
+
+        return None