Polygenic Risk Score Calculation

Polygenic means multiple genes are involved in the emergence of the given phenotype. Gene is a fuzzy concept with evolving definitions and boundaries that are used as a catch-all phrase. (Cipriano and Ballarino 2018) In our case what we are interested in are loci in the genome that carries the variation that effect the same phenotype. This variation can be in a gene that shows its effect through its impact on the protein like a change from one amino acid to another. However, regions outside of the genes also play a role in biological processes and the emergence of phenome which we will see in the association studies. Since association studies show correlations and not causation a variant like rs1421085 found on FTO affects ARID5B millions of bases away. (Claussnitzer et al. 2015)

Oxford English Dictionary defines risk as:

the possibility of something bad happening at some time in the future.

In statistics, there exist two types of risk. (“Understanding Statistics: Risk BMJ Best Practice” n.d.)

• Absolute Risk (AR) ratio of the events to the total number in a given group
• Relative Risk (RR) ratio of the event happening in the case group to an event happening in the control group.

Here the term risk refers to a relative risk since we are comparing the case group with the control group. (“Polygenic Risk Scores” n.d.) RR should be evaluated with absolute risk on an individual basis. (Andrade 2015)

We can look at two different equations that calculate PRS. They both calculate by summing the effect size (denoted as in Collister et al. and in PLINK formula) and dosage (denoted as in Collister et al. and in PLINK formula) of the alleles. Differently, the PLINK formula has in the denominator. Where is the ploidy and is the number of non-missing Single Nucleotide Polymorphisms (SNPs). If a sample missing SNPs their score would be lower than samples with full genotypes. Dividing by non-missing number of SNPs and ploidy makes the score comparable between individuals with missing genotypes. (Collister, Liu, and Clifton 2022) With enough samples, these scores would approach a normal distribution.

1

(Collister, Liu, and Clifton 2022)

2

(Choi, Mak, and O’Reilly 2020)

GWAS performs statistical tests on thousands of genomic loci in thousands of genomes to see if any of those regions are associated with the given trait. (Uffelmann et al. 2021) As a result, we get p-values for how significant that association is and the effect size of the variation.

Let’s say we have found 10 SNPs significant in a study with 50000 cases and 50000 controls. First, let’s go over the SNPs.

snp_number <- 10

gwas <- data.frame(
  SNP_id = paste("SNP", seq(snp_number), sep="")
)
gwas[c(
  "non_effect_allele",
  "effect_allele"
)] <- t(replicate(snp_number,
  sample(c("A", "T", "C", "G"), 2)
  ))

print(rbind(head(gwas, 3), tail(gwas, 3)))

I selected some genotypes randomly. However, we don’t need to know the non-effect allele and effect allele nucleotides, that’s why they will be simply named A or 0 for the non-effect allele and B or 1 for the effect allele.

gwas$non_effect_allele = "A"
gwas$effect_allele = "B"
print(rbind(head(gwas, 3), tail(gwas, 3)))

Another thing we need from GWAS summary statistics is the effect sizes of these SNPs. We are going to use the log of OR as the effect size.

gwas_control_number = 50000
gwas_case_number = 50000
increment_percent = 0.02
gwas_sample_number = gwas_control_number + gwas_case_number

step_case = round(gwas_case_number * increment_percent)
gwas$obs_case = seq(gwas_case_number / 2 + step_case,
                    gwas_case_number / 2 + step_case * nrow(gwas),
                    step_case)

step_control = round(gwas_control_number * increment_percent)
gwas$obs_control = seq(gwas_control_number / 2 - step_control,
                       gwas_control_number / 2 - step_control * nrow(gwas),
                       -step_control)

gwas$log_OR <- signif(log(
    (gwas$obs_case / (gwas_case_number - gwas$obs_case)) /
    (gwas$obs_control / (gwas_control_number - gwas$obs_control))
), 2)

print(rbind(head(gwas, 3), tail(gwas, 3)))

We get effect sizes for the first genotype and for the last genotype. A BB genotype would have score of for the first genotype and .

