Coverage for rest_api/data_entry/annotation_import.py: 20%
102 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-10 03:22 +0000
1from dataclasses import dataclass
2import gzip
3import re
5from django.db.models import Q
7from rest_api.models import AnnotationType
8from rest_api.models import NucleotideMutation
9from sonar_backend.settings import LOGGER
12@dataclass
13class VCFRaw:
14 chrom: str # mutation__replicon__accession
15 pos: int # mutation__start
16 id: str
17 ref: str # mutation__ref
18 alt: str | None
19 qual: str
20 filter: str
21 info: str
22 format: str
25@dataclass
26class VCFInfoANNRaw:
27 # Functional annotations
28 allele: str
29 annotation: str
30 annotation_impact: str
31 gene_name: str
32 gene_id: str
33 feature_type: str
34 feature_id: str
35 transcript_bio_type: str
36 rank: str
37 hgvs_c: str
38 hgvs_p: str
39 c_dna_pos_length: str
40 cds_pos_length: str
41 aa_pos_length: str
42 distance: str
43 errors_warnings_info: str
46@dataclass
47class MutationLookupToAnnotations:
48 start: str
49 end: str
50 ref: str
51 alt: str
52 replicon__accession: str
53 annotations: list[VCFInfoANNRaw]
56class AnnotationImport:
57 def __init__(self, path: str):
58 self.vcf_file_path = path
59 self.raw_lines = [line for line in self._import_vcf()]
60 self.mutation_lookups_to_annotations = self.convert_lines()
61 self.annotation_q_obj = Q()
63 def _import_vcf(self):
64 # Open compressed or uncompressed VCF based on file extension
65 open_func = gzip.open if self.vcf_file_path.endswith(".gz") else open
66 with open_func(self.vcf_file_path, "rt") as handle: # 'rt' for text mode
67 for line in handle:
68 if line.startswith(("#", "##")):
69 continue
70 line = line.strip("\r\n").split("\t")
71 vcf_raw = VCFRaw(
72 chrom=line.pop(0),
73 pos=int(line.pop(0)),
74 id=line.pop(0),
75 ref=line.pop(0),
76 alt=None if (alt := line.pop(0)) == "." else alt,
77 qual=line.pop(0),
78 filter=line.pop(0),
79 info=line.pop(0),
80 format=line.pop(0),
81 )
82 yield vcf_raw
84 def convert_lines(self) -> list[MutationLookupToAnnotations]:
85 mutation_lookups_to_annotations = []
86 for line in self.raw_lines:
87 annotations = self._parse_line_info(line.info)
88 allele_to_annotations = {}
89 for annotation in annotations:
90 if annotation.allele not in allele_to_annotations:
91 allele_to_annotations[annotation.allele] = []
92 allele_to_annotations[annotation.allele].append(annotation)
93 for alt in line.alt.split(","):
94 if alt not in allele_to_annotations:
95 continue
96 mutation_lookup_to_annotations = MutationLookupToAnnotations(
97 start=line.pos
98 - 1, # snp (we deduct by one because our database use 0-based)
99 end=int(line.pos),
100 ref=line.ref,
101 alt=alt,
102 replicon__accession=line.chrom,
103 annotations=allele_to_annotations[alt],
104 )
105 if len(alt) < len(mutation_lookup_to_annotations.ref):
106 # deletion and alt not null
107 # because in vcf it always show the positon before deletion
108 # for example; MN908947.3 506 . CATGGTCATGTTATGGTTG C
109 mutation_lookup_to_annotations.start += 1
110 mutation_lookup_to_annotations.end = (
111 mutation_lookup_to_annotations.end
112 + len(mutation_lookup_to_annotations.ref)
113 - 1
114 )
115 # add None here if we dont want to keep the deletion in ref
116 # column, however the program frozen once I change to
117 # mutation_lookup_to_annotations.ref = None
118 mutation_lookup_to_annotations.ref = ""
119 mutation_lookup_to_annotations.alt = ""
121 mutation_lookups_to_annotations.append(mutation_lookup_to_annotations)
122 return mutation_lookups_to_annotations
124 def _parse_line_info(self, info: str) -> list[VCFInfoANNRaw]:
125 """
126 Extracts the SNP effect annotation (ANN) substring from the 8th column of a tab-separated SnpEff annotation line.
128 Parameters:
129 info (str): The 8th column of a tab-separated SnpEff annotation line. Example:
130 "ANN=C|upstream_gene_variant|MODIFIER|ORF1a|Gene_265_13467|transcript|ORF1a|protein_coding||c.-217_-213delTCTTG|||||217|WARNING_TRANSCRIPT_NO_STOP_CODON,
131 C|intergenic_region|MODIFIER|CHR_START-ORF1a|CHR_START-Gene_265_13467|intergenic_region|CHR_START-Gene_265_13467|||n.49_53delTCTTG||||||"
133 Returns:
134 list[VCFInfoANNRaw]: list of extracted annotations, different annotations are separated by ','
135 """
136 for field in info.split(";"):
137 if field.startswith("ANN="):
138 ann_field = field.removeprefix("ANN=")
139 annotations = []
140 for annotation in ann_field.split(","):
141 annotation = annotation.split("|")
142 try:
143 annotations.append(VCFInfoANNRaw(*annotation))
144 except Exception:
145 LOGGER.warning(
146 f"Failed to parse annotation: {annotation}, from file {self.vcf_file_path}"
147 )
148 # Return after finding the first ANN= field. If there are
149 # multiple, all others will be ignored.
150 return annotations
151 return []
153 def get_annotation_objs(
154 self,
155 ) -> tuple[list[AnnotationType], list[AnnotationType.mutations.through]]:
156 annotation_objs = []
157 q_obj = Q()
158 # self.mutation_lookups_to_annotations read from vcf file.
159 for mutation_lookup_to_annotations in self.mutation_lookups_to_annotations:
160 q_obj |= Q(
161 start=mutation_lookup_to_annotations.start,
162 end=mutation_lookup_to_annotations.end,
163 ref=mutation_lookup_to_annotations.ref,
164 alt=mutation_lookup_to_annotations.alt,
165 replicon__accession=mutation_lookup_to_annotations.replicon__accession,
166 )
167 mutations = NucleotideMutation.objects.filter(q_obj).prefetch_related(
168 "replicon"
169 )
170 annotation_q_obj = Q()
171 relation_info = {}
172 for mutation in mutations:
173 # Problem 1: some samples got error 'StopIteration'
174 # I think because there are no items in the filtered iterable
175 # that match the conditions specified by the lambda function
176 # But How did this can happen?,
177 # Solution: because there are cds and nt at the same position we should filter only NT
179 try:
180 mut_lookup_to_annotation = self.mutation_lookups_to_annotations.pop(
181 self.mutation_lookups_to_annotations.index(
182 next(
183 filter(
184 lambda x: int(x.start) == int(mutation.start)
185 and int(x.end) == int(mutation.end)
186 and x.ref == mutation.ref
187 and x.alt == mutation.alt
188 and x.replicon__accession
189 == mutation.replicon.accession,
190 self.mutation_lookups_to_annotations,
191 )
192 )
193 )
194 )
195 except (ValueError, StopIteration) as e:
196 LOGGER.error(f"Error: {e}")
197 LOGGER.error(
198 f"Mutation details: start={mutation.start}, end={mutation.end}, ref={mutation.ref}, alt={mutation.alt}, replicon__accession={mutation.replicon.accession}"
199 )
200 raise
202 # print(mut_lookup_to_annotation)
203 # Problem 2: We have too many entries in mut_lookup_to_annotation.annotations.
204 # For example, if we have the mutation MN908947.3 26565 . A ANNNNNNN (insertion with ambiguous lots of Ns),
205 # snpEff tries to predict the effect for every possible combination:
206 # insCAAAAAA, insCAAAAAC, insCAAAAAG, insCAAAAAT, ..., insTAAAAAA, ..., insGAAAAAA.
207 # This can generate a lookup annotation size of 4 (A,T,C,G) to the power of N (position).
208 # In this case, 4 to the power of 7 equals 16384.
209 # Currently, we only use ontology (e.g., frameshift_variant) and annotation_impact (e.g., HIGH),
210 # which means a potentially highly redundant query, hence taking too long to finish the import process (> 10 mins).
212 # Temporary solution: reduce the redundancy in alleles, annotations, and impacts.
213 # we skip processing based on alleles, annotations, and impacts
214 seen = set()
215 # start to pass each VCFInfoANNRaw
216 for a in mut_lookup_to_annotation.annotations:
217 key = (a.allele, a.annotation, a.annotation_impact)
218 if key in seen:
219 continue
220 seen.add(key)
221 for ontology in a.annotation.split("&"):
222 annotation_obj = AnnotationType(
223 seq_ontology=ontology, impact=a.annotation_impact
224 )
225 annotation_q_obj |= Q(
226 seq_ontology=ontology, impact=a.annotation_impact
227 )
228 if ontology not in relation_info:
229 relation_info[ontology] = {}
230 if a.annotation_impact not in relation_info[ontology]:
231 relation_info[ontology][a.annotation_impact] = []
232 relation_info[ontology][a.annotation_impact].append(
233 {
234 "mutation": mutation,
235 }
236 )
237 annotation_objs.append(annotation_obj)
238 self.annotation_q_obj = annotation_q_obj
239 self.relation_info = relation_info
240 return annotation_objs
242 def get_annotation2mutation_objs(self) -> list[AnnotationType.mutations.through]:
243 annotations = AnnotationType.objects.filter(self.annotation_q_obj)
244 mutation2annotation_objs = []
245 for annotation in annotations:
246 for relation in self.relation_info[annotation.seq_ontology][
247 annotation.impact
248 ]:
249 mutation = relation["mutation"]
250 mutation2annotation_objs.append(
251 AnnotationType.mutations.through(
252 annotationtype_id=annotation.id,
253 nucleotidemutation_id=mutation.id,
254 )
255 )
256 return mutation2annotation_objs