Skip to content

Commit 339e0b1

Browse files
madhavajaymadhavajay
authored
Merge pull request #73 from OpenMined/madhava/fix-aliases
fixing rsid alias matching
2 parents b61fde6 + 8809137 commit 339e0b1

8 files changed

Lines changed: 170 additions & 22 deletions

File tree

‎rust/bioscript-formats/src/genotype.rs‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -712,6 +712,15 @@ mod tests {
712712
Some("--")
713713
);
714714
assert_eq!(genotype_from_vcf_gt("2/2", "A", &["G"]), None);
715+
assert_eq!(
716+
genotype_from_vcf_gt("0/1", "TTTTTTTTTTTTT", &["TTTTTTTTTTTT"]).as_deref(),
717+
Some("TTTTTTTTTTTTT/TTTTTTTTTTTT")
718+
);
719+
assert_eq!(
720+
genotype_from_vcf_gt("0/2", "TTTTTTTTTTTTT", &["TTTTTTTTTTTT", "TTTTTTTTTTTTTT"])
721+
.as_deref(),
722+
Some("TTTTTTTTTTTTT/TTTTTTTTTTTTTT")
723+
);
715724
assert_eq!(vcf_reference_token("AT", &["A"]), "I");
716725
assert_eq!(vcf_reference_token("A", &["AT"]), "D");
717726
assert_eq!(vcf_reference_token("A", &["<NON_REF>"]), "A");

