Skip to content

Commit f80f75b

Browse files
ahmigahmig
committed
Handle unset or empty GenBank feature filters
Clarifies feature iteration and filtering, and performs "exclude" filter after the "include" filter
1 parent 5c6cfd5 commit f80f75b

2 files changed

Lines changed: 110 additions & 48 deletions

File tree

‎workflow/scripts/N_S_sites_from_fasta.py‎

Lines changed: 54 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -3,11 +3,12 @@
33
import logging
44
import json
55
import itertools as it
6-
from typing import Dict
6+
from typing import Dict, Iterable
77

88
import pandas as pd
99
from Bio import SeqIO
1010
from Bio.SeqRecord import SeqRecord
11+
from Bio.SeqFeature import SeqFeature
1112
from Bio.Seq import Seq
1213

1314

@@ -68,6 +69,54 @@ def calculate_ns_sites(codons: dict, genetic_code: dict) -> pd.DataFrame:
6869
return pd.DataFrame({"feature": features, "N": N_sites, "S": S_sites})
6970

7071

72+
def iter_features_filtering(features: Iterable[SeqFeature], included: Dict[str, str], excluded: Dict[str, str]) -> Iterable[SeqFeature]:
73+
# No filters
74+
if len(included) == 0 and len(excluded) == 0:
75+
logging.debug("Selecting all features")
76+
return iter(features)
77+
# Only inclusion filter
78+
elif len(included) == 0 and len(excluded) != 0:
79+
logging.debug(f"Selecting features excluding all of {excluded}")
80+
return (
81+
feature for feature in features
82+
if all(
83+
(qualifier_value not in excluded.get(qualifier_key, []))
84+
for qualifier_key in excluded.keys()
85+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
86+
)
87+
)
88+
# Only exclusion filter
89+
elif len(included) != 0 and len(excluded) == 0:
90+
logging.debug(f"Selecting features including any of {included}")
91+
return (
92+
feature for feature in features
93+
if any(
94+
(qualifier_value in included.get(qualifier_key, []))
95+
for qualifier_key in included.keys()
96+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
97+
)
98+
)
99+
# Inclusion then exclusion filter
100+
else:
101+
logging.debug(f"Selecting features including any of {included} and then excluding all of {excluded}")
102+
included_features = (
103+
feature for feature in features
104+
if any(
105+
(qualifier_value in included.get(qualifier_key, []))
106+
for qualifier_key in included.keys()
107+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
108+
)
109+
)
110+
return (
111+
feature for feature in included_features
112+
if all(
113+
(qualifier_value not in excluded.get(qualifier_key, []))
114+
for qualifier_key in excluded.keys()
115+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
116+
)
117+
)
118+
119+
71120
def main():
72121

