ApplyPolygenicScore (version 4.0.1)

apply.polygenic.score: Apply polygenic score to VCF data

Description

Apply a polygenic score model to VCF data.

Usage

apply.polygenic.score(
  vcf.data,
  vcf.long.format = FALSE,
  pgs.weight.data,
  phenotype.data = NULL,
  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,
  missing.genotype.method = "mean.dosage",
  use.external.effect.allele.frequency = FALSE,
  n.percentiles = NULL,
  analysis.source.pgs = NULL,
  validate.inputs.only = FALSE
)

Value

A list containing per-sample PGS output and per-phenotype regression output if phenotype analysis columns are provided.

Output Structure

The outputed list contains the following elements:

• pgs.output: A data.frame containing the PGS per sample and optional phenotype data.

• regression.output: A data.frame containing the results of the regression analysis if phenotype.analysis.columns are provided, otherwise NULL.

pgs.output columns:

• Indiv: A character string indicating the sample ID.

• PGS: A numeric vector indicating the PGS per sample. (only if missing.genotype.method includes "none")

• PGS.with.normalized.missing: A numeric vector indicating the PGS per sample with missing genotypes normalized. (only if missing.genotype.method includes "normalize")

• PGS.with.replaced.missing: A numeric vector indicating the PGS per sample with missing genotypes replaced by mean dosage. (only if missing.genotype.method includes "mean.dosage")

• percentile: A numeric vector indicating the percentile rank of the PGS.

• decile: A numeric vector indicating the decile rank of the PGS.

• quartile: A numeric vector indicating the quartile rank of the PGS.

• percentile.X: A numeric vector indicating the user-specified percentile rank of the PGS where "X" is substituted by n.percentiles. (only if n.percentiles is specified)

• n.missing.genotypes: A numeric vector indicating the number of missing genotypes per sample.

• percent.missing.genotypes: A numeric vector indicating the percentage of missing genotypes per sample.

• All columns in phenotype.data if provided.

regression.output columns:

• phenotype: A character vector of phenotype names.

• model: A character vector indicating the regression model used. One of "logistic.regression" or "linear.regression".

• beta: A numeric vector indicating the beta coefficient of the regression analysis.

• se: A numeric vector indicating the standard error of the beta coefficient.

• p.value: A numeric vector indicating the p-value of the beta coefficient.

• r.squared: A numeric vector indicating the r-squared value of linear regression analysis. NA for logistic regression.

• AUC: A numeric vector indicating the area under the curve of logistic regression analysis. NA for linear regression.

PGS Calculation

PGS for each individual i is calculated as the sum of the product of the dosage and beta coefficient for each variant in the PGS: $$PGS_i = \sum_{m=1}^{M} \left( \beta_m \times dosage_{im} \right)$$

Where m is a PGS component variant out of a total M variants.

Missing Genotype Handling

VCF genotype data are matched to PGS data by chromosome and position. If a SNP cannot be matched by genomic coordinate, an attempt is made to match by rsID (if available). If a SNP from the PGS weight data is not found in the VCF data after these two matching attempts, it is considered a cohort-wide missing variant.

Missing genotypes (in individual samples) among successfully matched variants are handled by three methods:

none: Missing genotype dosages are excluded from the PGS calculation. This is equivalent to assuming that all missing genotypes are homozygous for the non-effect allele, resulting in a dosage of 0.

normalize: Missing genotypes are excluded from score calculation but the final score is normalized by the number of non-missing alleles. The calculation assumes a diploid genome: $$PGS_i = \dfrac{\sum \left( \beta_m \times dosage_{im} \right)}{P_i * M_{non-missing}}$$

Where P is the ploidy and has the value 2 and \(M_{non-missing}\) is the number of non-missing genotypes.

mean.dosage: Missing genotype dosages are replaced by the mean population dosage of the variant which is calculated as the product of the effect allele frequency EAF and the ploidy of a diploid genome: $$\overline{dosage_{k}} = EAF_k * P$$

