1919GWAS_CONFIG_DIR = REPO_ROOT / "config" / "gwas"
2020
2121
22- def _count_effect_alleles (genotype : str | None , effect_allele : str | None ) -> int | None :
22+ _COMPLEMENT = str .maketrans ("ACGT" , "TGCA" )
23+ _PALINDROMIC = {frozenset (("A" , "T" )), frozenset (("C" , "G" ))}
24+
25+
26+ def _count_effect_alleles (
27+ genotype : str | None ,
28+ effect_allele : str | None ,
29+ other_allele : str | None = None ,
30+ ) -> int | None :
2331 """Count copies of the effect allele in a diploid genotype like 'AG' or 'A/G'.
2432
25- Returns None if we can't determine (missing data, indel, etc).
33+ When ``other_allele`` is supplied, the genotype is validated against the
34+ GWAS {effect, other} allele pair so we can orient it to the reference strand:
35+
36+ * If the genotype alleles already fall within {effect, other}, count directly.
37+ * Otherwise try the complement strand; count only if it then matches.
38+ * Palindromic SNPs (A/T, C/G) are their own complement and therefore
39+ un-orientable from genotype alone — return None.
40+ * Genotypes matching neither orientation are unresolvable — return None.
41+
42+ Returns None if we can't determine (missing data, indel, strand-ambiguous, etc).
2643 """
2744 if not genotype or not effect_allele :
2845 return None
@@ -34,9 +51,52 @@ def _count_effect_alleles(genotype: str | None, effect_allele: str | None) -> in
3451 return None
3552 if len (ea ) != 1 :
3653 return None
54+
55+ oa = other_allele .upper () if other_allele else None
56+ if oa and len (oa ) != 1 :
57+ return None
58+
59+ if oa :
60+ expected = {ea , oa }
61+ # Palindromic SNPs are un-orientable from genotype alone.
62+ if frozenset (expected ) in _PALINDROMIC :
63+ return None
64+
65+ if set (g ) <= expected :
66+ return sum (1 for base in g if base == ea )
67+
68+ comp = g .translate (_COMPLEMENT )
69+ if set (comp ) <= expected :
70+ return sum (1 for base in comp if base == ea )
71+
72+ # Matches neither orientation — unresolvable.
73+ return None
74+
3775 return sum (1 for base in g if base == ea )
3876
3977
78+ def _count_traits_per_entry (
79+ user_genotype : str | None , traits_list : list [dict ]
80+ ) -> list [dict ]:
81+ """Attach a per-trait ``effect_allele_count`` to each pleiotropic trait entry.
82+
83+ Each trait may carry a different effect/other allele for the same SNP, so the
84+ dosage is computed independently per entry rather than reusing one trait's
85+ effect allele for all (audit finding #12).
86+ """
87+ return [
88+ {
89+ ** t ,
90+ "effect_allele_count" : _count_effect_alleles (
91+ user_genotype ,
92+ t .get ("effect_allele" ),
93+ t .get ("other_allele" ),
94+ ),
95+ }
96+ for t in traits_list
97+ ]
98+
99+
40100@router .get ("/traits" )
41101async def list_traits ():
42102 """List traits that have pre-computed GWAS hits available."""
@@ -108,12 +168,11 @@ async def get_gwas_overlap():
108168
109169 results : list [dict ] = []
110170 for rsid , traits_list in pleiotropic .items ():
111- snp = await genome_db .get_snp (rsid )
112- # Use the first trait's effect allele for counting (they may differ across studies)
113- ea = traits_list [0 ].get ("effect_allele" )
114- ea_count = _count_effect_alleles (
115- snp .get ("genotype" ) if snp else None , ea
116- )
171+ snp = await genome_db .get_snp (rsid , profile_id = "default" )
172+ user_genotype = snp .get ("genotype" ) if snp else None
173+ # Count the effect-allele dosage PER trait — a pleiotropic SNP can carry
174+ # different effect/other alleles across studies (audit finding #12).
175+ traits_with_counts = _count_traits_per_entry (user_genotype , traits_list )
117176
118177 avg_p = sum (t ["p_value" ] for t in traits_list if t ["p_value" ]) / max (
119178 sum (1 for t in traits_list if t ["p_value" ]), 1
@@ -125,9 +184,8 @@ async def get_gwas_overlap():
125184 "pos" : rsid_info [rsid ]["pos" ],
126185 "n_traits" : len (traits_list ),
127186 "avg_p_value" : avg_p ,
128- "user_genotype" : snp .get ("genotype" ) if snp else None ,
129- "effect_allele_count" : ea_count ,
130- "traits" : traits_list ,
187+ "user_genotype" : user_genotype ,
188+ "traits" : traits_with_counts ,
131189 })
132190
133191 # Sort: most traits first, then lowest average p-value
@@ -175,10 +233,10 @@ async def get_gwas_summary():
175233 rsid = hit .get ("rsid" )
176234 if not rsid :
177235 continue
178- snp = await genome_db .get_snp (rsid )
236+ snp = await genome_db .get_snp (rsid , profile_id = "default" )
179237 if not snp :
180238 continue
181- ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ))
239+ ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ), hit . get ( "other_allele" ) )
182240 if ea_count is None :
183241 continue
184242
@@ -250,10 +308,10 @@ async def _compute_trait_summary(trait: str) -> dict | None:
250308 rsid = hit .get ("rsid" )
251309 if not rsid :
252310 continue
253- snp = await genome_db .get_snp (rsid )
311+ snp = await genome_db .get_snp (rsid , profile_id = "default" )
254312 if not snp :
255313 continue
256- ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ))
314+ ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ), hit . get ( "other_allele" ) )
257315 if ea_count is None :
258316 continue
259317
@@ -358,11 +416,11 @@ async def get_addiction_gwas_summary():
358416 if not rsid :
359417 continue
360418
361- snp = await genome_db .get_snp (rsid )
419+ snp = await genome_db .get_snp (rsid , profile_id = "default" )
362420 if not snp :
363421 continue
364422
365- ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ))
423+ ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ), hit . get ( "other_allele" ) )
366424 if ea_count is None :
367425 continue
368426
@@ -456,11 +514,11 @@ async def get_gwas_matches(trait: str, clumped: bool = Query(False)):
456514 if not rsid :
457515 continue
458516
459- snp = await genome_db .get_snp (rsid )
517+ snp = await genome_db .get_snp (rsid , profile_id = "default" )
460518 if not snp :
461519 continue
462520
463- ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ))
521+ ea_count = _count_effect_alleles (snp .get ("genotype" ), hit .get ("effect_allele" ), hit . get ( "other_allele" ) )
464522 if ea_count is None :
465523 continue
466524
0 commit comments