Allow IUPAC codes on Match alleles - #67
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## development #67 +/- ##
===============================================
- Coverage 81.11% 80.92% -0.19%
===============================================
Files 19 19
Lines 2277 2349 +72
===============================================
+ Hits 1847 1901 +54
- Misses 430 448 +18 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Cristianetaniguti
added
enhancement
There was a problem hiding this comment.
Pull request overview
This PR updates BIGr’s MADC sanity checks and VCF conversion pipeline to allow/ignore IUPAC ambiguity codes in Match alleles (and to treat some Ref/Alt IUPAC cases differently), aligning BIGr behavior with HapApp’s current IUPAC-cleaning scope.
Changes:
- Expanded
check_madc_sanity()to distinguish (a) problematic Ref/Alt IUPAC differences vs (b) identical-position Ref/Alt IUPAC vs (c) IUPAC present in RefMatch/AltMatch alleles. - Updated
madc2vcf_all()to pass through additional “Other allele” filtering knobs (others_min_dist,others_max_close_snps) and adjust comparison logic to ignore non-ATCG bases in Match/Other alleles. - Updated tests and docs (Rd + NEWS) to reflect the new checks/behavior and revised error messages.
Reviewed changes
Copilot reviewed 5 out of 7 changed files in this pull request and generated 7 comments.
Show a summary per file
| File | Description |
|---|---|
R/madc2vcf_all.R |
Adds new “Other allele” filtering parameters, extends compare logic, and broadens markers_info ID column support. |
R/check_madc_sanity.R |
Splits IUPAC detection into multiple checks (Ref/Alt differing vs identical-position vs Match alleles) and updates messages/return shape. |
tests/testthat/test-madc2vcf_all.R |
Adjusts expectations around warnings and updates IUPAC-related error regex. |
tests/testthat/test-check_madc_sanity.R |
Updates expected check vector to match the expanded set of checks. |
man/madc2vcf_all.Rd |
Documents new arguments and default changes. |
man/check_madc_sanity.Rd |
Documents the expanded check set and new returned element. |
NEWS.md |
Adds release notes for the new IUPAC handling and “Other allele” filtering options. |
Files not reviewed (2)
- man/check_madc_sanity.Rd: Generated file
- man/madc2vcf_all.Rd: Generated file
Comments suppressed due to low confidence (3)
R/madc2vcf_all.R:175
markers_infonow supportsMarker_ID, but the validation error messages still only mentionCloneID/BI_markerID, which is misleading. Also, theif(id_col == ... | id_col == ...)check is scalar and would be clearer/safer as%in%(and avoids the vectorized|).
id_cols <- intersect(c("CloneID", "BI_markerID","Marker_ID"), colnames(mi_df))
if(!length(id_cols)) {
stop("markers_info must contain a marker ID column named either 'CloneID' or 'BI_markerID'.")
}
match_n <- vapply(id_cols, function(col) {
tests/testthat/test-madc2vcf_all.R:343
- This
expect_warning()call does not assert the expected warning message/class, so it can pass even if an unrelated warning starts being emitted. Please add aregexpmatching the warning you expect here.
expect_warning(madc2vcf_all(madc = alfalfa_lowercase,
botloci_file = alfalfa_botloci,
hap_seq_file = NULL,
n.cores = 1,
out_vcf = out,
tests/testthat/test-madc2vcf_all.R:363
- This
expect_warning()call does not assert the expected warning message/class, so it can pass even if an unrelated warning starts being emitted. Please add aregexpmatching the warning you expect here.
expect_warning(madc2vcf_all(madc = alfalfa_lowercase,
botloci_file = alfalfa_botloci,
hap_seq_file = NULL,
n.cores = 1,
markers_info = alfalfa_markers_info,
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 7 out of 9 changed files in this pull request and generated no new comments.
Files not reviewed (2)
- man/check_madc_sanity.Rd: Generated file
- man/madc2vcf_all.Rd: Generated file
Comments suppressed due to low confidence (7)
R/check_madc_sanity.R:159
IUPACcodes_IdenticalRefAltcurrently treats Ref/Alt as “identical IUPAC” when both alleles have any IUPAC code at the same positions, even if the actual ambiguity codes differ (e.g.,RvsY). This can cause differing Ref/Alt IUPAC codes to be misclassified as identical and bypass the intended stop condition.
# Check if IUPAC codes are at the same positions
ref_is_iupac <- !ref_chars %in% c("A", "T", "C", "G", "-")
alt_is_iupac <- !alt_chars %in% c("A", "T", "C", "G", "-")
# IUPAC codes must be at exactly the same positions
all(ref_is_iupac == alt_is_iupac)
tests/testthat/test-madc2vcf_all.R:365
- This test now accepts any warning from
madc2vcf_all()(noregexp), which can mask unrelated warnings and make failures harder to diagnose. It’s better to assert the expected warning text (or explicitly silence specific warnings).
expect_warning(madc2vcf_all(madc = alfalfa_lowercase,
botloci_file = alfalfa_botloci,
hap_seq_file = NULL,
n.cores = 1,
markers_info = alfalfa_markers_info,
out_vcf = out,
verbose = FALSE))
tests/testthat/test-madc2vcf_all.R:55
- This test now accepts any warning from
madc2vcf_all()(noregexp), which can mask unrelated warnings and make failures harder to diagnose. It’s better to assert the expected warning text (or explicitly silence specific warnings).
expect_warning(madc2vcf_all(madc = madc_file,
botloci_file = bot_file,
hap_seq_file = NULL,
n.cores = 2,
rm_multiallelic_SNP = FALSE,
tests/testthat/test-madc2vcf_all.R:344
- This test now accepts any warning from
madc2vcf_all()(noregexp), which can mask unrelated warnings and make failures harder to diagnose. It’s better to assert the expected warning text (or explicitly silence specific warnings).
expect_warning(madc2vcf_all(madc = alfalfa_lowercase,
botloci_file = alfalfa_botloci,
hap_seq_file = NULL,
n.cores = 1,
out_vcf = out,
verbose = FALSE))
R/check_madc_sanity.R:285
- Grammar in the new
IUPACcodes_MatchAllelesmessage: “This codes” → “These codes”, and add punctuation for readability (this string is surfaced to users viacheck_madc_sanity).
messages[["IUPACcodes_MatchAlleles"]] <- c("IUPAC (non-ATCG) codes found in RefMatch/AltMatch AlleleSequence. This codes are not currently supported by BIGr/BIGapp they will be ignored in the conversion to VCF",
"No IUPAC (non-ATCG) codes found in RefMatch/AltMatch AlleleSequence")
R/madc2vcf_all.R:704
- The “close SNPs” filter counts adjacent pairs of SNPs closer than
others_min_dist(sum(diff(pos_ref_idx) < ...)), but the parameter/message/NEWS describe a threshold in terms of SNPs that are close to another SNP. This can undercount by 1 (e.g., 4 clustered SNPs produce 3 close gaps) and keep alleles that should be discarded per the documented behavior.
# Discard Others with SNPs too close to each other
## Defined others_min_dist (default 5)
## Defined others_max_close_snps (default 3)
if(!is.null(others_min_dist) &&
!is.null(others_max_close_snps) &&
any(diff(pos_ref_idx) < others_min_dist)) {
if(sum(diff(pos_ref_idx) < others_min_dist) > others_max_close_snps) {
n_rm_others_close_snps <- n_rm_others_close_snps + 1L
next
}
}
R/madc2vcf_all.R:173
markers_infonow supportsMarker_ID, but this error message still only mentionsCloneIDandBI_markerID, which can mislead users when validation fails.
stop("markers_info must contain a marker ID column named either 'CloneID' or 'BI_markerID'.")
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 8 out of 10 changed files in this pull request and generated no new comments.
Files not reviewed (2)
- man/check_madc_sanity.Rd: Generated file
- man/madc2vcf_all.Rd: Generated file
Comments suppressed due to low confidence (3)
R/madc2vcf_all.R:648
- Match-allele polymorphisms at positions where the Ref has an IUPAC base (allowed when Ref/Alt share identical IUPAC) are still being reported because filtering only excludes non-ATCG bases in the Match sequence. This contradicts the intent to ignore polymorphisms at IUPAC positions.
# If Match sequences have N or IUPAC codes, do not consider the polymorphism
if(any(!alt_base_match %in% c("A", "T", "C", "G"))) {
ref_base_match <- ref_base_match[-which(!alt_base_match %in% c("A", "T", "C", "G"))]
pos_ref_idx <- pos_ref_idx[-which(!alt_base_match %in% c("A", "T", "C", "G"))]
alt_base_match <- alt_base_match[-which(!alt_base_match %in% c("A", "T", "C", "G"))]
}
R/madc2vcf_all.R:174
- The error message for missing marker ID columns doesn’t mention the newly supported
Marker_IDcolumn, which can confuse users whenMarker_IDis present but the message still claims onlyCloneID/BI_markerIDare allowed.
id_cols <- intersect(c("CloneID", "BI_markerID","Marker_ID"), colnames(mi_df))
if(!length(id_cols)) {
stop("markers_info must contain a marker ID column named either 'CloneID' or 'BI_markerID'.")
}
R/madc2vcf_all.R:704
- The close-SNP filter counts the number of adjacent gaps with distance <
others_min_dist(i.e.,sum(diff(pos_ref_idx) < ...)). That undercounts the number of SNPs in a dense cluster (e.g., 4 SNPs in a row produces 3 gaps), so alleles can be kept even when they exceedothers_max_close_snpsper the documentation/NEWS text.
# Discard Others with SNPs too close to each other
## Defined others_min_dist (default 5)
## Defined others_max_close_snps (default 3)
if(!is.null(others_min_dist) &&
!is.null(others_max_close_snps) &&
any(diff(pos_ref_idx) < others_min_dist)) {
if(sum(diff(pos_ref_idx) < others_min_dist) > others_max_close_snps) {
n_rm_others_close_snps <- n_rm_others_close_snps + 1L
next
}
}
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 8 out of 10 changed files in this pull request and generated no new comments.
Files not reviewed (2)
- man/check_madc_sanity.Rd: Generated file
- man/madc2vcf_all.Rd: Generated file
Suppressed comments (5)
R/madc2vcf_all.R:737
- Same issue as above for Other allele mismatches: the
if (...) { ... } else if (...) { ... }block only filters non-ATCG bases in the reference when there are no non-ATCG bases in the Other allele at all. Ifother_alt_basecontains any IUPAC/N, positions whereother_ref_baseis IUPAC/N (from allowed identical REF/ALT IUPAC) won't be filtered and may be reported as polymorphisms. Filter positions where either base is non-ATCG.
if(any(!other_alt_base %in% c("A", "T", "C", "G"))) {
rm_pos <- which(!other_alt_base %in% c("A", "T", "C", "G"))
other_ref_base <- other_ref_base[-rm_pos]
pos_ref_idx <- pos_ref_idx[-rm_pos]
other_alt_base <- other_alt_base[-rm_pos]
# Check if the reference sequence doesn't have N or IUPAC codes
# all previous checks will let pass identical IUPAC codes between REF and ALT
} else if(any(!other_ref_base %in% c("A", "T", "C", "G"))){
rm_pos <- which(!other_ref_base %in% c("A", "T", "C", "G"))
other_ref_base <- other_ref_base[-rm_pos]
pos_ref_idx <- pos_ref_idx[-rm_pos]
other_alt_base <- other_alt_base[-rm_pos]
R/check_madc_sanity.R:286
- Grammar in this user-facing message: use “These codes” (not “This codes”) and add punctuation so the warning reads clearly.
messages[["IUPACcodes_MatchAlleles"]] <- c("IUPAC (non-ATCG) codes found in RefMatch/AltMatch/Other AlleleSequence. This codes are not currently supported by BIGr/BIGapp they will be ignored in the conversion to VCF",
"No IUPAC (non-ATCG) codes found in RefMatch/AltMatch/Other AlleleSequence")
R/madc2vcf_all.R:652
- The IUPAC/N filtering for Match allele mismatches only checks
ref_base_matchfor non-ATCG bases when none ofalt_base_matchare non-ATCG (due to theif (...) { ... } else if (...) { ... }structure). Ifalt_base_matchcontains at least one IUPAC/N, mismatches where the REF base is IUPAC/N at other positions will not be filtered and can be incorrectly treated as polymorphisms. Filter out positions where either REF or Match base is non-ATCG.
This issue also appears on line 726 of the same file.
alt_base_match <- alt_base_match[-which(!alt_base_match %in% c("A", "T", "C", "G"))]
} else if(any(!ref_base_match %in% c("A", "T", "C", "G"))){
alt_base_match <- alt_base_match[-which(!ref_base_match %in% c("A", "T", "C", "G"))]
pos_ref_idx <- pos_ref_idx[-which(!ref_base_match %in% c("A", "T", "C", "G"))]
ref_base_match <- ref_base_match[-which(!ref_base_match %in% c("A", "T", "C", "G"))]
R/madc2vcf_all.R:186
- Use
||/%in%instead of the vectorized|here sinceid_colis scalar. This avoids relying on vector recycling and reads more clearly.
if(id_col == "BI_markerID" | id_col == "Marker_ID") {
vmsg("markers_info: 'BI_markerID' or 'Marker_ID' column copied to 'CloneID' for internal use", verbose = verbose, level = 1)
R/madc2vcf_all.R:636
- Typo/grammar in comment text: “discarted, they should have be” should be “discarded; they should have been”.
# Cases found where the AltMatch is another alternative for the target SNP - they are discarted, they should have be named Other
alex-sandercock
left a comment
alex-sandercock
left a comment
There was a problem hiding this comment.
I added some in-line comments and some file comments on madc2vcf_all()
| loop_though_dartag_report <- function(report, | ||
| botloci, | ||
| hap_seq, | ||
| n.cores=1, |
There was a problem hiding this comment.
n.cores = 1 is hardcoded here despite being a user option. Is there a reason for this, or should this be n.cores = n.cores?
There was a problem hiding this comment.
@alex-sandercock it should be fine. This line just define the default of the function, what really matters is this one here:
https://github.com/Breeding-Insight/BIGr/blame/2994469818454470758f01fc090d91e6896c2e49/R/madc2vcf_all.R#L259
There was a problem hiding this comment.
Noting some review findings:
-
[P2] Exclude IUPAC mismatches before applying Other filters — /Users/aherranssanderco/UF Dropbox/Alexander Sandercock/AS_Projects/Pipelines/BIGr/R/madc2vcf_all.R:702-704
When an Other sequence contains allowed IUPAC differences, pos_ref_idx includes those ambiguity mismatches and this proximity filter runs before lines 730–735 remove
them. Consequently, four nearby IUPAC mismatches plus one valid distant SNP causes the entire allele to be discarded even though the IUPAC positions are supposed to be
ignored; filter non-ATCG positions before both the proximity and maximum-SNP checks. -
[P2] Sanitize IUPAC alleles before passing data to polyRAD — /Users/aherranssanderco/UF Dropbox/Alexander Sandercock/AS_Projects/Pipelines/BIGr/R/madc2vcf_multi.R:91-92
For reports with identical IUPAC codes in Ref/Alt or IUPAC codes in Match/Other rows, the revised sanity check leaves IUPACcodes false, so this guard permits the report
unchanged. readDArTag retains those complete haplotype sequences and RADdata2VCF(..., asSNPs = FALSE) can emit ambiguity codes such as R or Y in REF/ALT rather than
ignoring them, producing unsupported VCF alleles; either reject these cases for this converter or sanitize them before loading polyRAD.
Additional review comments:
Cleaned up. Here's the review.
Chapter: Review findings
Setup note: the PR's real base is development, not main — git diff main...HEAD pulls in the whole grape/BIGpopA-removal merge and is misleading. Against development the change is 10 files, +289/−100.
Verified good: full suite passes (95 assertions in test-madc2vcf_all.R, 18 in test-check_madc_sanity.R, 0 failures across test_dir), man/ and NAMESPACE are exactly what roxygen regenerates, version bump and NEWS heading agree. The core mechanism works — I confirmed with pwalign that identical IUPAC at identical positions is not reported as a mismatch (so target-SNP detection is unaffected), and that the new rm_pos filter correctly suppresses the false off-target SNPs those positions would otherwise generate.
Note: pkgload::load_all() + parLapply will fail here — the cluster workers resolve compare() from the installed BIGr, not the dev namespace. I installed to a temp lib to test.
1. The Other-allele filters count the very positions the PR declares non-polymorphic
Both others_max_snps (madc2vcf_all.R:694) and the new close-SNP filter (madc2vcf_all.R:702) run on raw pos_ref_idx, before the non-ATCG filtering at madc2vcf_all.R:730. A masked position in Ref vs a resolved base in the Other allele is a mismatch, so it counts toward both filters — then gets thrown away as "not a polymorphism" a few lines later.
Two synthetic tags, each with exactly one genuine off-target SNP:
| tag | result with defaults |
|---|---|
| 4 masked positions 2bp apart + 1 real SNP | discarded, rm_close = 1 |
| 5 masked positions >5bp apart + 1 real SNP | discarded, rm_maxsnps = 1 |
With the filter disabled, the same tag yields the one correct off-target SNP. This bites exactly the VariantMasked MADCs the PR targets (the debug comment at madc2vcf_all.R:539 names one). Fix: prune non-ATCG positions from pos_ref_idx/pos_alt_idx right after they're computed at line 690–691, then apply both filters.
2. madc2vcf_multi now writes IUPAC codes into VCF REF/ALT
Relaxing the stop at madc2vcf_multi.R:91 lets identical-IUPAC files through, but unlike madc2vcf_all this path never resolves single bases — polyRAD::RADdata2VCF(asSNPs = FALSE) writes whole haplotype sequences. I injected an identical R into one Ref/Alt pair of alfalfa_madc.csv; development rejects the file, this branch succeeds and emits:
chr7.1 54162026 REF: GGTARTTTTTTCTTCAATTGTGACTTTGCTCGG...
ALT: GGTARTTTTTTCTTCAATTGTGACTTTACTCGG...
VCF 4.3 restricts REF/ALT bases to A,C,G,T,N. Either keep the stop for madc2vcf_multi (gate on IUPACcodes | IUPACcodes_IdenticalRefAlt), or replace IUPAC with N before handing off to polyRAD.
3. IUPAC detection silently lost for non-HapApp MADCs
check_madc_sanity.R:130 selects with grepl("Ref_", AlleleID) / grepl("Alt_", AlleleID), but unfixed MADCs use CloneID|Ref / CloneID|Alt — the function's own RefAltSeqs check handles both forms at lines 241–246. A |Ref/|Alt pair containing an N returns IUPACcodes = FALSE; the old check scanned every AlleleSequence. All three converters demand fixed IDs before inferring REF/ALT, so nothing downstream breaks today — but check_madc_sanity is exported and its documented contract is now wrong for raw files.
4. Neither new feature has a test
Every one of the 8 fixtures in test-check_madc_sanity.R expects IUPACcodes_IdenticalRefAlt = FALSE, so the headline pathway (and iupac_identical_clone_ids) is never exercised. Nothing touches others_min_dist/others_max_close_snps or the new Marker_ID column. My probes above are basically the missing tests.
Minor
others_min_distandothers_max_close_snpsskip the validation block at madc2vcf_all.R:127 that guardsothers_max_snps;others_max_close_snps = "3"silently string-compares. They're also absent from the verbose parameter dump at lines 105–108.others_min_distdocs don't mention thatNULLdisables the filter, though the code checks for it.- The
expect_warning()wrappers added intest-madc2vcf_all.Rhave noregexp. I checked — the warning ("1 tags discarded due to lack of Alt_0002 sequence") already fires ondevelopment, so this is just tightening, not new behavior. But bareexpect_warning()will absorb any future unrelated warning; worth pinning the message. - NEWS.md line 8 is a bare paragraph that breaks the bullet list;
“Other“has two opening curly quotes; "polymorphims" typo on line 9. @detailsforIUPACcodes_IdenticalRefAltsays "at exactly the same positions", but check_madc_sanity.R:159 also requires the codes themselves to match (which the inline comment states correctly).
Verdict: the design is sound and the mechanism does what it claims, but I'd hold the PR on #1 and #2 — both make the feature misfire on the input class it was built for. #3 and #4 are worth folding in while you're here.
There was a problem hiding this comment.
About:
[P2] Exclude IUPAC mismatches before applying Other filters — /Users/aherranssanderco/UF Dropbox/Alexander Sandercock/AS_Projects/Pipelines/BIGr/R/madc2vcf_all.R:702-704
When an Other sequence contains allowed IUPAC differences, pos_ref_idx includes those ambiguity mismatches and this proximity filter runs before lines 730–735 remove
them. Consequently, four nearby IUPAC mismatches plus one valid distant SNP causes the entire allele to be discarded even though the IUPAC positions are supposed to be
ignored; filter non-ATCG positions before both the proximity and maximum-SNP checks.
co-pilot pointed that previously, I thought about it and decided to keep current behavior to be more conservative. If IUPAC code is present, It think it means that the base has some "issue", like a variant or that sequencing found difficulty. It shows a regions instability, we probably don't want to have a SNPs close to it.
Same observation is made in The Other-allele filters count the very positions the PR declares non-polymorphic
There was a problem hiding this comment.
About:
[P2] Sanitize IUPAC alleles before passing data to polyRAD — /Users/aherranssanderco/UF Dropbox/Alexander Sandercock/AS_Projects/Pipelines/BIGr/R/madc2vcf_multi.R:91-92
For reports with identical IUPAC codes in Ref/Alt or IUPAC codes in Match/Other rows, the revised sanity check leaves IUPACcodes false, so this guard permits the report
unchanged. readDArTag retains those complete haplotype sequences and RADdata2VCF(..., asSNPs = FALSE) can emit ambiguity codes such as R or Y in REF/ALT rather than
ignoring them, producing unsupported VCF alleles; either reject these cases for this converter or sanitize them before loading polyRAD.
I tested the grape MADC that has this scenario, when both REF and ALT has a R in the same position. polyRAD doesn't complain about it, so the sanity check doesn't block it.
Same observation reported in madc2vcf_multi now writes IUPAC codes into VCF REF/ALT
There was a problem hiding this comment.
About 3. IUPAC detection silently lost for non-HapApp MADCs
The only case non-HapApp MADCs is being processed in with Targets SNPs and it requires the markers info files to recover Positions, REF and ALT. Therefore, the IUPAC codes won't do any harm to it.
There was a problem hiding this comment.
About 4. Neither new feature has a test
I need to make a simulated MADC to test the new features. I plan to do it later card SAT-189
There was a problem hiding this comment.
Minor comments were addressed in the most recent commit
34c1f7e to
a851de7
Compare
HapApp pipelines only replaces IUPAC codes present in the Ref and Alt alleles. IUPAC codes present in Match alleles should be ignored by the BIGr SNP calling algorithms.