where k is a PGS component variant that is missing in between 1 and n-1 individuals in the cohort and P = ploidy = 2 This dosage calculation holds under assumptions of Hardy-Weinberg equilibrium. By default, the effect allele frequency is calculated from the provided VCF data. For variants that are missing in all individuals (cohort-wide), dosage is assumed to be zero (homozygous non-reference) for all individuals. An external allele frequency can be provided in the pgs.weight.data as a column named allelefrequency_effect and by setting use.external.effect.allele.frequency to TRUE.

Multiallelic Site Handling

If a PGS weight file provides weights for multiple effect alleles, the appropriate dosage is calculated for the alleles that each individual carries. It is assumed that multiallelic variants are encoded in the same row in the VCF data. This is known as "merged" format. Split multiallelic sites are not accepted. VCF data can be formatted to merged format using external tools for VCF file manipulation.

Allele Mismatch Handling

Variants from the PGS weight data are merged with records in the VCF data by genetic coordinate. After the merge is complete, there may be cases where the VCF reference (REF) and alternative (ALT) alleles do not match their conventional counterparts in the PGS weight data (other allele and effect allele, respectively). This is usually caused by a strand flip: the variant in question was called against opposite DNA reference strands in the PGS training data and the VCF data. Strand flips can be detected and corrected by flipping the affected allele to its reverse complement. apply.polygenic.score uses assess.pgs.vcf.allele.match to assess allele concordance, and is controlled through the following arguments:

• correct.strand.flips: When TRUE, detected strand flips are corrected by flipping the affected value in the effect_allele column prior to dosage calling.

• remove.ambiguous.allele.matches: Corresponds to the return.ambiguous.as.missing argument in assess.pgs.vcf.allele.match. When 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 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.

• max.strand.flips: This argument only applies when 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 max.strand.flips.

• remove.mismatched.indels: Corresponds to the return.indels.as.missing argument in assess.pgs.vcf.allele.match. When TRUE, INDEL allele mismatches (which cannot be assessed for strand flips) are removed by marking the affected value in the 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.

Note that an allele match assessment requires the presence of both the other_allele and effect_allele in the PGS weight data. The other_allele column is not required by the PGS Catalog, and so is not always available.

Arguments

vcf.data

VCF genotype data as formatted by import.vcf(). Two formats are accepted: wide format (a list of elements named genotyped.alleles and vcf.fixed.fields) or long format(a data frame). See vcf.import for more details.

vcf.long.format

A logical indicating whether vcf.data is provided in long format. Default is FALSE.

pgs.weight.data

A data.frame containing PGS weight data as formatted by import.pgs.weight.file().

phenotype.data

A data.frame containing phenotype data. Must have an Indiv column matching vcf.data. Default is NULL.

phenotype.analysis.columns

A character vector of phenotype columns from phenotype.data to analyze in a regression analsyis. Default is NULL. Phenotype variables are automatically classified as continuous, binary, or neither based on data type and number of unique values. The calculated PGS is associated with each phenotype variable using linear or logistic regression for continuous or binary phenotypes, respectively. See run.pgs.regression for more details. If no phenotype.analysis.columns are provided, no regression analysis is performed.

correct.strand.flips

A logical indicating whether to check PGS weight data/VCF genotype data matches for strand flips and correct them. Default is TRUE. The PGS catalog standard column other_allele in pgs.weight.data is required for this check.

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 FALSE. The PGS catalog standard column other_allele in pgs.weight.data is required for this check.

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 0 which means that all ambiguous variants are removed regardless of the status of any other variant.

remove.mismatched.indels

A logical indicating whether to remove indel variants that are mismatched between PGS weight data and VCF genotype data. Default is FALSE. The PGS catalog standard column other_allele in pgs.weight.data is required for this check.

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. Files are tab-separate .txt files. Default is NULL in which case no files are written.

