Merge pull request #82 from uclahs-cds/nzeltser-add-strand-flip-smart-check

add strand flip smart check

Author: alkaZeltser

NEWS.md

@@ -11,6 +11,7 @@
 * Added minimum sample size check for grouped density curves
 * Added new plotting function `create.pgs.boxplot`
 * Added option for user to provide custom PGS source column(s) for plotting functions
+* Added option to `assess.pgs.vcf.allele.match` to condition the handling of ambiguous strand flips on the total number of unambiguous strand flips.
 
 # ApplyPolygenicScore 3.0.2
 

R/apply-pgs.R

@@ -100,6 +100,8 @@ validate.phenotype.data.input <- function(phenotype.data, phenotype.analysis.col
 #' The PGS catalog standard column \code{other_allele} in \code{pgs.weight.data} is required for this check.
 #' @param remove.ambiguous.allele.matches A logical indicating whether to remove PGS variants with ambiguous allele matches between PGS weight data and VCF genotype data. Default is \code{FALSE}.
 #' The PGS catalog standard column \code{other_allele} in \code{pgs.weight.data} is required for this check.
+#' @param max.strand.flips An integer indicating the number of unambiguous strand flips that need to be detected in order to discard all variants with ambiguous allele matches. Only applies if {return.ambiguous.as.missing == TRUE}.
+#' Default is \code{0} which means that all ambiguous variants are removed regardless of the status of any other variant.
 #' @param remove.mismatched.indels A logical indicating whether to remove indel variants that are mismatched between PGS weight data and VCF genotype data. Default is \code{FALSE}.
 #' The PGS catalog standard column \code{other_allele} in \code{pgs.weight.data} is required for this check.
 #' @param output.dir A character string indicating the directory to write output files. Separate files are written for per-sample pgs results and optional regression results.
@@ -198,6 +200,9 @@ validate.phenotype.data.input <- function(phenotype.data, phenotype.analysis.col
 #' \item \code{remove.ambiguous.allele.matches}: Corresponds to the \code{return.ambiguous.as.missing} argument in \code{assess.pgs.vcf.allele.match}. When \code{TRUE}, non-INDEL allele
 #' mismatches that cannot be resolved (due to palindromic alleles or causes other than strand flips) are removed by marking the affected value in the \code{effect_allele} column as missing
 #' prior to dosage calling and missing genotype handling. The corresponding dosage is set to NA and the variant is handled according to the chosen missing genotype method.
+#' \item \code{max.strand.flips}: This argument only applies when \code{remove.ambiguous.allele.matches} is on and modifies its behavior. In cases where none or very few unambiguous strand flips are detected,
+#' it is likely that all ambiguous allele matches are simply palindromic effect size flips. This option facilitates handling of ambiguous allele matches conditional on a maximum number of unambiguous strand flips.
+#' Variants with ambiguous strand flips will be marked as missing only if the number of unambiguous strand flips is greater than or equal to \code{max.strand.flips}.
 #' \item \code{remove.mismatched.indels}: Corresponds to the \code{return.indels.as.missing} argument in \code{assess.pgs.vcf.allele.match}. When \code{TRUE}, INDEL allele mismatches
 #' (which cannot be assessed for strand flips) are removed by marking the affected value in the \code{effect_allele} column as missing prior to dosage calling and missing genotype handling.
 #' The corresponding dosage is set to NA and the variant is handled according to the chosen missing genotype method.
@@ -280,6 +285,7 @@ apply.polygenic.score <- function(
     phenotype.analysis.columns = NULL,
     correct.strand.flips = TRUE,
     remove.ambiguous.allele.matches = FALSE,
+    max.strand.flips = 0,
     remove.mismatched.indels = FALSE,
     output.dir = NULL,
     file.prefix = NULL,
@@ -345,12 +351,16 @@ apply.polygenic.score <- function(
 
     ### Start Allele Match Check ###
     if (remove.ambiguous.allele.matches || correct.strand.flips) {
+        # adjust max.strand.flips threshold to account for the long-form of the variant x sample matrix
+        # each allele match appears once per sample, so the threshold is multiplied by the number of samples
+        adjusted.max.strand.flips <- max.strand.flips * length(unique(merged.vcf.with.pgs.data$Indiv));
         match.assessment <- ApplyPolygenicScore::assess.pgs.vcf.allele.match(
             vcf.ref.allele = merged.vcf.with.pgs.data$REF,
             vcf.alt.allele = merged.vcf.with.pgs.data$ALT,
             pgs.ref.allele = merged.vcf.with.pgs.data$other_allele,
             pgs.effect.allele = merged.vcf.with.pgs.data$effect_allele,
             return.ambiguous.as.missing = remove.ambiguous.allele.matches,
+            max.strand.flips = adjusted.max.strand.flips,
             return.indels.as.missing = remove.mismatched.indels
             );
         merged.vcf.with.pgs.data$effect_allele <- match.assessment$new.pgs.effect.allele;

R/assess-strand-flip.R

@@ -79,6 +79,8 @@ flip.DNA.allele <- function(alleles, return.indels.as.missing = FALSE) {
 #' @param return.indels.as.missing A logical value indicating whether to return NA for INDEL alleles with detected mismatches. Default is \code{FALSE}.
 #' @param return.ambiguous.as.missing A logical value indicating whether to return NA for ambiguous cases where both a strand flip and effect switch are possible,
 #' or no strand flip is detected and a mismatch cannot be resolved. Default is \code{FALSE}.
+#' @param max.strand.flips An integer indicating the number of non-ambiguous strand flips that must be present to implement the discarding all allele matches labeled "ambiguous_flip". Only applies if {return.ambiguous.as.missing == TRUE}.
+#' Defaults to \code{0}, meaning that no strand flips are allowed. Allele matches labeled "unresolved_mismatch" are not affected by this parameter.
 #' @return A list containing the match assessment, a new PGS effect allele, and a new PGS other allele.
 #'
 #' \strong{Output Structure}
@@ -88,7 +90,7 @@ flip.DNA.allele <- function(alleles, return.indels.as.missing = FALSE) {
 #' \item \code{match.status}: A character vector indicating the match status for each pair of allele pairs. Possible values are \code{default_match}, \code{effect_switch}, \code{strand_flip}, \code{effect_switch_with_strand_flip}, \code{ambiguous_flip}, \code{indel_mismatch}, and \code{unresolved_mismatch}.
 #' \item \code{new.pgs.effect.allele}: A character vector of new PGS effect alleles based on the match status. If the match status is \code{default_match}, \code{effect_switch} or \code{missing_allele}, the original PGS effect allele is returned.
 #' If the match status is \code{strand_flip} or \code{effect_switch_with_strand_flip} the flipped PGS effect allele is returned. If the match status is \code{ambiguous_flip}, \code{indel_mismatch}, or \code{unresolved_mismatch},
-#' the return value is either the original allele or NA as dictated by the \code{return.indels.as.missing} and \code{return.ambiguous.as.missing} parameters.
+#' the return value is either the original allele or NA as dictated by the \code{return.indels.as.missing}, \code{return.ambiguous.as.missing}, and \code{max.strand.flips} parameters.
 #' \item \code{new.pgs.other.allele}: A character vector of new PGS other alleles based on the match status, following the same logic as \code{new.pgs.effect.allele}.
 #' }
 #'
@@ -123,7 +125,8 @@ assess.pgs.vcf.allele.match <- function(
     pgs.ref.allele,
     pgs.effect.allele,
     return.indels.as.missing = FALSE,
-    return.ambiguous.as.missing = FALSE
+    return.ambiguous.as.missing = FALSE,
+    max.strand.flips = 0
     ) {
 
     # check that all inputs are one dimensional character vectors
@@ -144,6 +147,20 @@ assess.pgs.vcf.allele.match <- function(
         stop('vcf.ref.allele, vcf.alt.allele, pgs.ref.allele, and pgs.effect.allele must be the same length.');
         }
 
+    # check that all logical parameters are TRUE or FALSE
+    if (!is.logical(return.indels.as.missing) || length(return.indels.as.missing) != 1) {
+        stop('return.indels.as.missing must be a single logical value.');
+        }
+
+    if (!is.logical(return.ambiguous.as.missing) || length(return.ambiguous.as.missing) != 1) {
+        stop('return.ambiguous.as.missing must be a single logical value.');
+        }
+
+    # check that max.strand.flips is a single integer
+    if (!is.numeric(max.strand.flips) || length(max.strand.flips) != 1 || max.strand.flips < 0 || max.strand.flips != floor(max.strand.flips)) {
+        stop('max.strand.flips must be a single non-negative integer.');
+        }
+
     # verify valid alleles
     validate.allele.input(vcf.ref.allele, na.allowed = TRUE);
     validate.allele.input(pgs.ref.allele, na.allowed = TRUE);
@@ -283,16 +300,11 @@ assess.pgs.vcf.allele.match <- function(
         if (effect.switch.candidate && strand.flip.candidate) {
             # not possible to determine whether effect switch or strand flip has occured
             # This is an ambiguous case caused by palindromic SNPs
+            # Default is to assume that this is an effect switch until strand flip threshold is exceeded
             flip.designation[i] <- 'ambiguous_flip';
-            if (return.ambiguous.as.missing) {
-                flipped.effect.allele[i] <- NA;
-                flipped.other.allele[i] <- NA;
-                next;
-                } else {
-                flipped.effect.allele[i] <- current.pgs.effect.allele;
-                flipped.other.allele[i] <- current.pgs.ref.allele;
-                next;
-                }
+            flipped.effect.allele[i] <- current.pgs.effect.allele;
+            flipped.other.allele[i] <- current.pgs.ref.allele;
+            next;
             } else if (effect.switch.candidate) {
             # if this is a clear-cut effect switch, return the default PGS alleles
             # apply.polygenic.score automatically handles effect switches during PGS application
@@ -338,6 +350,26 @@ assess.pgs.vcf.allele.match <- function(
 
         }
 
+    # implement max strand flip threshold for ambiguous cases
+    # count number of unambiguous strand flips
+    if (return.ambiguous.as.missing) {
+        if (max.strand.flips > 0) {
+            # if a max strand flips threshold is set, check if it is exceeded
+            unambiguous.strand.flips <- sum(flip.designation == 'strand_flip' | flip.designation == 'effect_switch_with_strand_flip');
+            if (unambiguous.strand.flips >= max.strand.flips) {
+                # if the number of unambiguous strand flips equals or exceeds the threshold, return NA for all ambiguous cases
+                flipped.effect.allele[flip.designation == 'ambiguous_flip'] <- NA;
+                flipped.other.allele[flip.designation == 'ambiguous_flip'] <- NA;
+                }
+                # if max threshold is not exceeded, return all ambiguous cases with no flips (default)
+            } else {
+            # if max.strand.flips is 0 there is no tolerance for ambiguous cases, return all as NA
+            flipped.effect.allele[flip.designation == 'ambiguous_flip'] <- NA;
+            flipped.other.allele[flip.designation == 'ambiguous_flip'] <- NA;
+            }
+        # If return.ambiguous.as.missing is FALSE, keep all ambiguous cases as is (no flips)
+        }
+
     return(
         list(
             match.status = flip.designation,

man/apply.polygenic.score.Rd

@@ -737,6 +737,7 @@ test_that(
 test_that(
     'apply.polygenic.score correctly handles strand flipping', {
         load('data/strand.flip.test.data.Rda');
+        # no strand flip correction
         strand.flip.raw.data <- apply.polygenic.score(
             vcf.data = strand.flip.test.data$strand.flip.vcf.data,
             pgs.weight.data = strand.flip.test.data$strand.flip.pgs.weight.data,
@@ -754,6 +755,7 @@ test_that(
             c(0, 0, 0)
             );
 
+        # strand flip correction
         strand.flip.correct.flips <- apply.polygenic.score(
             vcf.data = strand.flip.test.data$strand.flip.vcf.data,
             pgs.weight.data = strand.flip.test.data$strand.flip.pgs.weight.data,
@@ -771,6 +773,7 @@ test_that(
             c(0, 0, 0)
             );
 
+        # ambiguous allele matches removed
         strand.flip.remove.ambiguous <- apply.polygenic.score(
             vcf.data = strand.flip.test.data$strand.flip.vcf.data,
             pgs.weight.data = strand.flip.test.data$strand.flip.pgs.weight.data,
@@ -789,6 +792,28 @@ test_that(
             c(2, 2, 2)
             );
 
+        # ambiguous allele matches not removed due to threshold
+        strand.flip.threshold.ambiguous <- apply.polygenic.score(
+            vcf.data = strand.flip.test.data$strand.flip.vcf.data,
+            pgs.weight.data = strand.flip.test.data$strand.flip.pgs.weight.data,
+            correct.strand.flips = TRUE,
+            missing.genotype.method = c('none', 'mean.dosage', 'normalize'),
+            remove.ambiguous.allele.matches = TRUE,
+            max.strand.flips = 3 # there are two unambiguous strand flips in the test data, setting threshold above that
+            );
+
+        # expect equal to not removing ambiguous allele matches
+        expect_equal(
+            strand.flip.threshold.ambiguous$pgs.output$PGS, # dosage is unaffected because unresolved mismatch does not contain an effect allele so dosage is 0
+            strand.flip.correct.flips$pgs.output$PGS
+            );
+
+        expect_equal(
+            strand.flip.threshold.ambiguous$pgs.output$n.missing.genotypes,
+            c(1, 1, 1) # the unresolved mismatch is marked missing
+            );
+
+        # indels removed
         strand.flip.remove.indel.mismatches <- apply.polygenic.score(
             vcf.data = strand.flip.test.data$strand.flip.vcf.data,
             pgs.weight.data = strand.flip.test.data$strand.flip.pgs.weight.data,