From 2dd502b805e12ea7e3e8f2122038a821cf35495d Mon Sep 17 00:00:00 2001 From: binho Date: Tue, 2 Apr 2024 13:18:49 -0500 Subject: [PATCH] Add code for paleogeo event annotation, update unit tests --- .../functionality/event_series.py | 183 ++++++++++++++---- src/phylojunction/functionality/feature_io.py | 58 +++++- src/phylojunction/functionality/stoch_map.py | 25 +++ .../test_event_series_four_regions.py | 112 +++++------ .../test_event_series_two_regions.py | 45 ++++- 5 files changed, 321 insertions(+), 102 deletions(-) diff --git a/src/phylojunction/functionality/event_series.py b/src/phylojunction/functionality/event_series.py index 56f4932..c6e6e0e 100644 --- a/src/phylojunction/functionality/event_series.py +++ b/src/phylojunction/functionality/event_series.py @@ -457,6 +457,12 @@ def __init__(self, if verbose: print(" ... finished reading paleogeographic events.") + # annotate paleogeographic events + self.annotate_paleogeo_events() + + if verbose: + print(" ... finished annotating paleogeographic events.") + self.populate_trunc_event_series() if verbose: @@ -803,7 +809,7 @@ def parse_range_contraction(ch1_set, # if range is stable before the contraction and unstable # after, then this is 'speciation by extinction' - def do_anagenetic_smap(ch1_set, ch2_set, ana_smap): + def do_anagenetic_smap(ch1_set, ch2_set, ana_smap) -> None: # gathering initial information to annotate # anagenetic stochastic map @@ -844,11 +850,11 @@ def do_anagenetic_smap(ch1_set, ch2_set, ana_smap): ana_conn_graph) def recursively_populate_event_series_dict(nd: dp.Node, - it_idx: int): + it_idx: int) -> None: """Populate self._event_series_dict recursively. - Disambiguation of range expansion (dispersal) events - happens here. + Leaf nodes are ignored. Disambiguation of range expansion + (dispersal) events happens here. """ # do current node @@ -897,18 +903,20 @@ def recursively_populate_event_series_dict(nd: dp.Node, conn_graph_list[clado_smap_time_slice_idx] # child 1 - ch1_bp = clado_smap.to_state_bit_patt + # ch1_bp = clado_smap.to_state_bit_patt # child 2 (will be None if no range split at speciation) - ch2_bp = clado_smap.to_state2_bit_patt + # ch2_bp = clado_smap.to_state2_bit_patt # get sets of region indices for the two mutually # exclusive ranges - ch1_set = \ - set([idx for idx, b in enumerate(ch1_bp) if b == "1"]) + ch1_set = clado_smap.to_state_idx_set + # ch1_set = \ + # set([idx for idx, b in enumerate(ch1_bp) if b == "1"]) - ch2_set = ch1_set - if ch2_bp is not None: - ch2_set = \ - set([idx for idx, b in enumerate(ch2_bp) if b == "1"]) + ch2_set = clado_smap.to_state2_idx_set + # ch2_set = ch1_set + # if ch2_bp is not None: + # ch2_set = \ + # set([idx for idx, b in enumerate(ch2_bp) if b == "1"]) # only speciation events with range splitting # if clado_smap.to_state2_bit_patt is not None: @@ -969,16 +977,17 @@ def recursively_populate_event_series_dict(nd: dp.Node, # iterating over each MCMC iteration when stochastic maps were logged for it_idx, smap in self._smap_collection.stoch_maps_tree_dict.items(): - ann_tr = smap.ann_tr seed_nd = smap.ann_tr.origin_node if smap.ann_tr.with_origin \ else smap.ann_tr.root_node # populate self._self._event_series_dict + # (terminal nodes are ignored) # # key1: node name - # value1 iteration dict: - # key2: iteration index - # value2 event series object + # value1 iteration dict + # key2: iteration index + # value2 event series object (or empty dictionary if internal + # node but no events, like origin node) recursively_populate_event_series_dict(seed_nd, it_idx) @@ -1036,9 +1045,9 @@ def insort_barrier(from_region_idx, for nd_label, it_event_series_dict in self._event_series_dict.items(): for it_idx, event_series in it_event_series_dict.items(): - # event_series will be an empty dictionary if range - # did not split at node -- we do not care about those! - # (even if subtending branch has maps!) + # node at the start of process may have an empty event series + # if stoch mapping file for some reason did not have an entry + # those nodes if isinstance(event_series, EvolRelevantEventSeries): event_list = event_series.event_list @@ -1080,6 +1089,7 @@ def insort_barrier(from_region_idx, else range((from_region_idx+1), n_regions) for to_region_idx in inner_loop_idxs: + # ignore diagonal elements if directed and from_region_idx == to_region_idx: continue @@ -1103,6 +1113,109 @@ def insort_barrier(from_region_idx, to_region_idx, event_list) + def annotate_paleogeo_events(self) -> None: + """Annotate paleogeo events as (de)stabilizing or not. + + After paleogeographic events have been added to the event series + of each internal node, we visit each one of them and annotate + them according to their (de)stabilizing status with respect to + the range splitting event happening at speciation. If the range + does not split at speciation, the annotation is 'N/A'. + """ + + # only internal nodes in _event_series_dict + for nd_label, it_event_series_dict in self._event_series_dict.items(): + for it_idx, event_series in it_event_series_dict.items(): + # node at the start of process may have an empty event series + # if stoch mapping file for some reason did not have an entry + # those nodes + if isinstance(event_series, EvolRelevantEventSeries): + event_list = event_series.event_list + n_events = len(event_list) + range_split_smap = event_list[-1] + assert (isinstance(range_split_smap, pjsmap.RangeSplitOrBirth)) + + # gathering splitting range info + ch1_set = range_split_smap.to_state_idx_set + ch2_set = range_split_smap.to_state2_idx_set + + for ev_idx, ev in enumerate(event_list): + if isinstance(ev, (pjfio.BarrierAppearance, + pjfio.BarrierDisappearance)): + + # gathering paleogeo event info + # + # pair of regions dis/reconnected by barrier + # appearing or disappearing + reg1_idx = ev.from_node_idx + reg2_idx = ev.to_node_idx + + # if the range splits at the internal node, then we + # need to do things + # + # grab range bit pattern at paleogeographic event + range_bp_at_paleogeo_ev = "" + + for next_ev_idx in range((ev_idx + 1), n_events): + next_ev = event_list[next_ev_idx] + + if isinstance(next_ev, pjsmap.StochMap): + range_bp_at_paleogeo_ev = \ + next_ev.from_state_bit_patt + # now get region indices of occupied range + range_bp_at_paleogeo_set = \ + set([idx for idx, b in \ + enumerate(range_bp_at_paleogeo_ev) \ + if b == "1"]) + break + + assert (range_bp_at_paleogeo_ev != "") + + # annotate paleogeo event with occupied range + ev.range_bit_patt = range_bp_at_paleogeo_ev + + ################################################## + # Scenario 1 in which (de)stabilization does not # + # make sense: if there is no range split # + ################################################## + if ch1_set == ch2_set: + if isinstance(ev, pjfio.BarrierDisappearance): + ev.restabilized_range = "N/A" + + elif isinstance(ev, pjfio.BarrierAppearance): + ev.destabilized_range = "N/A" + + continue + + ################################################## + # Scenario 2 in which (de)stabilization does not # + # make sense: when the regions (dis)reconnected # + # are not occupied by lineage # + ################################################## + if reg1_idx not in range_bp_at_paleogeo_set or \ + reg2_idx not in range_bp_at_paleogeo_set: + if isinstance(ev, pjfio.BarrierDisappearance): + ev.restabilized_range = "N/A" + + elif isinstance(ev, pjfio.BarrierAppearance): + ev.destabilized_range = "N/A" + + continue + + # if paleogeo eventis split-relevant + if (reg1_idx in ch1_set and reg2_idx in ch2_set) or \ + (reg1_idx in ch2_set and reg2_idx in ch1_set) and \ + ch1_set != ch2_set: + + # if barrier disappearing is split-relevant + # (if split is within-region, we will never get inside + # this if block) + if isinstance(ev, pjfio.BarrierDisappearance): + ev.restabilized_range = "stab" + + elif isinstance(ev, pjfio.BarrierAppearance): + ev.destabilized_range = "destab" + def populate_trunc_event_series(self) -> None: """Populate _truncated_event_series_dict @@ -1133,12 +1246,14 @@ def populate_trunc_event_series(self) -> None: assert (isinstance(range_split_smap, pjsmap.RangeSplitOrBirth)) # gathering splitting range info - ch1_bp = range_split_smap.to_state_bit_patt - ch2_bp = range_split_smap.to_state2_bit_patt - ch1_set = set([idx for idx, b in enumerate(ch1_bp) if b == "1"]) - ch2_set = ch1_set - if ch2_bp is not None: - ch2_set = set([idx for idx, b in enumerate(ch2_bp) if b == "1"]) + ch1_set = range_split_smap.to_state_idx_set + ch2_set = range_split_smap.to_state2_idx_set + # ch1_bp = range_split_smap.to_state_bit_patt + # ch2_bp = range_split_smap.to_state2_bit_patt + # ch1_set = set([idx for idx, b in enumerate(ch1_bp) if b == "1"]) + # ch2_set = ch1_set + # if ch2_bp is not None: + # ch2_set = set([idx for idx, b in enumerate(ch2_bp) if b == "1"]) # now preparing truncated event list # @@ -1163,7 +1278,7 @@ def populate_trunc_event_series(self) -> None: ############################################ if isinstance(ev, pjsmap.RangeExpansion) and \ ev.stabilized_range_wrt_split and \ - ch2_bp is not None: + ch1_set != ch2_set: truncate_here = True ############################################# @@ -1173,18 +1288,20 @@ def populate_trunc_event_series(self) -> None: # regions are on opposing sides of the # # split # ############################################# - elif isinstance(ev, pjfio.BarrierDisappearance): + elif isinstance(ev, pjfio.BarrierDisappearance) and \ + ev.restabilized_range == "stab" and \ + ch1_set != ch2_set: # pair of regions reconnected by barrier disappearing - reg1_idx = ev.from_node_idx - reg2_idx = ev.to_node_idx + # reg1_idx = ev.from_node_idx + # reg2_idx = ev.to_node_idx # if barrier disappearing is split-relevant # (if split is within-region, we will never get inside # this if block) - if (reg1_idx in ch1_set and reg2_idx in ch2_set) or \ - (reg1_idx in ch2_set and reg2_idx in ch1_set) and \ - ch2_bp is not None: - truncate_here = True + # if (reg1_idx in ch1_set and reg2_idx in ch2_set) or \ + # (reg1_idx in ch2_set and reg2_idx in ch1_set) and \ + # ch2_bp is not None: + truncate_here = True # actually truncate! if truncate_here: diff --git a/src/phylojunction/functionality/feature_io.py b/src/phylojunction/functionality/feature_io.py index 640fa9c..95df943 100644 --- a/src/phylojunction/functionality/feature_io.py +++ b/src/phylojunction/functionality/feature_io.py @@ -251,6 +251,11 @@ class BarrierAppearance(pjev.EvolRelevantEvent): in barrier. """ + from_node_idx: int + to_node_idx: int + destabilized_range: ty.Optional[str] + range_bit_patt: ty.Optional[str] + def __init__(self, n_chars: int, age: ty.Optional[float], @@ -263,17 +268,37 @@ def __init__(self, self.from_node_idx = from_node_idx self.to_node_idx = to_node_idx + self.destabilized_range = None + self.range_bit_patt = None def short_str(self) -> str: - short_str = "b+(" + str(self.age) + ")" + short_str = "b+(" + str(self.age) + ")" \ + + "_reg(" + str(self.from_node_idx) + "|" \ + + str(self.to_node_idx) + ")" + + if self.destabilized_range is not None and \ + self.range_bit_patt is not None: + short_str += "_" + self.destabilized_range \ + + "(" + self.range_bit_patt + ")" return short_str def __str__(self) -> str: - return "Barrier (at age = " + str(self.age) \ + str_representation = \ + "Barrier (at age = " + str(self.age) \ + ") between region " + str(self.from_node_idx) \ + " and " + str(self.to_node_idx) + if self.destabilized_range is not None and \ + self.range_bit_patt is not None: + str1 = "\n Destabilized range (" + self.range_bit_patt + "): " + str2 = "True" if self.destabilized_range == "destab" \ + else "N/A" + suffix = str1 + str2 + str_representation += suffix + + return str_representation + def __lt__(self, other) -> bool: return super().__lt__(other) @@ -288,6 +313,11 @@ class BarrierDisappearance(pjev.EvolRelevantEvent): in barrier. """ + from_node_idx: int + to_node_idx: int + restabilized_range: ty.Optional[str] + range_bit_patt: ty.Optional[str] # range at time of event + def __init__(self, n_chars: int, age: ty.Optional[float], @@ -300,18 +330,38 @@ def __init__(self, self.from_node_idx = from_node_idx self.to_node_idx = to_node_idx + self.restabilized_range = None + self.range_bit_patt = None def short_str(self) -> str: - short_str = "b-(" + str(self.age) + ")" + short_str = "b-(" + str(self.age) + ")" \ + + "_reg(" + str(self.from_node_idx) + "|" \ + + str(self.to_node_idx) + ")" + + if self.restabilized_range is not None and \ + self.range_bit_patt is not None: + short_str += "_" + self.restabilized_range \ + + "(" + self.range_bit_patt + ")" return short_str def __str__(self) -> str: - return "Connectivity restablished (at age = " \ + str_representation = \ + "Connectivity restablished (at age = " \ + str(self.age) + ") between region " \ + str(self.from_node_idx) + " and " \ + str(self.to_node_idx) + if self.restabilized_range is not None and \ + self.range_bit_patt is not None: + str1 = "\n Restabilized range (" + self.range_bit_patt + "): " + str2 = "True" if self.restabilized_range == "stab" \ + else "N/A" + suffix = str1 + str2 + str_representation += suffix + + return str_representation + def __lt__(self, other) -> bool: return super().__lt__(other) diff --git a/src/phylojunction/functionality/stoch_map.py b/src/phylojunction/functionality/stoch_map.py index 9235b17..d0180e9 100644 --- a/src/phylojunction/functionality/stoch_map.py +++ b/src/phylojunction/functionality/stoch_map.py @@ -27,12 +27,16 @@ class StochMap(pjev.EvolRelevantEvent): target biogeographic range). to_state_bit_pattern (str): Bit pattern representation of target state (i.e., target biogeographic range). + to_state_idx_set (set): Set of indices corresponding to the + position of '1's in 'to_state_bit_patt', for child 1. to_state2 (int): Integer representation of second target state if stochastic map is cladogenetic (i.e., second target biogeographic range). to_state2_bit_pattern (str): Bit pattern representation of second target state if stochastic map is cladogenetic (i.e., second target biogeographic range). + to_state2_idx_set (set): Set of indices corresponding to the + position of '1's in 'to_state_bit_patt', for child 2. focal_node_name (str): Name of the node subtending the branch along which the stochastic map was placed. parent_node_name (str): Name of the node parent to the focal @@ -49,8 +53,10 @@ class StochMap(pjev.EvolRelevantEvent): from_state_bit_patt: str to_state: int to_state_bit_patt: str + to_state_idx_set: ty.Set[int] to_state2: int to_state2_bit_patt: str + to_state2_idx_set: ty.Set[int] focal_node_name: str parent_node_name: str child_node_name: str @@ -69,8 +75,10 @@ def __init__(self, child_node_name: str, map_type: ty.Optional[str] = None, time: ty.Optional[int] = None, + to_state_idx_set: ty.Optional[ty.Set[int]] = None, to_state2: ty.Optional[int] = None, to_state2_bit_patt: ty.Optional[str] = None, + to_state2_idx_set: ty.Optional[ty.Set[int]] = None, child2_node_name: ty.Optional[str] = None) -> None: # initializing EvolRelevantEvent members @@ -81,6 +89,7 @@ def __init__(self, self.from_state_bit_patt = from_state_bit_patt self.to_state = to_state self.to_state_bit_patt = to_state_bit_patt + self.to_state_idx_set = to_state_idx_set self.focal_node_name = focal_node_name self.parent_node_name = parent_node_name self.child_node_name = child_node_name @@ -89,6 +98,7 @@ def __init__(self, self.child2_node_name = child2_node_name self.to_state2 = to_state2 self.to_state2_bit_patt = to_state2_bit_patt + self.to_state2_idx_set = to_state2_idx_set def short_str(self) -> str: super().short_str() @@ -520,6 +530,16 @@ def __init__(self, to_state2_bit_patt: ty.Optional[str] = None, time: ty.Optional[float] = None) -> None: + to_state_idx_set = \ + set([idx for idx, b in \ + enumerate(to_state_bit_patt) if b == "1"]) + + to_state2_idx_set = to_state_idx_set + if to_state2_bit_patt is not None: + to_state2_idx_set = \ + set([idx for idx, b in \ + enumerate(to_state2_bit_patt) if b == "1"]) + super().__init__(n_regions, age, from_state, @@ -529,9 +549,11 @@ def __init__(self, focal_node_name, parent_node_name, child_node_name, + to_state_idx_set=to_state_idx_set, child2_node_name=child2_node_name, to_state2=to_state2, to_state2_bit_patt=to_state2_bit_patt, + to_state2_idx_set=to_state2_idx_set, map_type="range_split", time=time) @@ -550,6 +572,9 @@ def add_child2_info(self, self.to_state2 = to_state2 self.to_state2_bit_patt = to_state2_bit_patt + self.to_state2_idx_set = \ + set([idx for idx, b in enumerate(to_state2_bit_patt) \ + if b == "1"]) self.size_of_child2_range = \ sum(int(i) for i in to_state2_bit_patt) self._update_str_representation_child2() diff --git a/tests/functionality/test_event_series_four_regions.py b/tests/functionality/test_event_series_four_regions.py index 8cda459..3f4a1ad 100644 --- a/tests/functionality/test_event_series_four_regions.py +++ b/tests/functionality/test_event_series_four_regions.py @@ -89,21 +89,22 @@ def setUpClass(cls) -> None: def test_event_list_4reg_it1(self) -> None: """Test four-region, 4-taxon tree event list, iteration 1.""" - nd7_exp = ("b+(5.0)" - "b+(4.0)" + nd7_exp = ("b+(5.0)_reg(0|2)_N/A(1001)" + "b+(4.0)_reg(1|3)_N/A(1001)" "s(3.25)_st(1001>1001_1001)") - nd5_exp = ("b+(5.0)" - "b+(4.0)" + nd5_exp = ("b+(5.0)_reg(0|2)_N/A(1001)" + "b+(4.0)_reg(1|3)_N/A(1001)" "s(3.25)_st(1001>1001_1001)" - "b+(3.0)" + "b+(3.0)_reg(0|3)_destab(1001)" "d(2.5)_si|nb|wc|un(1001)>un(1101)" - "b+(2.0)" + "b+(2.0)_reg(1|2)_N/A(1101)" "d(0.75)_si|nb|wc|un(1101)>un(1111)" "s(0.25)_un(1111>1100_0011)") - nd6_exp = ("b+(5.0)""b+(4.0)" + nd6_exp = ("b+(5.0)_reg(0|2)_N/A(1001)" + "b+(4.0)_reg(1|3)_N/A(1001)" "s(3.25)_st(1001>1001_1001)" - "b+(3.0)" - "b+(2.0)" + "b+(3.0)_reg(0|3)_destab(1001)" + "b+(2.0)_reg(1|2)_N/A(1001)" "d(1.5)_sr|ob|oc|un(1001)>un(1011)" "d(0.75)_si|nb|wc|un(1011)>un(1111)" "s(0.25)_un(1111>1100_0011)") @@ -136,24 +137,24 @@ def test_event_list_4reg_it1(self) -> None: def test_event_list_4reg_it2(self) -> None: """Test four-region, 4-taxon tree event list, iteration 2.""" - nd7_exp = ("b+(5.0)" - "b+(4.0)" + nd7_exp = ("b+(5.0)_reg(0|2)_N/A(1101)" + "b+(4.0)_reg(1|3)_N/A(1101)" "s(3.25)_st(1101>1101_1101)") - nd5_exp = ("b+(5.0)" - "b+(4.0)" + nd5_exp = ("b+(5.0)_reg(0|2)_N/A(1101)" + "b+(4.0)_reg(1|3)_destab(1101)" "s(3.25)_st(1101>1101_1101)" - "b+(3.0)" + "b+(3.0)_reg(0|3)_destab(1101)" "e(2.5)_|st(1101)>un(1100)" - "b+(2.0)" + "b+(2.0)_reg(1|2)_N/A(1100)" "d(1.5)_sr|ob|oc|st(1100)>un(1110)" "d(0.75)_si|nb|wc|un(1110)>un(1111)" "s(0.25)_un(1111>1100_0011)") - nd6_exp = ("b+(5.0)" - "b+(4.0)" + nd6_exp = ("b+(5.0)_reg(0|2)_N/A(1101)" + "b+(4.0)_reg(1|3)_destab(1101)" "s(3.25)_st(1101>1101_1101)" - "b+(3.0)" + "b+(3.0)_reg(0|3)_destab(1101)" "e(2.5)_|st(1101)>st(1001)" - "b+(2.0)" + "b+(2.0)_reg(1|2)_N/A(1001)" "d(1.5)_si|nb|wc|un(1001)>un(1011)" "d(0.75)_si|nb|wc|un(1011)>un(1111)" "s(0.25)_un(1111>1100_0011)") @@ -186,25 +187,25 @@ def test_event_list_4reg_it2(self) -> None: def test_event_list_4reg_it3(self) -> None: """Test four-region, 4-taxon tree event list, iteration 3.""" - nd7_exp = ("b+(5.0)" + nd7_exp = ("b+(5.0)_reg(0|2)_N/A(1000)" "d(4.5)_si|ob|wc|st(1000)>st(1010)" - "b+(4.0)" + "b+(4.0)_reg(1|3)_N/A(1010)" "s(3.25)_st(1010>1010_1010)") - nd5_exp = ("b+(5.0)" + nd5_exp = ("b+(5.0)_reg(0|2)_N/A(1000)" "d(4.5)_sr|ob|wc|st(1000)>un(1010)" - "b+(4.0)" + "b+(4.0)_reg(1|3)_N/A(1010)" "s(3.25)_st(1010>1010_1010)" - "b+(3.0)" + "b+(3.0)_reg(0|3)_N/A(1010)" "d(2.5)_si|nb|wc|un(1010)>st(1110)" - "b+(2.0)" + "b+(2.0)_reg(1|2)_destab(1110)" "d(0.75)_si|nb|wc|un(1110)>un(1111)" "s(0.25)_un(1111>1100_0011)") - nd6_exp = ("b+(5.0)" + nd6_exp = ("b+(5.0)_reg(0|2)_N/A(1000)" "d(4.5)_sr|ob|wc|st(1000)>un(1010)" - "b+(4.0)" + "b+(4.0)_reg(1|3)_N/A(1010)" "s(3.25)_st(1010>1010_1010)" - "b+(3.0)" - "b+(2.0)" + "b+(3.0)_reg(0|3)_N/A(1010)" + "b+(2.0)_reg(1|2)_N/A(1010)" "d(1.5)_si|nb|wc|un(1010)>un(1110)" "d(0.75)_si|nb|wc|un(1110)>un(1111)" "s(0.25)_un(1111>1100_0011)") @@ -237,21 +238,22 @@ def test_event_list_4reg_it3(self) -> None: def test_trunc_event_list_4reg_it1(self) -> None: """Test four-region, 4-taxon tree truncated event list, iteration 1.""" - nd7_exp = ("b+(5.0)" - "b+(4.0)" + nd7_exp = ("b+(5.0)_reg(0|2)_N/A(1001)" + "b+(4.0)_reg(1|3)_N/A(1001)" "s(3.25)_st(1001>1001_1001)") - nd5_exp = ("b+(5.0)" - "b+(4.0)" + nd5_exp = ("b+(5.0)_reg(0|2)_N/A(1001)" + "b+(4.0)_reg(1|3)_N/A(1001)" "s(3.25)_st(1001>1001_1001)" - "b+(3.0)" + "b+(3.0)_reg(0|3)_destab(1001)" "d(2.5)_si|nb|wc|un(1001)>un(1101)" - "b+(2.0)" + "b+(2.0)_reg(1|2)_N/A(1101)" "d(0.75)_si|nb|wc|un(1101)>un(1111)" "s(0.25)_un(1111>1100_0011)") - nd6_exp = ("b+(5.0)""b+(4.0)" + nd6_exp = ("b+(5.0)_reg(0|2)_N/A(1001)" + "b+(4.0)_reg(1|3)_N/A(1001)" "s(3.25)_st(1001>1001_1001)" - "b+(3.0)" - "b+(2.0)" + "b+(3.0)_reg(0|3)_destab(1001)" + "b+(2.0)_reg(1|2)_N/A(1001)" "d(1.5)_sr|ob|oc|un(1001)>un(1011)" "d(0.75)_si|nb|wc|un(1011)>un(1111)" "s(0.25)_un(1111>1100_0011)") @@ -284,24 +286,24 @@ def test_trunc_event_list_4reg_it1(self) -> None: def test_trunc_event_list_4reg_it2(self) -> None: """Test four-region, 4-taxon tree truncated event list, iteration 2.""" - nd7_exp = ("b+(5.0)" - "b+(4.0)" + nd7_exp = ("b+(5.0)_reg(0|2)_N/A(1101)" + "b+(4.0)_reg(1|3)_N/A(1101)" "s(3.25)_st(1101>1101_1101)") - nd5_exp = ("b+(5.0)" - "b+(4.0)" + nd5_exp = ("b+(5.0)_reg(0|2)_N/A(1101)" + "b+(4.0)_reg(1|3)_destab(1101)" "s(3.25)_st(1101>1101_1101)" - "b+(3.0)" + "b+(3.0)_reg(0|3)_destab(1101)" "e(2.5)_|st(1101)>un(1100)" - "b+(2.0)" + "b+(2.0)_reg(1|2)_N/A(1100)" "d(1.5)_sr|ob|oc|st(1100)>un(1110)" "d(0.75)_si|nb|wc|un(1110)>un(1111)" "s(0.25)_un(1111>1100_0011)") - nd6_exp = ("b+(5.0)" - "b+(4.0)" + nd6_exp = ("b+(5.0)_reg(0|2)_N/A(1101)" + "b+(4.0)_reg(1|3)_destab(1101)" "s(3.25)_st(1101>1101_1101)" - "b+(3.0)" + "b+(3.0)_reg(0|3)_destab(1101)" "e(2.5)_|st(1101)>st(1001)" - "b+(2.0)" + "b+(2.0)_reg(1|2)_N/A(1001)" "d(1.5)_si|nb|wc|un(1001)>un(1011)" "d(0.75)_si|nb|wc|un(1011)>un(1111)" "s(0.25)_un(1111>1100_0011)") @@ -334,19 +336,19 @@ def test_trunc_event_list_4reg_it2(self) -> None: def test_trunc_event_list_4reg_it3(self) -> None: """Test four-region, 4-taxon tree truncated event list, iteration 3.""" - nd7_exp = ("b+(5.0)" + nd7_exp = ("b+(5.0)_reg(0|2)_N/A(1000)" "d(4.5)_si|ob|wc|st(1000)>st(1010)" - "b+(4.0)" + "b+(4.0)_reg(1|3)_N/A(1010)" "s(3.25)_st(1010>1010_1010)") - nd5_exp = ("b+(2.0)" + nd5_exp = ("b+(2.0)_reg(1|2)_destab(1110)" "d(0.75)_si|nb|wc|un(1110)>un(1111)" "s(0.25)_un(1111>1100_0011)") - nd6_exp = ("b+(5.0)" + nd6_exp = ("b+(5.0)_reg(0|2)_N/A(1000)" "d(4.5)_sr|ob|wc|st(1000)>un(1010)" - "b+(4.0)" + "b+(4.0)_reg(1|3)_N/A(1010)" "s(3.25)_st(1010>1010_1010)" - "b+(3.0)" - "b+(2.0)" + "b+(3.0)_reg(0|3)_N/A(1010)" + "b+(2.0)_reg(1|2)_N/A(1010)" "d(1.5)_si|nb|wc|un(1010)>un(1110)" "d(0.75)_si|nb|wc|un(1110)>un(1111)" "s(0.25)_un(1111>1100_0011)") diff --git a/tests/functionality/test_event_series_two_regions.py b/tests/functionality/test_event_series_two_regions.py index c106ff5..757ccd5 100644 --- a/tests/functionality/test_event_series_two_regions.py +++ b/tests/functionality/test_event_series_two_regions.py @@ -89,22 +89,24 @@ def setUpClass(cls) -> None: def test_event_list_it1(self) -> None: """Test two-region, 5-taxon tree event list.""" + nd9_exp = ("s(5.0)_st(10>10_10)") nd6_exp = ("s(5.0)_st(10>10_10)" "d(4.75)_sr|nb|wc|st(10)>st(11)" - "b+(4.25)" + "b+(4.25)_reg(0|1)_destab(11)" "s(4.0)_un(11>10_01)") nd8_exp = ("s(5.0)_st(10>10_10)" - "b+(4.25)" + "b+(4.25)_reg(0|1)_N/A(10)" "d(3.25)_sr|ob|oc|st(10)>un(11)" "s(3.0)_un(11>10_01)") nd7_exp = ("s(5.0)_st(10>10_10)" - "b+(4.25)" + "b+(4.25)_reg(0|1)_N/A(10)" "d(3.25)_sr|ob|oc|st(10)>un(11)" "s(3.0)_un(11>10_01)" "d(2.75)_sr|ob|oc|st(01)>un(11)" - "b-(1.5)" + "b-(1.5)_reg(0|1)_stab(11)" "s(1.0)_st(11>10_01)") + nd9_obs = str() nd6_obs = str() nd8_obs = str() nd7_obs = str() @@ -117,7 +119,10 @@ def test_event_list_it1(self) -> None: if isinstance(event_series, pjev.EvolRelevantEventSeries): if it_idx in it_to_look_at: for ev in event_series.event_list: - if nd_label == "nd6": + if nd_label == "nd9": + nd9_obs += ev.short_str() + + elif nd_label == "nd6": nd6_obs += ev.short_str() elif nd_label == "nd7": @@ -126,6 +131,7 @@ def test_event_list_it1(self) -> None: elif nd_label == "nd8": nd8_obs += ev.short_str() + self.assertEqual(nd9_exp, nd9_obs) self.assertEqual(nd6_exp, nd6_obs) self.assertEqual(nd7_exp, nd7_obs) self.assertEqual(nd8_exp, nd8_obs) @@ -133,14 +139,16 @@ def test_event_list_it1(self) -> None: def test_trunc_event_list_it1(self) -> None: """Test two-region, 5-taxon tree truncated event list.""" - nd6_exp = ("b+(4.25)" + nd9_exp = "s(5.0)_st(10>10_10)" + nd6_exp = ("b+(4.25)_reg(0|1)_destab(11)" "s(4.0)_un(11>10_01)") nd8_exp = ("s(5.0)_st(10>10_10)" - "b+(4.25)" + "b+(4.25)_reg(0|1)_N/A(10)" "d(3.25)_sr|ob|oc|st(10)>un(11)" "s(3.0)_un(11>10_01)") nd7_exp = "s(1.0)_st(11>10_01)" + nd9_obs = str() nd6_obs = str() nd8_obs = str() nd7_obs = str() @@ -153,7 +161,10 @@ def test_trunc_event_list_it1(self) -> None: if isinstance(event_series, pjev.EvolRelevantEventSeries): if it_idx in it_to_look_at: for ev in event_series.trunc_event_list: - if nd_label == "nd6": + if nd_label == "nd9": + nd9_obs += ev.short_str() + + elif nd_label == "nd6": nd6_obs += ev.short_str() elif nd_label == "nd7": @@ -162,6 +173,7 @@ def test_trunc_event_list_it1(self) -> None: elif nd_label == "nd8": nd8_obs += ev.short_str() + self.assertEqual(nd9_exp, nd9_obs) self.assertEqual(nd6_exp, nd6_obs) self.assertEqual(nd8_exp, nd8_obs) self.assertEqual(nd7_exp, nd7_obs) @@ -169,10 +181,12 @@ def test_trunc_event_list_it1(self) -> None: def test_hyp_ann(self) -> None: """Test two-region, 5-taxon tree hypothesis annotation.""" + nd9_exp = str("Within-region speciation") nd6_exp = str("Vicariance") nd8_exp = str("Founder event") nd7_exp = str("Ambiguous") + nd9_obs = str() nd6_obs = str() nd8_obs = str() nd7_obs = str() @@ -182,7 +196,10 @@ def test_hyp_ann(self) -> None: # event_series of root will be an empty dictionary, need # to check if isinstance(event_series, pjev.EvolRelevantEventSeries): - if nd_label == "nd6": + if nd_label == "nd9": + nd9_obs = str(event_series.supported_hyp) + + elif nd_label == "nd6": nd6_obs = str(event_series.supported_hyp) elif nd_label == "nd7": @@ -191,6 +208,7 @@ def test_hyp_ann(self) -> None: elif nd_label == "nd8": nd8_obs = str(event_series.supported_hyp) + self.assertEqual(nd9_exp, nd9_obs) self.assertEqual(nd6_exp, nd6_obs) self.assertEqual(nd8_exp, nd8_obs) self.assertEqual(nd7_exp, nd7_obs) @@ -237,6 +255,12 @@ def test_node_count_supporting_hyp_dict(self) -> None: def test_hyp_support_by_node_dict(self) -> None: """Test two-region, 5-taxon tree, node hypothesis.""" + nd9_exp = \ + {'Vicariance': 0, + 'Founder event': 0, + 'Speciation by extinction': 0, + 'Ambiguous': 0, + 'Within-region speciation': 1} nd6_exp = \ {'Vicariance': 1, 'Founder event': 0, @@ -262,11 +286,12 @@ def test_hyp_support_by_node_dict(self) -> None: 'Ambiguous': 0, 'Within-region speciation': 1} + nd9_obs = self.est.hyp_support_by_node_dict["nd9"] nd6_obs = self.est.hyp_support_by_node_dict["nd6"] nd7_obs = self.est.hyp_support_by_node_dict["nd7"] nd8_obs = self.est.hyp_support_by_node_dict["nd8"] - nd9_obs = self.est.hyp_support_by_node_dict["nd9"] + self.assertEqual(nd9_exp, nd9_obs) self.assertEqual(nd6_exp, nd6_obs) self.assertEqual(nd7_exp, nd7_obs) self.assertEqual(nd8_exp, nd8_obs)