25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220 | class TypingResult:
"""
This is a class to store the results of a typing analysis for a single sample. It is designed to be flexible
enough to store results from both read and assembly typing, and can be easily reconstructed from JSON to use
with the `convert` utility. It should not store any information that is not directly related to the typing
such as contig sequences or read alignments.
"""
def __init__(
self, sample_name: str | None, db: Database | None, best_match: Locus = None,
pieces: list[LocusPiece] = None, expected_genes_inside_locus: list[GeneResult] = None,
expected_genes_outside_locus: list[GeneResult] = None, missing_genes: list[str] = None,
unexpected_genes_inside_locus: list[GeneResult] = None,
unexpected_genes_outside_locus: list[GeneResult] = None,
extra_genes: list[GeneResult] = None):
self.sample_name = sample_name or ""
self.db = db
self.best_match = best_match
self.pieces = pieces or [] # Pieces of locus reconstructed from alignments
self.expected_genes_inside_locus = expected_genes_inside_locus or [] # Genes from best_match
self.expected_genes_outside_locus = expected_genes_outside_locus or [] # Genes from best_match
self.missing_genes = missing_genes or [] # Genes from best_match that were not found
self.unexpected_genes_inside_locus = unexpected_genes_inside_locus or [] # Genes from other loci
self.unexpected_genes_outside_locus = unexpected_genes_outside_locus or [] # Genes from other loci
self.extra_genes = extra_genes or [] # in db.extra_genes, ALWAYS outside locus (gene_result.piece == None)
# Properties to cache the values, these are protected to prevent accidental modification
self._percent_identity = None
self._percent_coverage = None
self._phenotype = None
self._problems = None
self._confidence = None
def __repr__(self):
return f"{self.sample_name} {self.best_match.name}"
def __len__(self):
return sum(len(i) for i in self.pieces) if self.pieces else 0
def __iter__(self):
return chain(
self.expected_genes_inside_locus, self.unexpected_genes_inside_locus,
self.expected_genes_outside_locus, self.unexpected_genes_outside_locus, self.extra_genes)
def add_gene_result(self, gene_result: GeneResult):
if gene_result.piece: # If gene_result.piece is not None, the gene is inside the locus
gene_result.piece.add_gene_result(gene_result)
gene_type = f"{gene_result.gene_type}{'_inside_locus' if gene_result.gene_type.startswith(('expected', 'unexpected')) else ''}"
else: # If gene_result.piece is None, the gene is outside the locus
gene_type = f"{gene_result.gene_type}{'_outside_locus' if gene_result.gene_type.startswith(('expected', 'unexpected')) else ''}"
getattr(self, gene_type).append(gene_result) # Add gene result to the appropriate list
@property
def percent_identity(self) -> float:
if self._percent_identity is None:
self._percent_identity = (sum(i.percent_identity for i in x) / len(x)) if (
x := self.expected_genes_inside_locus) else 0
return self._percent_identity
@property
def percent_coverage(self) -> float:
if self._percent_coverage is None:
self._percent_coverage = sum(len(i) for i in x) / sum(len(i) for i in self.best_match.genes.values()) * 100 \
if (x := self.expected_genes_inside_locus) else 0
return self._percent_coverage
@property
def phenotype(self) -> str:
if self._phenotype is None:
gene_phenotypes = set() # Init set to store gene phenotypes to be used as a key in the phenotypes dict
for gene in self:
if gene.gene_type in {'expected_genes', 'extra_genes'}: # The reported phenotype only considers these
gene_phenotypes.add((gene.gene.name, gene.phenotype)) # or extra genes
# NOTE: The best_match.phenotypes MUST be sorted from largest to smallest gene set to make sure any sets with
# extra genes are tested first.
self._phenotype = next(
(p for gs, p in self.best_match.phenotypes if len(gs) == len(gene_phenotypes.intersection(gs))),
self.best_match.type_label) # If no phenotype is found, return the type label
return self._phenotype
@property
def problems(self) -> str:
if self._problems is None:
self._problems = f'?{x}' if (x := len(self.pieces)) != 1 else ''
self._problems += '-' if self.missing_genes else ''
self._problems += '+' if self.unexpected_genes_inside_locus else ''
self._problems += '*' if any(
i.percent_coverage >= 90 and i.below_threshold for i in self.expected_genes_inside_locus) else ''
self._problems += '!' if any(i.phenotype == "truncated" for i in self) else ''
return self._problems
@property
def confidence(self) -> str:
return self._confidence if self._confidence is not None else "Not calculated"
def get_confidence(self, allow_below_threshold: bool, max_other_genes: int, percent_expected_genes: float):
p = len(set(i.gene.name for i in self.expected_genes_inside_locus)) / len(self.best_match.genes) * 100
other_genes = len(
set(i.gene.name for i in self.unexpected_genes_inside_locus if not i.phenotype == "truncated"))
if not allow_below_threshold and "*" in self.problems:
self._confidence = "Untypeable"
else:
if len(self.pieces) == 1 and not self.missing_genes and not other_genes:
self._confidence = "Typeable"
elif other_genes <= max_other_genes and p >= percent_expected_genes:
self._confidence = "Typeable"
else:
self._confidence = "Untypeable"
@classmethod
def from_dict(cls, d: dict, db: Database) -> TypingResult:
if not (best_match := db.loci.get(d['best_match'])):
raise TypingResultError(f"Best match {d['best_match']} not found in database")
self = TypingResult(sample_name=d['sample_name'], db=db, best_match=best_match,
missing_genes=d['missing_genes'])
# Set the cached properties
self._percent_identity = float(d['percent_identity'])
self._percent_coverage = float(d['percent_coverage'])
self._phenotype = d['phenotype']
self._problems = d['problems']
self._confidence = d['confidence']
# Add the pieces and create the gene results
self.pieces = [LocusPiece.from_dict(i, result=self) for i in d['pieces']]
pieces = {i.__repr__(): i for i in self.pieces}
gene_results = {} # This was previously a dict comp, but we need to check the gene is in the database, see #31
for r in chain(d['expected_genes_inside_locus'], d['unexpected_genes_inside_locus'],
d['expected_genes_outside_locus'], d['unexpected_genes_outside_locus'],
d['extra_genes']):
if not (gene := db.genes.get(r['gene'])) and not (gene := db.extra_genes.get(r['gene'])):
raise TypingResultError(f"Gene {r['gene']} not found in database")
x = GeneResult.from_dict(r, result=self, piece=pieces.get(r['piece']), gene=gene)
gene_results[x.__repr__()] = x
for gene_result in gene_results.values():
self.add_gene_result(gene_result)
return self
def format(self, format_spec) -> str | dict:
if format_spec == 'tsv':
return '\t'.join(
[
self.sample_name,
self.best_match.name,
self.phenotype,
self.confidence,
self.problems,
f"{self.percent_identity:.2f}%",
f"{self.percent_coverage:.2f}%",
f"{self.__len__() - len(self.best_match)} bp" if len(self.pieces) == 1 else 'n/a',
f"{(n_inside := len({i.gene.name for i in self.expected_genes_inside_locus}))} / {(n_expected := len(self.best_match.genes))} ({100 * n_inside / n_expected:.2f}%)",
';'.join(map(str, self.expected_genes_inside_locus)),
';'.join(self.missing_genes),
f"{len(self.unexpected_genes_inside_locus)}",
';'.join(map(str, self.unexpected_genes_inside_locus)),
f"{(n_outside := len({i.gene.name for i in self.expected_genes_outside_locus}))} / {n_expected} ({100 * n_outside / n_expected:.2f}%)",
';'.join(map(str, self.expected_genes_outside_locus)),
f"{len(self.unexpected_genes_outside_locus)}",
';'.join(map(str, self.unexpected_genes_outside_locus)),
';'.join(map(str, filter(lambda z: z.phenotype == "truncated", self))),
';'.join(map(str, self.extra_genes))
]
) + "\n"
if format_spec == 'fna': # Return the nucleotide sequence of the locus
return "".join([i.format(format_spec) for i in self.pieces])
if format_spec in {'faa', 'ffn'}: # Return the protein or nucleotide sequence of gene results
return "".join([i.format(format_spec) for i in self])
if format_spec == 'json':
return dumps(
{
'sample_name': self.sample_name, 'best_match': self.best_match.name, 'confidence': self.confidence,
'phenotype': self.phenotype, 'problems': self.problems,
'percent_identity': str(self.percent_identity),
'percent_coverage': str(self.percent_coverage), 'missing_genes': self.missing_genes
} | {
attr: [i.format(format_spec) for i in getattr(self, attr)] for attr in {
'pieces', 'expected_genes_inside_locus', 'unexpected_genes_inside_locus',
'expected_genes_outside_locus', 'unexpected_genes_outside_locus', 'extra_genes'
}
}) + "\n"
raise ValueError(f"Unknown format specifier {format_spec}")
def write(self,
tsv: TextIO = None,
json: TextIO = None,
fna: str | PathLike | TextIO = None,
ffn: str | PathLike | TextIO = None,
faa: str | PathLike | TextIO = None,
):
"""Write the typing result to files or file handles."""
[f.write(self.format(fmt)) for f, fmt in [(tsv, 'tsv'), (json, 'json')] if isinstance(f, TextIOBase)]
for f, fmt in [(fna, 'fna'), (ffn, 'ffn'), (faa, 'faa')]:
if f:
if isinstance(f, TextIOBase):
f.write(self.format(fmt))
elif isinstance(f, PathLike) or isinstance(f, str):
with open(path.join(f, f'{self.sample_name}_kaptive_results.{fmt}'), 'wt') as handle:
handle.write(self.format(fmt))
|