73122
logging.basicConfig(
@@ -85,29 +134,11 @@ def main():
85134
logging.info("Reading input FASTA")
86135
record = read_monofasta(snakemake.input.fasta)
87136

88-
if len(snakemake.params.features) == 0:
89-
logging.debug("Selecting all features")
90-
feature_iterator = iter(gb.features)
91-
else:
92-
included = snakemake.params.features.get("INCLUDE", {})
93-
excluded = snakemake.params.features.get("EXCLUDE", {})
94-
logging.debug(f"Selecting features including any of {included} and excluding all of {excluded}")
95-
feature_iterator = (
96-
feature for feature in gb.features
97-
if any(
98-
(qualifier_value in included.get(qualifier_key, []))
99-
for qualifier_key in included.keys()
100-
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
101-
) and all(
102-
(qualifier_value not in excluded.get(qualifier_key, []))
103-
for qualifier_key in excluded.keys()
104-
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
105-
)
106-
)
107-
108137
logging.info("Extracting CDS")
109138
coding_records = {}
110-
for feature in feature_iterator:
139+
included = snakemake.params.features.get("INCLUDE", {})
140+
excluded = snakemake.params.features.get("EXCLUDE", {})
141+
for feature in iter_features_filtering(gb.features, included, excluded):
111142
logging.debug(
112143
"Processing SeqFeature: "
113144
f"ID={feature.id} type={feature.type} location={feature.location} "
@@ -121,6 +152,7 @@ def main():
121152
elif identifier in coding_records:
122153
logging.warning(f"Identifier '{identifier}' is already among coding records and will not be replaced by feature at {feature.location}")
123154
else:
155+
logging.debug(f"Adding feature")
124156
coding_records[identifier] = feature.extract(record)
125157

126158
logging.info(f"Splitting {len(coding_records)} records into codons")

‎workflow/scripts/window.py‎

Lines changed: 56 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,62 @@
11
#!/usr/bin/env python3
22

33
import logging
4-
import pandas as pd
54
from typing import NewType, Dict, Iterable, List
65

6+
import pandas as pd
77
from Bio import SeqIO
8-
from Bio.SeqFeature import SeqFeature
9-
from Bio.SeqFeature import CompoundLocation, FeatureLocation
8+
from Bio.SeqFeature import SeqFeature, SimpleLocation, CompoundLocation
109

1110

12-
FeatureIndex = NewType("FeatureIndex", Dict[str, FeatureLocation | CompoundLocation])
11+
FeatureIndex = NewType("FeatureIndex", Dict[str, SimpleLocation | CompoundLocation])
12+
13+
14+
def iter_features_filtering(features: Iterable[SeqFeature], included: Dict[str, str], excluded: Dict[str, str]) -> Iterable[SeqFeature]:
15+
# No filters
16+
if len(included) == 0 and len(excluded) == 0:
17+
logging.debug("Selecting all features")
18+
return iter(features)
19+
# Only inclusion filter
20+
elif len(included) == 0 and len(excluded) != 0:
21+
logging.debug(f"Selecting features excluding all of {excluded}")
22+
return (
23+
feature for feature in features
24+
if all(
25+
(qualifier_value not in excluded.get(qualifier_key, []))
26+
for qualifier_key in excluded.keys()
27+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
28+
)
29+
)
30+
# Only exclusion filter
31+
elif len(included) != 0 and len(excluded) == 0:
32+
logging.debug(f"Selecting features including any of {included}")
33+
return (
34+
feature for feature in features
35+
if any(
36+
(qualifier_value in included.get(qualifier_key, []))
37+
for qualifier_key in included.keys()
38+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
39+
)
40+
)
41+
# Inclusion then exclusion filter
42+
else:
43+
logging.debug(f"Selecting features including any of {included} and then excluding all of {excluded}")
44+
included_features = (
45+
feature for feature in features
46+
if any(
47+
(qualifier_value in included.get(qualifier_key, []))
48+
for qualifier_key in included.keys()
49+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
50+
)
51+
)
52+
return (
53+
feature for feature in included_features
54+
if all(
55+
(qualifier_value not in excluded.get(qualifier_key, []))
56+
for qualifier_key in excluded.keys()
57+
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
58+
)
59+
)
1360

1461

1562
def index_features(features: Iterable[SeqFeature]) -> FeatureIndex:
@@ -69,29 +116,12 @@ def main():
69116
logging.info("Reading GenBank file")
70117
gb = SeqIO.read(snakemake.input.gb, format="gb")
71118

72-
# Build feature iterator given the provided filters
73-
if len(snakemake.params.features) == 0:
74-
logging.debug("Selecting all features")
75-
feature_iterator = iter(gb.features)
76-
else:
77-
included = snakemake.params.features.get("INCLUDE", {})
78-
excluded = snakemake.params.features.get("EXCLUDE", {})
79-
logging.debug(f"Selecting features including any of {included} and excluding all of {excluded}")
80-
feature_iterator = (
81-
feature for feature in gb.features
82-
if any(
83-
(qualifier_value in included.get(qualifier_key, []))
84-
for qualifier_key in included.keys()
85-
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
86-
) and all(
87-
(qualifier_value not in excluded.get(qualifier_key, []))
88-
for qualifier_key in excluded.keys()
89-
for qualifier_value in feature.qualifiers.get(qualifier_key, [])
90-
)
91-
)
92-
93119
logging.info("Indexing feature locations")
94-
feature_index = index_features(feature_iterator)
120+
included = snakemake.params.features.get("INCLUDE", {})
121+
excluded = snakemake.params.features.get("EXCLUDE", {})
122+
feature_index = index_features(
123+
iter_features_filtering(gb.features, included, excluded)
124+
)
95125

96126
logging.info("Calculating polimorphic sites sliding window")
97127
windows = window_calculation(

0 commit comments

Comments
 (0)