diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e62ee3ac..5ff8c11d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -47,7 +47,7 @@ jobs: python -m site python -m pip install -U pip setuptools wheel python -m pip install -e ".[tests]" - coverage run --omit "src/readfish/read_until/*.py" -pm pytest + coverage run --omit "src/readfish/read_until/*.py,src/readfish/entry_points/targets.py" -pm pytest - name: Upload coverage data uses: actions/upload-artifact@v3 with: diff --git a/README.md b/README.md index c0a0978d..cb8585cc 100644 --- a/README.md +++ b/README.md @@ -555,6 +555,7 @@ And for our Awesome Logo please checkout out [@tim_bassford](https://twitter.com # Changelog ## Unreleased changes 1. Change the default `unblock_duration` on the `Analysis` class to use `DEFAULT_UNBLOCK` value defined in `_cli_args.py`. Change type on the Argparser for `--unblock-duration` to float. (#313) +1. Big dog Duplex feature - adds ability to select duplex reads that cover a target region. See pull request for details [(#324)](https://github.com/LooseLab/readfish/pull/324) ## 2023.1.1 1. Fix Readme Logo link 🥳 (#296) 1. Fix bug where we had accidentally started requiring barcoded TOMLs to specify a region. Thanks to @jamesemery for catching this. (#299) diff --git a/src/readfish/__about__.py b/src/readfish/__about__.py index 4e896e6f..9ba8864e 100644 --- a/src/readfish/__about__.py +++ b/src/readfish/__about__.py @@ -1,4 +1,4 @@ """__about__.py Version of the read until software """ -__version__ = "2023.1.1" +__version__ = "2024.2.0a" diff --git a/src/readfish/_cli_args.py b/src/readfish/_cli_args.py index 2102afa0..898d92b8 100644 --- a/src/readfish/_cli_args.py +++ b/src/readfish/_cli_args.py @@ -5,8 +5,21 @@ ``BASE_ARGS`` are the minimal required arguments for _all_ entry points as they used for initialising loggers. ``DEVICE_BASE_ARGS`` are the set of arguments that are used for connecting to a sequencer (device) and some other related settings for selective sequencing scripts. """ + +from enum import Enum, unique from readfish._utils import nice_join + +@unique +class Chemistry(Enum): + #: For the "smarter" version of duplex - does this read map to the previous reads opposite strand on the same contig. Won't work for no map based decisions + DUPLEX = "duplex" + #: Normal simplex chemistry - no duplex override shenanigans + SIMPLEX = "simplex" + #: Simple duplex - if we are going to unblock a read given the previous read on the same channel was stop receiving, sequence the current read instead. + DUPLEX_SIMPLE = "duplex_simple" + + DEFAULT_SERVER_HOST = "127.0.0.1" DEFAULT_SERVER_PORT = None DEFAULT_LOG_FORMAT = "%(asctime)s %(name)s %(message)s" @@ -133,4 +146,15 @@ type=int, ), ), + ( + "--chemistry", + dict( + help="**EXPERIMENTAL** Choose between duplex and simplex chemistry mode. duplex_simple accepts a read if the previous channels read was stop receiving," + "duplex checks that the previous reads alignment was on the same contig and opposite strand. default: SIMPLEX", + required=False, + type=str, + default=Chemistry.SIMPLEX, + choices=[chemistry.value for chemistry in Chemistry], + ), + ), ) + BASE_ARGS diff --git a/src/readfish/_loggers.py b/src/readfish/_loggers.py index c77a03f1..30897882 100644 --- a/src/readfish/_loggers.py +++ b/src/readfish/_loggers.py @@ -4,6 +4,7 @@ from logging.handlers import QueueHandler, QueueListener import queue from typing import Callable +from readfish.__about__ import __version__ def setup_logger( @@ -109,3 +110,4 @@ def print_args( for attr in dir(args): if attr[0] != "_" and attr not in exclude and attr.lower() == attr: printer(f"{attr}={getattr(args, attr)!r}") + printer(f"Version={__version__}") diff --git a/src/readfish/entry_points/targets.py b/src/readfish/entry_points/targets.py index 09da4973..a8ba1ea8 100644 --- a/src/readfish/entry_points/targets.py +++ b/src/readfish/entry_points/targets.py @@ -11,6 +11,10 @@ These chunks of raw data are processed by the basecaller to produce FASTA, which is then aligned against the chosen reference genome. The result of the alignment is used, along with the targets provided in the :doc:`TOML ` file, to make a decision on each read. +In the new **experimental** Duplex mode, it is possible to override a decision for a read based on the action taken for a previous read. +This is done by passing `--chemistry` and setting either duplex, or duplex simple. +`duplex_simple` accepts a read if the previous channels read was stop receiving, `duplex` checks that the previous reads alignment was on the same contig and opposite strand. +The default chemistry is simplex. Running this should result in a very short (<1kb, ideally 400-600 bases) unblock peak at the start of a read length histogram and longer sequenced reads. @@ -22,6 +26,15 @@ --log-file rf.log \\ --debug-log chunks.tsv +Example experimental duplex command:: + + readfish targets --device X3 \\ + --experiment-name "test" \\ + --toml my_exp.toml \\ + --log-file rf.log \\ + --debug-log chunks.tsv + --chemistry duplex + In the debug_log chunks.tsv file, if this argument is passed, each line represents detailed information about a batch of read signal that has been processed in an iteration. @@ -66,6 +79,7 @@ """ + # Core imports from __future__ import annotations import argparse @@ -81,7 +95,7 @@ from minknow_api import protocol_service # Library -from readfish._cli_args import DEVICE_BASE_ARGS, DEFAULT_UNBLOCK +from readfish._cli_args import DEVICE_BASE_ARGS, Chemistry from readfish._read_until_client import RUClient from readfish._config import Action, Conf, make_decision, _Condition from readfish._statistics import ReadfishStatistics @@ -92,7 +106,13 @@ Severity, ) from readfish.plugins.abc import AlignerABC, CallerABC -from readfish.plugins.utils import Decision, PreviouslySentActionTracker, Result +from readfish.plugins.utils import ( + Decision, + PreviouslySentActionTracker, + Result, + DuplexTracker, + Strand, +) _help = "Run targeted sequencing" @@ -123,6 +143,9 @@ ), ), ) +# When sequencing in duplex mode, overriding a decided `Action` on a currently sequenced molecule +# is not allowed if the previous molecules decision was one of these. +DISALLOWED_DUPLEX_DECISIONS = {Decision.first_read_override, Decision.duplex_override} class Analysis: @@ -131,14 +154,16 @@ class Analysis: function that is run threaded in the run function at the base of this file. Arguments listed in the __init__ docs. - :param client: An instance of the ReadUntilClient object - :param conf: An instance of the Conf object - :param logger: The command level logger for this module - :param debug_log: Whether to output the Debug Log. log Name is generated. Defaults to False. - :param throttle: The number of seconds interval between requests to the ReadUntilClient, defaults to 0.1 - :param unblock_duration: Time, in seconds, to apply unblock voltage, defaults to 0.5 - :param dry_run: If True unblocks are replaced with `stop_receiving` commands, defaults to False - :param toml: The path to the toml file containing experiment conf. Used for reloading, defaults to "a.toml" + :param client: An instance of the ReadUntilClient object. + :param conf: An instance of the Conf object. + :param logger: The command level logger for this module. + :param debug_log: Whether to output the Debug Log. log Name is generated. + :param throttle: The time interval (seconds) between requests to the ReadUntilClient. + :param unblock_duration: Time, in seconds, to apply unblock voltage. + :param dry_run: If True unblocks are replaced with `stop_receiving` commands. + :param toml: The path to the toml file containing experiment conf. Used as the path for checking if the TOML needs reloading. + :param chemistry: Instance of Chemistry Enum, representing the chemistry of the run (Simplex/Duplex). Used for + decision making on strands that may be part of a duplex pair. """ def __init__( @@ -147,10 +172,11 @@ def __init__( conf: Conf, logger: logging.Logger, debug_log: bool, - throttle: float = 0.1, - unblock_duration: float = DEFAULT_UNBLOCK, - dry_run: bool = False, - toml: str = "a.toml", + throttle: float, + unblock_duration: float, + dry_run: bool, + toml: str, + chemistry: Chemistry, ): self.client = client self.conf = conf @@ -161,7 +187,7 @@ def __init__( self.dry_run = dry_run self.live_toml = Path(f"{toml}_live").resolve() self.run_information = self.client.connection.protocol.get_run_info() - + self.chemistry = chemistry # Generate a run specific read log read_log_name = ( f"{self.run_information.run_id}_readfish.tsv" if debug_log else None @@ -194,6 +220,8 @@ def __init__( # This is an object to keep track of the last action sent to the client for each channel self.previous_action_tracker = PreviouslySentActionTracker() + # Keep track of previous alignments + self.duplex_tracker = DuplexTracker() # We assume that sequencing is already running. # If the run is not in sequencing phase when the read until loop starts will @@ -263,9 +291,11 @@ def check_override_action( Checks include: 1. If the read is in a control region, the action is always stop_receiving. 1. If the read is below the minimum chunks, use value in toml or default to proceed - 1. If the read is above the maximum chunks, use value in toml or default unblock + 1. If the read is above the maximum chunks, use value in toml or default unblock--throttle 1. First read seen for channel and readfish started during sequencing, override to stop_receiving 1. If action is unblock and we are dry-running, override to stop_receiving + 1. If we are running in duplex chemistry, check the previous reads final decision and Action, and potentially sequence + the current read, instead of unblocking it. :param control: Indicates read from a channel in a control region :param action: What action was decided for this read before any meddling @@ -306,18 +336,66 @@ def check_override_action( # previous_action will be None if the read has not been seen before. previous_action = self.previous_action_tracker.get_action(result.channel) action_overridden = False + # If --duplex flag override decisions made based on the strand and contig alignment of the previous read. + # Unfinished bruv + if ( + self.chemistry is Chemistry.DUPLEX + # Easy checks first, so wdon't do more complex processing unless we have to + and action == Action.unblock + and previous_action is Action.stop_receiving + ): + # Check if we think this read is possibly duplex + possible_duplex = any( + self.duplex_tracker.possible_duplex( + result.channel, result.read_id, al.ctg, al.strand + ) + for al in result.alignment_data + ) + # Check the previous decision for this channel was not already an override + previous_decision_allowed = ( + self.duplex_tracker.get_previous_decision(result.channel) + not in DISALLOWED_DUPLEX_DECISIONS + ) + if possible_duplex and previous_decision_allowed: + self.logger.debug( + f"Overriding read {result.read_id} as it is possibly second half of a duplex" + f"- previous read action {previous_action}, current_action: {action}," + f" previous_decision: {self.duplex_tracker.get_previous_decision(result.channel)}" + ) + action_overridden = True + result.decision = Decision.duplex_override + action = Action.stop_receiving + # Duplex + elif ( + self.chemistry is Chemistry.DUPLEX_SIMPLE + and previous_action is Action.stop_receiving + and action is Action.unblock + ): + previous_decision_allowed = ( + self.duplex_tracker.get_previous_decision(result.channel) + not in DISALLOWED_DUPLEX_DECISIONS + ) + if previous_decision_allowed: + self.logger.debug( + f"Overriding to duplex - previous read action {previous_action}, current_action: {action}," + f" previous_decision: {self.duplex_tracker.get_previous_decision(result.channel)}" + ) + action = Action.stop_receiving + action_overridden = True + result.decision = Decision.duplex_override + # Override to stop receiving if this is the first read ona channel and we started mid sequencing if previous_action is None and self.readfish_started_during_sequencing: self.logger.debug( f"This is the first suitable read chunk from channel {result.channel}. Translocated read length unknown, sequencing." ) action_overridden = True + result.decision = Decision.first_read_override action = Action.stop_receiving if action is Action.stop_receiving: stop_receiving_action_list.append((result.channel, result.read_number)) - # Add decided Action - self.previous_action_tracker.add_action(result.channel, action) + elif action is Action.unblock: if self.dry_run: # Log an 'unblock' action to previous action, but send a 'stop receiving' to prevent further read processing. @@ -327,8 +405,21 @@ def check_override_action( unblock_batch_action_list.append( (result.channel, result.read_number, result.read_id) ) + + # If we have made a final decision for this read and we shouldn't see it again! + if action is Action.unblock or action is Action.stop_receiving: # Add decided Action self.previous_action_tracker.add_action(result.channel, action) + # Add duplex based tracking if we are in duplex mode + if self.chemistry is Chemistry.DUPLEX_SIMPLE: + self.duplex_tracker.set_decision(result.channel, result.decision) + elif self.chemistry is Chemistry.DUPLEX: + self.duplex_tracker.set_decision(result.channel, result.decision) + self.duplex_tracker.set_alignments( + result.channel, + [(al.ctg, Strand(al.strand)) for al in result.alignment_data], + ) + return ( previous_action, action_overridden, @@ -530,6 +621,7 @@ def run( throttle=args.throttle, dry_run=args.dry_run, toml=args.toml, + chemistry=Chemistry(args.chemistry), ) # begin readfish function diff --git a/src/readfish/plugins/utils.py b/src/readfish/plugins/utils.py index df964722..b07684c7 100644 --- a/src/readfish/plugins/utils.py +++ b/src/readfish/plugins/utils.py @@ -30,6 +30,16 @@ class Strand(Enum): #: Reverse strand reverse = "-" + def __invert__(self): + """Flip the strand, using bitwise NOT (~) + + >>> ~Strand.forward + + >>> ~Strand.reverse + + """ + return Strand.forward if self is Strand.reverse else Strand.reverse + STRANDS = { 1: Strand.forward, @@ -160,11 +170,16 @@ def count_dict_elements(d: dict[Any]) -> int: """ return sum( ( - count_dict_elements(v) if isinstance(v, dict) - # If v is a list, tuple, sequence, dict etc., return the length of the container filtering out any empty sun elements, - else len(list(filterfalse(partial(is_empty), v))) - if isinstance(v, Container) - else 1 + ( + count_dict_elements(v) + if isinstance(v, dict) + # If v is a list, tuple, sequence, dict etc., return the length of the container filtering out any empty sun elements, + else ( + len(list(filterfalse(partial(is_empty), v))) + if isinstance(v, Container) + else 1 + ) + ) for v in d.values() ) ) @@ -420,6 +435,10 @@ class Decision(Enum): above_max_chunks = "above_max_chunks" #: Fewer signal chunks for this read collected than required below_min_chunks = "below_min_chunks" + #: Potential second half of a duplex read + duplex_override = "duplex_override" + #: Read sequenced as translocated portion was of unknown length at start of readfish + first_read_override = "first_read_override" @unique @@ -792,3 +811,111 @@ def get_action(self, channel: int) -> Optional[Action]: :return: The last action sent for the channel, or None if no action has been sent. """ return self.last_actions.get(channel, None) + + +@attrs.define +class DuplexTracker: + """ + Wrapper class to keep track the alignment location of the latest read seen on a channel, + and previous decision made, tracking whether we made a duplex override on the last read + Specifically, we store a list of tuples of any target contig names and strands that were aligned to, + keyed to channel number and the previous decision for a read made on that channel. + The decision should only be updated when a read has been finalised and should not be seen again, + i.e a Stop receiving or Unblock has been sent to MinKNOW + No maps are specified as (*, *) + """ + + # Note - `readfish.src.plugins.utils.Alignment` could be used here, instead of tuple[str, Strand] + # We could then use the results of the ALignment directly, if we ever wanted to do something more complex + # for duplex + previous_alignments: Dict[int, list[tuple[str, Strand]]] = attrs.Factory(dict) + previous_decision: Dict[int, Decision] = attrs.Factory(dict) + + def get_previous_decision(self, channel: int) -> Decision: + """ + Get the previous decision seen on this channel. + + :param channel: The channel number. + :return: Previously seen decision + >>> dt = DuplexTracker() + >>> dt.get_previous_decision(1) is None + True + >>> dt.set_decision(1, Decision.duplex_override) + >>> dt.get_previous_decision(1) + + """ + return self.previous_decision.get(channel, None) + + def set_decision(self, channel: int, decision: Decision) -> None: + """ + Set the previous decision for a given channel number. + + :param channel: The channel number. + :param decision: The decision taken. Should be the final decision, + i.e we won't see the read again. + >>> dt = DuplexTracker() + >>> dt.set_decision(1, Decision.no_map) + >>> dt.previous_decision[1] + + """ + self.previous_decision[channel] = decision + + def get_previous_alignments(self, channel: int) -> list[tuple[str, Strand]]: + """ + Retrieves last alignments, including no maps seen on the given channel. + + :param channel: The channel number to lookup the previous action for + :param read_id: Read of ID of the current alignment + :return: Returns a tuple of (contig_name, strand), for the last alignment seen on this channel + + >>> dt = DuplexTracker() + >>> dt.get_previous_alignments(1) is None + True + >>> dt.set_alignments(1, [("contig1", Strand.forward), ("contig2", Strand.reverse)]) + >>> dt.get_previous_alignments(1) + [('contig1', ), ('contig2', )] + """ + return self.previous_alignments.get(channel, None) + + def set_alignments( + self, channel: int, alignments: list[tuple[str, Strand]] + ) -> None: + """ + Add an alignment that has been seen for a channel. + + :param channel: The channel number to set the alignment for. + :param target_name: The name of the target contig aligned to + :param strand: The strand we have aligned to. + + >>> dt = DuplexTracker() + >>> dt.set_alignments(1, [("contig3", Strand.forward), ("contig4", Strand.reverse)]) + >>> dt.previous_alignments[1] + [('contig3', ), ('contig4', )] + """ + self.previous_alignments[channel] = alignments + + def possible_duplex( + self, channel: int, target_name: str, strand: Strand | str | int + ) -> bool: + """ + Compare the current alignment target_name and strand for a given channel + with the previous alignment target_name and strand. + + If the strand is opposite and the target is the same, return True, else False. + :param channel: Channel number to fetch alignment for + :param target_name: The name of the target contig for the current alignment + :param strand: The strand of the current alignment + :return: True if the strand is opposite and target contig the same + + >>> dt = DuplexTracker() + >>> dt.set_alignments(1, [("contig5", Strand.forward)]) + >>> dt.possible_duplex(1, "contig5", Strand.reverse) + True + >>> dt.possible_duplex(1, "contig6", Strand.reverse) + False + """ + strand = Strand(strand) + return any( + prev_alignment == (target_name, ~strand) + for prev_alignment in self.get_previous_alignments(channel) + ) diff --git a/tests/static/mappy_validation_test/fail/001_wrong_reference_type.txt b/tests/static/mappy_validation_test/fail/001_wrong_reference_type.txt index 0193d70c..52f62746 100644 --- a/tests/static/mappy_validation_test/fail/001_wrong_reference_type.txt +++ b/tests/static/mappy_validation_test/fail/001_wrong_reference_type.txt @@ -1 +1 @@ -Provided index file appears to be of an incorrect type - should be one of ['.fasta', '.fna', '.fsa', '.fa', '.fastq', '.fq', '.fasta.gz', '.fna.gz', '.fsa.gz', '.fa.gz', '.fastq.gz', '.fq.gz', '.mmi'] +Provided index file appears to be of an incorrect type - should be one of ['.fasta', '.fna', '.fsa', '.fa', '.fastq', '.fq', '.fasta.gz', '.fna.gz', '.fsa.gz', '.fa.gz', '.fastq.gz', '.fq.gz', '.mmi'] \ No newline at end of file diff --git a/tests/static/mappy_validation_test/fail/003_wrong_reference_type_mmi_gz.txt b/tests/static/mappy_validation_test/fail/003_wrong_reference_type_mmi_gz.txt index 0193d70c..52f62746 100644 --- a/tests/static/mappy_validation_test/fail/003_wrong_reference_type_mmi_gz.txt +++ b/tests/static/mappy_validation_test/fail/003_wrong_reference_type_mmi_gz.txt @@ -1 +1 @@ -Provided index file appears to be of an incorrect type - should be one of ['.fasta', '.fna', '.fsa', '.fa', '.fastq', '.fq', '.fasta.gz', '.fna.gz', '.fsa.gz', '.fa.gz', '.fastq.gz', '.fq.gz', '.mmi'] +Provided index file appears to be of an incorrect type - should be one of ['.fasta', '.fna', '.fsa', '.fa', '.fastq', '.fq', '.fasta.gz', '.fna.gz', '.fsa.gz', '.fa.gz', '.fastq.gz', '.fq.gz', '.mmi'] \ No newline at end of file