apply.polygenic.score | R Documentation |
Apply a polygenic score to VCF data.
apply.polygenic.score(
vcf.data,
pgs.weight.data,
phenotype.data = NULL,
phenotype.analysis.columns = NULL,
correct.strand.flips = TRUE,
remove.ambiguous.allele.matches = FALSE,
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
)
vcf.data |
A data.frame containing VCF genotype data as formatted by |
pgs.weight.data |
A data.frame containing PGS weight data as formatted by |
phenotype.data |
A data.frame containing phenotype data. Must have an Indiv column matching vcf.data. Default is |
phenotype.analysis.columns |
A character vector of phenotype columns from phenotype.data to analyze in a regression analsyis. Default is |
correct.strand.flips |
A logical indicating whether to check PGS weight data/VCF genotype data matches for strand flips and correct them. Default is |
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 |
remove.mismatched.indels |
A logical indicating whether to remove indel variants that are mismatched between PGS weight data and VCF genotype data. Default is |
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 |
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 |
n.percentiles |
An integer indicating the number of percentiles to calculate for the PGS. Default is |
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 |
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 |
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.
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.
# Example VCF
vcf.path <- system.file(
'extdata',
'HG001_GIAB.vcf.gz',
package = 'ApplyPolygenicScore',
mustWork = TRUE
);
vcf.import <- import.vcf(vcf.path);
# 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$dat,
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$dat,
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$dat,
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
phenotype.data <- data.frame(
Indiv = unique(vcf.import$dat$Indiv),
continuous.phenotype = rnorm(length(unique(vcf.import$dat$Indiv))),
binary.phenotype = sample(
c('a', 'b'),
length(unique(vcf.import$dat$Indiv)),
replace = TRUE
)
);
pgs.data <- apply.polygenic.score(
vcf.data = vcf.import$dat,
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$dat,
pgs.weight.data = pgs.import$pgs.weight.data,
validate.inputs.only = TRUE
);
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.