file.prefix

A character string to prepend to the output file names. Default is NULL.

missing.genotype.method

A character string indicating the method to handle missing genotypes. Options are "mean.dosage", "normalize", or "none". Default is "mean.dosage".

use.external.effect.allele.frequency

A logical indicating whether to use an external effect allele frequency for calculating mean dosage when handling missing genotypes. Default is FALSE. Provide allele frequency as a column is pgs.weight.data named allelefrequency_effect.

n.percentiles

An integer indicating the number of percentiles to calculate for the PGS. Default is NULL.

analysis.source.pgs

A character string indicating the source PGS for percentile calculation and regression analyses. Options are "mean.dosage", "normalize", or "none". When not specified, defaults to missing.genotype.method choice and if more than one PGS missing genotype method is chosen, calculation defaults to the first selection.

validate.inputs.only

A logical indicating whether to only perform input data validation checks without running PGS application. If no errors are triggered, a message is printed and TRUE is returned. Default is FALSE.

Examples

# Example VCF
vcf.path <- system.file(
    'extdata',
    'HG001_GIAB.vcf.gz',
    package = 'ApplyPolygenicScore',
    mustWork = TRUE
    );
vcf.import <- import.vcf(vcf.path, long.format = TRUE);

# Example pgs weight file
pgs.weight.path <- system.file(
    'extdata',
    'PGS000662_hmPOS_GRCh38.txt.gz',
    package = 'ApplyPolygenicScore',
    mustWork = TRUE
    );
pgs.import <- import.pgs.weight.file(pgs.weight.path);

pgs.data <- apply.polygenic.score(
    vcf.data = vcf.import$split.wide.vcf.matrices,
    pgs.weight.data = pgs.import$pgs.weight.data,
    missing.genotype.method = 'none'
    );

# Use long format
pgs.data <- apply.polygenic.score(
    vcf.data = vcf.import$combined.long.vcf.df$dat,
    vcf.long.format = TRUE,
    pgs.weight.data = pgs.import$pgs.weight.data,
    missing.genotype.method = 'none'
    );

# Specify different methods for handling missing genotypes
pgs.import$pgs.weight.data$allelefrequency_effect <- rep(0.5, nrow(pgs.import$pgs.weight.data));
pgs.data <- apply.polygenic.score(
    vcf.data = vcf.import$split.wide.vcf.matrices,
    pgs.weight.data = pgs.import$pgs.weight.data,
    missing.genotype.method = c('none', 'mean.dosage', 'normalize'),
    use.external.effect.allele.frequency = TRUE
    );

# Specify allele mismatch handling
pgs.data <- apply.polygenic.score(
   vcf.data = vcf.import$split.wide.vcf.matrices,
   pgs.weight.data = pgs.import$pgs.weight.data,
   correct.strand.flips = TRUE,
   remove.ambiguous.allele.matches = TRUE,
   remove.mismatched.indels = FALSE
   );

# Provide phenotype data for basic correlation analysis
n.samples <- length(colnames(vcf.import$split.wide.vcf.matrices$genotyped.alleles))
phenotype.data <- data.frame(
    Indiv = colnames(vcf.import$split.wide.vcf.matrices$genotyped.alleles),
    continuous.phenotype = rnorm(n.samples),
    binary.phenotype = sample(
        c('a', 'b'),
        n.samples,
        replace = TRUE
        )
    );

pgs.data <- apply.polygenic.score(
    vcf.data = vcf.import$split.wide.vcf.matrices,
    pgs.weight.data = pgs.import$pgs.weight.data,
    phenotype.data = phenotype.data
    );

# Only run validation checks on input data and report back
apply.polygenic.score(
    vcf.data = vcf.import$split.wide.vcf.matrices,
    pgs.weight.data = pgs.import$pgs.weight.data,
    validate.inputs.only = TRUE
    );