We are screening our population concerning a phenotype. We aim to determine risk scores using the GWAS summary statistics and see if we can make a correlation between PRS and phenotype. We have genotyped the same 10 SNPs in 2000 people for our study.

set.seed(42)

sample_number <- 2000
snp_number <- 10
samples <- data.frame(
  sample_id = paste("Sample", seq(sample_number), sep=""),
  status = c(
    rep("Control", sample_number / 2),
    rep("Case", sample_number / 2))
)

get_genotype <- function(prob=c(.5, .5)) {
  paste(sample(c("A", "B"), 2, replace=T, prob = prob), collapse="")
}

control_prob <- c(0.55, 0.45)
case_prob <- c(0.45, 0.55)
samples[gwas$SNP_id] <- rbind(
  matrix(
    replicate(
      snp_number * sample_number / 2,
      get_genotype(control_prob)),
    nrow=sample_number / 2),
  matrix(
    replicate(
      snp_number * sample_number / 2,
      get_genotype(case_prob)),
    nrow=sample_number / 2)
)
print(rbind(head(samples, 3), tail(samples, 3)))

We get random genotypes for controls for alleles A, and B with probabilities 0.55, 0.45 ; cases for alleles A, and B with probabilities 0.45, 0.55 This way cases are more prone to having the effect allele.

We will be using the simpler version of the equations. Mainly because we don’t have any missing genotypes and secondly because it’s simpler.

We implement the formula by using an apply function to multiply dosage and the effect size then doing a row sum would give us the score.

dosage <- list("AA" = 0, "AB" = 1, "BA" = 1, "BB" = 2)

samples$PRS <- colSums(as.data.frame(
  apply(
    samples[gwas$SNP_id],
    1,
    function(x) as.numeric(dosage[x]) * gwas$log_OR)))

print(rbind(head(samples, 3), tail(samples, 3)))

What we will be seeing around in papers and around the web is a plot of normal distributions comparing cases to controls. To create the plot we will need ggplot2.

library(ggplot2)

We shift the scores by subtracting the mean of the controls to make the controls center around 0.

samples$status <- as.factor(with(samples, reorder(status, PRS)))
samples$PRS <- samples$PRS - mean(samples[which(samples$status=="Control"), "PRS"])
mean_control <- mean(samples[which(samples$status=="Control"), "PRS"])
sd_control <- sd(samples[which(samples$status=="Control"), "PRS"])

mean_case <- mean(samples[which(samples$status=="Case"), "PRS"])
sd_case <- sd(samples[which(samples$status=="Case"), "PRS"])

ggplot(samples, aes(x=PRS, color=status, fill=status)) +
  geom_histogram(aes(y=..density..), position="dodge", alpha=.5, bins=30) +
  stat_function(fun = dnorm,
                color="darkred",
                args = list(mean = mean_control, sd = sd_control)) +
  stat_function(fun = dnorm,
                color="darkblue",
                args = list(mean = mean_case, sd = sd_case )) +
  ggtitle(sprintf(
    "PRS distributions Controls μ %.2f σ %.2f Cases μ %.2f σ %.2f",
    mean_control, sd_control, mean_case, sd_case))

When you visit www.pgscatalog.org or image search “polygenic risk score” you would be seeing the plot of these two peaks. Cases have more effect alleles with higher totals on average this is why their curve is slightly shifted to the right. The case curve also has a larger standard deviation because effect alleles have increasing effect sizes.

Logistic regression is a model used in statistics that fit an S curve to data based on probabilities. Details of the subject are out of the scope of this post. I would suggest StatQuest’s logistic regression playlist to anyone interested.

Logistic regressions are used for prediction. However, one thing I would note in context to PRS calculations is that, unlike this common use, it’s not used for prediction in genetic studies. Because the scores and the risk calculated here are relative. The effect becomes tangible only in the tails of the curves but as the number of cases and controls drops the significance of the effect also drops.

