Merge pull request #70 from uclahs-cds/nzeltser-update-strand-flip

small functionality patches to strand flip checks, dosage calculator, automated statistics

Author: alkaZeltser

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: ApplyPolygenicScore
 Type: Package
 Title: Utilities for the Application of a Polygenic Score to a VCF
-Version: 3.0.0
+Version: 3.0.1
 Authors@R: c(
     person('Paul', 'Boutros', role = 'cre', email = '[email protected]'),
     person('Nicole', 'Zeltser', role = 'aut', comment = c(ORCID = '0000-0001-7246-2771')),

NEWS.md

@@ -1,5 +1,14 @@
 # Unreleased
 
+# ApplyPolygenicScore 3.0.1
+
+## Added
+* Added hemizygous allele handling to dosage calculation
+
+## Changed
+* Updated INDEL effect switch reporting by strand flip checker
+* Updated data structuring for automated statistical analysis in `apply.polygenic.score`
+
 # ApplyPolygenicScore 3.0.0 (2024-12-02)
 
 ## Added

R/apply-pgs.R

@@ -561,9 +561,12 @@ apply.polygenic.score <- function(
         ### Begin Phenotype Analysis ###
 
         if (!is.null(phenotype.analysis.columns)) {
+            # post merge data selection
+            pgs.for.phenotype.stats <- pgs.output[ , missing.method.to.colname.ref[analysis.source.pgs]];
+
             regression.output <- run.pgs.regression(
-                pgs = pgs.for.stats,
-                phenotype.data = subset(phenotype.data, select = phenotype.analysis.columns)
+                pgs = pgs.for.phenotype.stats,
+                phenotype.data = subset(pgs.output, select = phenotype.analysis.columns)
                 );
             }
         ### End Phenotype Analysis ###

R/assess-strand-flip.R

@@ -216,6 +216,14 @@ assess.pgs.vcf.allele.match <- function(
         # identify insertion/deletion alleles in PGS
         if (nchar(current.pgs.ref.allele) > 1 || nchar(current.pgs.effect.allele) > 1) {
             pgs.indel <- TRUE;
+            if (effect.switch.candidate) {
+                # if an INDEL effect switch is detected, do not continue to strand flip assessment
+                # and return effect switch designation
+                flip.designation[i] <- 'effect_switch';
+                flipped.effect.allele[i] <- current.pgs.effect.allele;
+                flipped.other.allele[i] <- current.pgs.ref.allele;
+                next;
+                }
             # no INDEL flipping supported, return either NA or the original allele
             warning('Mismatch detected in INDEL PGS allele. Skipping strand flip assessment.');
             flip.designation[i] <- 'indel_mismatch';
@@ -234,6 +242,14 @@ assess.pgs.vcf.allele.match <- function(
 
         # identify insertion/deletion alleles in VCF
         if (nchar(current.vcf.ref.allele) > 1 || all(nchar(current.vcf.alt.allele.split) > 1)) {
+            if (effect.switch.candidate) {
+                # if an INDEL effect switch is detected, do not continue to strand flip assessment
+                # and return effect switch designation
+                flip.designation[i] <- 'effect_switch';
+                flipped.effect.allele[i] <- current.pgs.effect.allele;
+                flipped.other.allele[i] <- current.pgs.ref.allele;
+                next;
+                }
             # no INDEL flipping supported, return either NA or the original allele
             warning('Mismatch detected in INDEL VCF allele. Skipping strand flip assessment.');
             flip.designation[i] <- 'indel_mismatch';

R/calculate-dosage.R

@@ -2,7 +2,7 @@
 #' @description Convert genotype calls in the form of witten out alleles (e.g. 'A/T') to dosages (0, 1, 2) based on provided risk alleles from a PGS.
 #' @param called.alleles A vector of genotypes in allelic notation separated by a slash or pipe.
 #' @param risk.alleles A vector of risk alleles from a polygenic score corresponding to each genotype (by locus) in called.alleles.
-#' @return A vector of dosages corresponding to each genotype in called.alleles.
+#' @return A vector of dosages corresponding to each genotype in called.alleles. Hemizygous genotypes (one allele e.g. 'A') are counted as 1.
 #' @examples
 #' called.alleles <- c('A/A', 'A/T', 'T/T');
 #' risk.alleles <- c('T', 'T', 'T');
@@ -31,12 +31,16 @@ convert.alleles.to.pgs.dosage <- function(called.alleles, risk.alleles) {
         } else {
             # check that called.alleles is a vector of genotypes in allelic notation or '.' separated by a slash or pipe
             # "*" characters represent overlapping deletions from an upstream indel and are accepted VCF format
-            allowed.pattern <- '^((([A-Z]+|\\.|\\*)[/\\|]([A-Z]+|\\.|\\*))|\\.)$' # '|' are special chars in regular expressions
+            allowed.pattern <- '^((([A-Z]+|\\.|\\*)[/\\|]([A-Z]+|\\.|\\*))|\\.|[A-Z]+)$' # '|' are special chars in regular expressions
             passing.alleles <- grepl(allowed.pattern, called.alleles);
             passing.alleles[is.na(called.alleles)] <- TRUE; # NA allowed
             if (!all(passing.alleles)) {
                 stop('unrecognized called.alleles format, must be capitalized letters, "." or "*" separated by a slash or pipe.');
                 }
+            # replace hemizygous genotypes with a placeholder for easier splitting
+            # index for non-NA alleles that are missing allele separators:
+            no.sep.index <- (!grepl('/|\\|', called.alleles) & !is.na(called.alleles) & called.alleles != '.');
+            called.alleles[no.sep.index] <- paste0(called.alleles[no.sep.index], '/-');
             split.alleles <- data.table::tstrsplit(called.alleles, split = c('/|\\|'), keep = c(1,2)); # '|' are special chars in regular expressions
             }
     names(split.alleles) <- c('called.allele.a', 'called.allele.b');

R/run-pgs-statistics.R

@@ -62,7 +62,7 @@ classify.variable.type <- function(data, continuous.threshold = 4) {
     continuous.vars.index <- sapply(
         X = data,
         FUN = function(x) {
-            'numeric' == class(x) & continuous.threshold < length(unique(na.omit(x)));
+            ('numeric' == class(x) | 'integer' == class(x)) & continuous.threshold < length(unique(na.omit(x)));
             }
         );
     binary.vars.index <- sapply(

man/convert.alleles.to.pgs.dosage.Rd

@@ -34,13 +34,6 @@ test_that(
             );
 
         # check for correct called.alleles format
-        expect_error(
-            convert.alleles.to.pgs.dosage(
-                called.alleles = c('A/A', 'A'),
-                risk.alleles = c('A', 'T')
-                ),
-            'unrecognized called.alleles format, must be capitalized letters, "." or "\\*" separated by a slash or pipe.'
-            );
         expect_error(
             convert.alleles.to.pgs.dosage(
                 called.alleles = c('A/A', 'A,'),
@@ -101,8 +94,8 @@ test_that(
         # check that correct input is accepted
         expect_silent(
             convert.alleles.to.pgs.dosage(
-                called.alleles = c('A/A', 'A|T', 'TA/T', 'A/ATTTT', './.', '.', '*/T', 'T/*', '*/*'),
-                risk.alleles = c('A', 'T', 'A', 'T', 'A', 'T', 'A', 'T', 'A')
+                called.alleles = c('A/A', 'A|T', 'TA/T', 'A/ATTTT', './.', '.', '*/T', 'T/*', '*/*', 'A', 'T'),
+                risk.alleles = c('A', 'T', 'A', 'T', 'A', 'T', 'A', 'T', 'A', 'T', 'A')
                 )
             );
     }
@@ -132,6 +125,18 @@ test_that(
         }
     );
 
+test_that(
+    'convert.alleles.to.pgs.dosage calculates dosage correctly from hemizygous genotypes', {
+        expect_equal(
+            convert.alleles.to.pgs.dosage(
+                called.alleles = c('A', 'T'),
+                risk.alleles = c('A', 'A')
+                ),
+            c(1, 0)
+            );
+        }
+    );
+
 test_that(
     'convert.alleles.to.pgs.dosage calculates dosage correctly from short indel alleles', {
         expect_equal(

tests/testthat/test-dosage-calculator.R

@@ -206,6 +206,12 @@ test_that(
             assess.pgs.vcf.allele.match('ATCG', 'A', 'G', 'A'),
             'Mismatch detected in INDEL VCF allele. Skipping strand flip assessment.'
             );
+
+        # indel effect switch case
+        expect_silent(
+            assess.pgs.vcf.allele.match('A', 'ATCG', 'ATCG', 'A')
+            );
+
         # indel but everything matches
         expect_silent(
             assess.pgs.vcf.allele.match('ATCG', 'A', 'ATCG', 'A')
@@ -214,26 +220,26 @@ test_that(
         # check indel handling when return.indels.as.missing == FALSE
         # VCF ALT allele is an INDEL
         test.leave.indels.vcf.alt <- assess.pgs.vcf.allele.match(
-            c('A', 'A'),
-            c('G', 'ATCG'),
-            c('A', 'A'),
-            c('G', 'G'),
+            c('A', 'A', 'A'),
+            c('G', 'ATCG', 'ATCG'),
+            c('A', 'A', 'ATCG'),
+            c('G', 'G', 'A'),
             return.indels.as.missing = FALSE
             );
 
         expect_equal(
             test.leave.indels.vcf.alt$match.status,
-            c('default_match', 'indel_mismatch')
+            c('default_match', 'indel_mismatch', 'effect_switch')
             );
 
         expect_equal(
             test.leave.indels.vcf.alt$new.pgs.effect.allele,
-            c('G', 'G')
+            c('G', 'G', 'A')
             );
 
         expect_equal(
             test.leave.indels.vcf.alt$new.pgs.other.allele,
-            c('A', 'A')
+            c('A', 'A', 'ATCG')
             );
 
         # VCF REF allele is an INDEL