add strand flip unit tests

Author: alkaZeltser

tests/testthat/test-pgs-application.R

@@ -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,

tests/testthat/test-strand-flip-handling.R

@@ -442,6 +442,61 @@ test_that(
             c('A', NA)
             );
 
+        # palindromic case, return.ambiguous.as.missing == TRUE, max.strand.flips = 1
+        # expected output: all ambiguous returned as-is due to lack of unambiguous strand flips
+
+        test.ambiguous.below.threshold <- assess.pgs.vcf.allele.match(
+            c('A', 'A'),
+            c('G', 'T'),
+            c('A', 'T'),
+            c('G', 'A'),
+            return.ambiguous.as.missing = TRUE,
+            max.strand.flips = 1
+            );
+
+        expect_equal(
+            test.ambiguous.below.threshold$match.status,
+            c('default_match', 'ambiguous_flip')
+            );
+
+        expect_equal(
+            test.ambiguous.below.threshold$new.pgs.effect.allele,
+            c('G', 'A')
+            );
+
+        expect_equal(
+            test.ambiguous.below.threshold$new.pgs.other.allele,
+            c('A', 'T')
+            );
+
+        # palindromic case, return.ambiguous.as.missing == TRUE, max.strand.flips = 2, strand flips exceed threshold
+        # expected output: all ambiguous returned as missing due to number of unambiguous strand flips exceeding threshold
+
+        test.ambiguous.above.threshold <- assess.pgs.vcf.allele.match(
+            # last three alleles are unambiguous strand flips
+            c('A', 'A', 'A', 'A', 'A'),
+            c('G', 'T', 'G', 'C', 'G'),
+            c('A', 'T', 'T', 'T', 'T'),
+            c('G', 'A', 'C', 'G', 'C'),
+            return.ambiguous.as.missing = TRUE,
+            max.strand.flips = 2
+            );
+
+        expect_equal(
+            test.ambiguous.above.threshold$match.status,
+            c('default_match', 'ambiguous_flip', 'strand_flip', 'strand_flip', 'strand_flip')
+            );
+
+        expect_equal(
+            test.ambiguous.above.threshold$new.pgs.effect.allele,
+            c('G', NA, 'G', 'C', 'G')
+            );
+
+        expect_equal(
+            test.ambiguous.above.threshold$new.pgs.other.allele,
+            c('A', NA, 'A', 'A', 'A')
+            );
+
         # unresolved case, return.ambiguous.as.missing == TRUE
         test.unresolved.missing <- assess.pgs.vcf.allele.match(
             c('A', 'A'),