ggplot(samples, aes(x=PRS, y=as.numeric(samples$status) -1, color=status)) +
  geom_point(shape = "|", position = position_jitter(w = 0.5, h = 0)) +
  geom_hline(yintercept = c(0,1), linetype = "dashed", color = "grey") +
  scale_y_discrete(name ="Status", labels=c("Control","Case"), limits=c(0,1))

model <- glm(status~PRS, family="binomial", data=samples)
summary(model)

Call:
glm(formula = status ~ PRS, family = "binomial", data = samples)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.30654    0.07353  -4.169 3.06e-05 ***
PRS          0.34170    0.03149  10.852  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1386.3  on 999  degrees of freedom
Residual deviance: 1242.5  on 998  degrees of freedom
AIC: 1246.5

Number of Fisher Scoring iterations: 4

In the model summary, we can find coefficients

for the control group and

for the case group. We can write our equation as

• .

Inserting the sample’s score into this equation would give us the log odds of having the phenotype. One unit increase in the would increase the log odds by units.

We can translate the log of odds to probabilities. Because it’s log we can turn them back to odds with exponential function. Odds are the event happening (let’s say ) ratio to the event not happening (let’s say ) which would be . gives us the log of the odds; gives us the odds. We can write the odds as ; this would be which would simplify to . It would give us the probability of the event happening over the total events.

3

(“An Introduction to Statistical Learning” n.d.)

We can plot our model’s S-shaped curve that shows the probabilities.

ggplot(samples, aes(x=PRS, y=as.numeric(samples$status) -1)) +
  geom_point(shape = "|", position = position_jitter(w = 0.5, h = 0), aes(color=status)) +
  geom_smooth(method = "glm",
    method.args = list(family = "binomial"),
    se = FALSE) +
  geom_hline(yintercept = c(0,1), linetype = "dashed", color = "grey") +
  scale_y_discrete(name ="Status", labels=c("Control","Case"), limits=c(0,1))

One of the measures of how good or model is . in linear regression is calculated using the sum of squares of residuals. This calculation does not apply to the logistic regression because we have a probability curve. There are several approaches to calculating and the most used one I have seen is the Nagelkerke method (also the most optimistic). The closer it gets to 1 better the fit.

library(rcompanion)

nagelkerke(model)

$Models

Model: "glm, status ~ PRS, binomial, samples"
Null:  "glm, status ~ 1, binomial, samples"

$Pseudo.R.squared.for.model.vs.null
                             Pseudo.R.squared
McFadden                             0.103711
Cox and Snell (ML)                   0.133917
Nagelkerke (Cragg and Uhler)         0.178556

$Likelihood.ratio.test
 Df.diff LogLik.diff  Chisq    p.value
      -1     -71.887 143.77 3.9798e-33

$Number.of.observations

Model: 1000
Null:  1000

$Messages
[1] "Note: For models fit with REML, these statistics are based on refitting with ML"

$Warnings
[1] "None"

This is the table-2 from the paper (Leu et al. 2019). Here it looks at the most extreme tails of the curve to see if the effect is larger.

We get the samples bigger and smaller than our percent thresholds with the quantile function. We then pull the values used in the table from the summary function. Differently, we use the confusionMatrix from the caret library to calculate the sensitivity and specificity.

What we are looking at in this table is the increase in the odds in the more extreme percentages. One odd thing I notice is that in my simulations sensitivity goes up and specificity goes down as we get to the edge of the tails where the opposite is seen in the paper.

library(caret)