‎rust/bioscript-formats/src/genotype/vcf.rs‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -309,7 +309,10 @@ fn resolve_vcf_row(
309309
backend: backend.backend_name().to_owned(),
310310
matched_rsid: Some(rsid.clone()),
311311
assembly: targets.detected_assembly,
312-
genotype: Some(row.genotype.clone()),
312+
genotype: Some(matching::vcf_row_genotype_for_variant(
313+
row,
314+
&targets.variants[target_idx],
315+
)),
313316
evidence: vec![
314317
format!("resolved by rsid {rsid}"),
315318
format!("source line: {}", row.raw_line),
@@ -339,7 +342,10 @@ fn resolve_vcf_row(
339342
backend: backend.backend_name().to_owned(),
340343
matched_rsid: row.rsid.clone(),
341344
assembly: targets.detected_assembly,
342-
genotype: Some(row.genotype.clone()),
345+
genotype: Some(matching::vcf_row_genotype_for_variant(
346+
row,
347+
&targets.variants[target_idx],
348+
)),
343349
evidence: vec![
344350
format!("resolved by locus {}:{}", row.chrom, row.position),
345351
format!("source line: {}", row.raw_line),

‎rust/bioscript-formats/src/genotype/vcf/matching.rs‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,51 @@ pub(crate) fn vcf_row_matches_variant(
142142
}
143143
}
144144

145+
pub(crate) fn vcf_row_genotype_for_variant(row: &ParsedVcfRow, variant: &VariantSpec) -> String {
146+
if !matches!(
147+
variant.kind,
148+
Some(VariantKind::Deletion | VariantKind::Insertion | VariantKind::Indel)
149+
) {
150+
return row.genotype.clone();
151+
}
152+
153+
let parts: Vec<&str> = row.genotype.split('/').collect();
154+
if parts.len() <= 1 {
155+
return row.genotype.clone();
156+
}
157+
158+
let alternates: Vec<&str> = row.alternates.iter().map(String::as_str).collect();
159+
let mut tokens = Vec::with_capacity(parts.len());
160+
for part in parts {
161+
if part.eq_ignore_ascii_case(&row.reference) {
162+
tokens.push(super::super::vcf_tokens::vcf_reference_token(
163+
&row.reference,
164+
&alternates,
165+
));
166+
} else if let Some(alternate) = row
167+
.alternates
168+
.iter()
169+
.find(|alternate| alternate.eq_ignore_ascii_case(part))
170+
{
171+
tokens.push(super::super::vcf_tokens::vcf_alt_token(
172+
&row.reference,
173+
alternate,
174+
));
175+
} else {
176+
return row.genotype.clone();
177+
}
178+
}
179+
180+
if tokens
181+
.iter()
182+
.all(|token| token.chars().count() == 1 && token != "--")
183+
{
184+
return super::super::normalize_genotype(&tokens.join(""));
185+
}
186+
187+
row.genotype.clone()
188+
}
189+
145190
fn snp_row_has_catalog_allele(row: &ParsedVcfRow, variant: &VariantSpec) -> bool {
146191
let Some(alternate) = variant.alternate.as_ref() else {
147192
return true;

‎rust/bioscript-formats/src/genotype/vcf/reader.rs‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ use noodles::tabix;
77

88
use bioscript_core::{Assembly, GenomicLocus, RuntimeError, VariantObservation, VariantSpec};
99

10-
use super::{parse_vcf_record, vcf_row_matches_variant};
10+
use super::{matching::vcf_row_genotype_for_variant, parse_vcf_record, vcf_row_matches_variant};
1111

1212
/// Observe a SNP at `locus` over an already-built tabix-indexed bgzipped VCF
1313
/// reader. Caller builds `csi::io::IndexedReader::new(reader, tabix_index)`
@@ -208,7 +208,7 @@ where
208208
backend: "vcf".to_owned(),
209209
matched_rsid: matched_rsid.or_else(|| row.rsid.clone()),
210210
assembly,
211-
genotype: Some(row.genotype.clone()),
211+
genotype: Some(vcf_row_genotype_for_variant(&row, variant)),
212212
evidence: vec![format!("{label}: resolved by indexed locus {locus_label}")],
213213
..VariantObservation::default()
214214
});

‎rust/bioscript-formats/src/genotype/vcf_tokens.rs‎

Lines changed: 24 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -15,21 +15,40 @@ pub(crate) fn genotype_from_vcf_gt(
1515
return Some("--".to_owned());
1616
}
1717

18-
let ref_token = vcf_reference_token(reference, alternates);
19-
let mut out = String::new();
18+
let mut indexes = Vec::with_capacity(parts.len());
2019
for part in parts {
2120
let Ok(idx) = part.parse::<usize>() else {
2221
return Some("--".to_owned());
2322
};
23+
indexes.push(idx);
24+
}
25+
26+
let alternate_only = indexes.iter().all(|idx| *idx > 0);
27+
let mut tokens = Vec::with_capacity(indexes.len());
28+
for idx in indexes {
2429
if idx == 0 {
25-
out.push_str(&ref_token);
30+
tokens.push(normalize_sequence_token(reference));
2631
} else {
2732
let alt = alternates.get(idx - 1)?;
28-
out.push_str(&vcf_alt_token(reference, alt));
33+
if is_symbolic_vcf_alt(alt) {
34+
return Some("--".to_owned());
35+
}
36+
if alternate_only {
37+
tokens.push(vcf_alt_token(reference, alt));
38+
} else {
39+
tokens.push(normalize_sequence_token(alt));
40+
}
2941
}
3042
}
3143

32-
Some(normalize_genotype(&out))
44+
if tokens
45+
.iter()
46+
.all(|token| token.chars().count() == 1 && token != "--")
47+
{
48+
return Some(normalize_genotype(&tokens.join("")));
49+
}
50+
51+
Some(tokens.join("/"))
3352
}
3453

3554
pub(crate) fn vcf_reference_token(reference: &str, alternates: &[&str]) -> String {

‎rust/bioscript-reporting/src/manifest_catalogue.rs‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -87,8 +87,13 @@ fn catalogue_row_task(
8787
columns.value(row, "identifier.aliases"),
8888
columns.separator("identifier.aliases"),
8989
));
90-
rsids.sort();
91-
rsids.dedup();
90+
let mut deduped_rsids = Vec::with_capacity(rsids.len());
91+
for rsid in rsids {
92+
if !deduped_rsids.iter().any(|seen| seen == &rsid) {
93+
deduped_rsids.push(rsid);
94+
}
95+
}
96+
let rsids = deduped_rsids;
9297
let alternates = split_list(
9398
columns.value(row, "alleles.alts"),
9499
columns.separator("alleles.alts"),

‎rust/bioscript-reporting/src/observation.rs‎

Lines changed: 38 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -338,6 +338,12 @@ fn render_app_observation_json(input: AppObservationJson) -> serde_json::Value {
338338
} else {
339339
gene
340340
};
341+
let canonical_rsid = manifest.spec.rsids.first().filter(|rsid| !rsid.is_empty());
342+
let matched_rsid = row.get("matched_rsid").filter(|rsid| !rsid.is_empty());
343+
let alias_match_note = matched_rsid
344+
.filter(|matched| canonical_rsid.is_some_and(|canonical| canonical != *matched))
345+
.filter(|matched| manifest.spec.rsids.iter().any(|rsid| rsid == *matched))
346+
.map(|matched| format!("matched alias: {matched}"));
341347
let source = if source.is_null() {
342348
manifest_default_source(&row, &manifest)
343349
} else {
@@ -349,7 +355,7 @@ fn render_app_observation_json(input: AppObservationJson) -> serde_json::Value {
349355
"assay_version": "1.0",
350356
"variant_key": manifest.name,
351357
"variant_path": row_path,
352-
"rsid": row.get("matched_rsid").filter(|value| !value.is_empty()).cloned().or_else(|| manifest.spec.rsids.first().cloned()),
358+
"rsid": canonical_rsid.cloned().or_else(|| matched_rsid.cloned()),
353359
"gene": gene,
354360
"assembly": if assembly.is_empty() { serde_json::Value::Null } else { serde_json::Value::String(assembly.to_uppercase()) },
355361
"chrom": chrom,
@@ -377,6 +383,8 @@ fn render_app_observation_json(input: AppObservationJson) -> serde_json::Value {
377383
"match_quality": if weak_indel_match { serde_json::Value::String("weak".to_owned()) } else { serde_json::Value::Null },
378384
"match_notes": if weak_indel_match {
379385
serde_json::Value::String("consumer genotype file reported an insertion/deletion token at the marker, not sequence-resolved evidence for the exact deletion allele".to_owned())
386+
} else if let Some(note) = alias_match_note {
387+
serde_json::Value::String(note)
380388
} else {
381389
serde_json::Value::Null
382390
},
@@ -397,10 +405,10 @@ fn manifest_default_source(
397405
manifest: &VariantManifest,
398406
) -> serde_json::Value {
399407
let rsid = row
400-
.get("matched_rsid")
408+
.get("rsid")
401409
.filter(|rsid| !rsid.is_empty())
402-
.or_else(|| row.get("rsid").filter(|rsid| !rsid.is_empty()))
403-
.or_else(|| manifest.spec.rsids.first().filter(|rsid| !rsid.is_empty()));
410+
.or_else(|| manifest.spec.rsids.first().filter(|rsid| !rsid.is_empty()))
411+
.or_else(|| row.get("matched_rsid").filter(|rsid| !rsid.is_empty()));
404412
let Some(rsid) = rsid else {
405413
return serde_json::Value::Null;
406414
};
@@ -445,6 +453,32 @@ mod tests {
445453
);
446454
}
447455

456+
#[test]
457+
fn repeat_indel_insertion_deletion_tokens_remain_ambiguous_without_sequence_alleles() {
458+
assert_eq!(
459+
normalize_app_genotype(
460+
"II",
461+
"TTTTTTTTTTTTT",
462+
"TTTTTTTTTTTT",
463+
Some(VariantKind::Indel),
464+
"19",
465+
None,
466+
),
467+
("II".to_owned(), "unknown".to_owned())
468+
);
469+
assert_eq!(
470+
normalize_app_genotype(
471+
"ID",
472+
"TTTTTTTTTTTTT",
473+
"TTTTTTTTTTTT",
474+
Some(VariantKind::Indel),
475+
"19",
476+
None,
477+
),
478+
("ID".to_owned(), "unknown".to_owned())
479+
);
480+
}
481+
448482
#[test]
449483
fn displays_cram_long_deletion_copy_number_as_insertion_deletion_tokens() {
450484
let manifest = VariantManifest {

‎rust/bioscript-reporting/src/observation/genotype_display.rs‎

Lines changed: 37 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -58,14 +58,17 @@ pub(super) fn normalize_app_genotype(
5858
if display.is_empty() {
5959
return ("./.".to_owned(), "unknown".to_owned());
6060
}
61-
if matches!(kind, Some(VariantKind::Deletion))
62-
&& ref_allele.len() != 1
63-
&& display
64-
.chars()
65-
.filter(char::is_ascii_alphabetic)
66-
.all(|allele| matches!(allele.to_ascii_uppercase(), 'I' | 'D'))
61+
if let Some((reference_token, alternate_token)) =
62+
indel_display_tokens(display, ref_allele, alt_allele, kind)
6763
{
68-
return normalize_app_genotype(display, "I", "D", None, chrom, inferred_sex);
64+
return normalize_app_genotype(
65+
display,
66+
&reference_token,
67+
&alternate_token,
68+
None,
69+
chrom,
70+
inferred_sex,
71+
);
6972
}
7073
if let Some(normalized) = normalize_long_allele_genotype(display, ref_allele, alt_allele) {
7174
return normalized;
@@ -109,6 +112,33 @@ pub(super) fn normalize_app_genotype(
109112
}
110113
}
111114

115+
fn indel_display_tokens(
116+
display: &str,
117+
ref_allele: &str,
118+
_alt_allele: &str,
119+
kind: Option<VariantKind>,
120+
) -> Option<(String, String)> {
121+
if ref_allele.len() <= 1
122+
|| !matches!(kind, Some(VariantKind::Deletion | VariantKind::Insertion))
123+
{
124+
return None;
125+
}
126+
if !display
127+
.chars()
128+
.filter(char::is_ascii_alphabetic)
129+
.all(|allele| matches!(allele.to_ascii_uppercase(), 'I' | 'D'))
130+
{
131+
return None;
132+
}
133+
if matches!(kind, Some(VariantKind::Deletion)) {
134+
return Some(("I".to_owned(), "D".to_owned()));
135+
}
136+
if matches!(kind, Some(VariantKind::Insertion)) {
137+
return Some(("D".to_owned(), "I".to_owned()));
138+
}
139+
None
140+
}
141+
112142
fn normalize_long_allele_genotype(
113143
display: &str,
114144
ref_allele: &str,

0 commit comments

Comments
 (0)