my_sum_rows = list()
percents = c(20, 10, 5, 1, 0.5)
for (percent in percents) {
    samples$percentile <- ifelse(
      samples$PRS >= quantile(samples$PRS,prob=1-percent/100),
      "Above", "Below")
    samples$percentile <- as.factor(with(samples, reorder(percentile, PRS)))

    model <- glm(status~percentile, family="binomial", data=samples)

    threshold <- 0.5
    cm <- confusionMatrix(
      as.factor(ifelse(predict(model,type="response")> threshold, 1, 0)),
      as.factor(model$y)
    )

    above_control <- nrow(
      samples[which(samples$status=="Control" & samples$percentile=="Above"),]
    )
    below_control <- nrow(
      samples[which(samples$status=="Control" & samples$percentile=="Below"),]
    )
    above_case <- nrow(
      samples[which(samples$status=="Case" & samples$percentile=="Above"),]
    )
    below_case <- nrow(
      samples[which(samples$status=="Case" & samples$percentile=="Below"),]
    )

    my_sum_row <- data.frame(
    "Reference Group" = sprintf("Remaining %s%% ", 100-percent),
    "Cases/Controls Above PRS%" = sprintf("%s/%s", above_case, above_control),
    "Cases/Controls Below PRS%" = sprintf("%s/%s", below_case, below_control),
    "Odds Ratio" = round(exp(coef(model)[2]), 2),
    "95% CI" = do.call(sprintf, c(
                                  "%.2f - %.2f",
                                  as.list(exp(confint(model, level=.90)[2,])))),
    "P value" = sprintf("%.2e",
                        summary(model)$coefficients['percentileAbove', 'Pr(>|z|)']
                        ),
    "Sensitivity / Specificity" = sprintf(
        "%.3f / %.3f",
        as.list(cm$byClass)$Sensitivity,
        as.list(cm$byClass)$Specificity
      )
    )
    my_sum_rows[[sprintf("Top %s%% of distribution", percent)]] <- my_sum_row
}
my_sum_table <- do.call(rbind, my_sum_rows)
colnames(my_sum_table) <- c(
  "Reference group",
  "Cases/Controls above PRS%",
  "Cases/Controls below PRS%",
  "Odds ratio",
  "95% CI",
  "P value",
  "Sensitivity / Specificity"
)
my_sum_table

This graph is from the paper (Kloeve-Mogensen et al. 2021). Here they compare scores from each decile. The fifth decile is used as the reference that’s why we are seeing a stair pattern where the decile before the fifth gets closer and the deciles after go up incrementally. I have found this answer which sets the fifth decile as a reference differently. (socialscientist 2022) I just set the levels so the fifth is the first one which the generalized linear models (glm) uses as the reference. Resulting coefficients are considered as odds ratio.(Noah 2018)

samples$deciles <- factor(.bincode(samples$PRS,
        breaks = quantile(samples$PRS, seq(0, 1, 0.1)),
        include.lowest = TRUE
      ), levels = c(5, 1, 2, 3, 4, 6, 7, 8, 9, 10))

model <- glm(status~deciles, family="binomial", data=samples)
summary(model)

Call:
glm(formula = status ~ deciles, family = "binomial", data = samples)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.19867    0.14141  -1.405  0.16006
deciles1    -1.32376    0.23203  -5.705 1.16e-08 ***
deciles2    -0.88014    0.21511  -4.092 4.28e-05 ***
deciles3    -0.36885    0.20435  -1.805  0.07108 .
deciles4    -0.03352    0.20092  -0.167  0.86749
deciles6     0.50095    0.20114   2.491  0.01276 *
deciles7     0.58733    0.20242   2.902  0.00371 **
deciles8     1.02227    0.20877   4.897 9.75e-07 ***
deciles9     1.07696    0.20974   5.135 2.82e-07 ***
deciles10    1.26406    0.21532   5.871 4.34e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2772.6  on 1999  degrees of freedom
Residual deviance: 2484.1  on 1990  degrees of freedom
AIC: 2504.1

Number of Fisher Scoring iterations: 4

my_sum_data <- as.data.frame(cbind(
  exp(confint(model, level=.90)[2:10,]),
  "Odds ratio"=round(exp(coef(model)[2:10]), 2),
  "P value"=summary(model)$coefficients[2:10, 'Pr(>|z|)']
))
colnames(my_sum_data) <- c(
  "lower",
  "upper",
  "odds_ratio",
  "p_val"
)
my_sum_data <- rbind(my_sum_data, c(1,1,1,1))
my_sum_data$decile <- factor(c(1,2,3,4,6,7,8,9,10,5))

ggplot(my_sum_data, aes(decile, odds_ratio)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width=.2) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey") +
  geom_point(shape=22, size=6, color="black", fill="orange", stroke=.8)