PGS Lesson

Datasets

Coronary Artery PGS

First, we need to download the coronary artery polygenic score data from https://www.pgscatalog.org/score/PGS000010/. This particular PGS is composed of 27 SNPs identified from a GWAS of 86,995 European individuals.

# download the coronary artery PGS data
URL = "https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/PGS000010/ScoringFiles/PGS000010.txt.gz"
download.file(URL, destfile="PGS000010.txt.gz", quiet=T)
pgs = fread("PGS000010.txt.gz") %>% mutate(effect_weight=round(effect_weight, 4))
datatable(pgs, rownames=F) 

23 and Me Examples

Then, we will load in two sample genotypes from 23 and Me and merge them with the PGS data. Additionally, we need to calculate the allele counts (AC) by adding up the number of PGS effect alleles in each genotype. You can see below that the columns geno1, geno2, AC1, and AC2 have been added to the PGS data. You can also see that 23 and Me information is missing for two of the SNPs (rs3825807 and rs9411489).

# read in the 23 and Me example data
ex = read.csv("23andMeExamples.csv")

# merge the pgs and example data
pgs2 = merge(pgs, ex, by="rsID")

# count the number of effect alleles for genotype 1
geno1 = sapply(strsplit(pgs2$geno1, "/"), function(x) x[1:2])
pgs2$AC1 = (geno1[1,]==pgs2$effect_allele) + (geno1[2,]==pgs2$effect_allele)

# count the number of effect alleles for genotype 2
geno2 = sapply(strsplit(pgs2$geno2, "/"), function(x) x[1:2])
pgs2$AC2 = (geno2[1,]==pgs2$effect_allele) + (geno2[2,]==pgs2$effect_allele)
datatable(pgs2, rownames=F) 

1000 Genomes

Now, we need to load in the individual-level 1000 Genomes data. This data will be used to calculate distributions of polygenic scores that we will compare the 23 and Me examples to. A subset of the data is shown below, where the first five columns provide information related to each SNP and the next four columns beginning with HG correspond to individual genotypes. The entire dataset contains genotypes for 2,504 individuals. The ID column below corresponds to the rsID column from the PGS data above.

# read in the 1000 genomes data
data = read.table("1000G.phase3.PRS.vcf", sep="\t", header=T)
samples = colnames(data)[10:ncol(data)] # subset the sample names
datatable(data[,c(1:5, 10:13)], rownames=F)

Notice how the genotypes in the 1000 Genomes data are designated by 0’s and 1’s instead of by alleles (A, T, C, or G) as in the 23 and Me data. Here, a 0 represents the reference allele (REF) and a 1 represents the alternate allele (ALT). So, looking at the first SNP, the genotype 0|0 for the first individual (HG00096) corresponds to T/T, whereas the genotype 1|0 for the second individual (HG00097) corresponds to C/T. These genotypes can easily be converted to alternate allele counts by just adding up the 0’s and 1’s for each genotype as seen below.

The allele frequencies for the European (EUR_AF) and African (AFR_AF) populations were also added to the data. These frequencies describe how often the alternate allele is seen in the respective population. Notice how for some SNPs, the two frequencies are very similar, but for others, such as rs17465637, they can be very different. Also notice that the 1000 Genomes data only has 26 rows while the PGS data has 27. One SNP (rs9411489) is missing from the 1000 Genomes data.

# convert the genotypes into allele counts
data2 = geno(data)

# extract the European allele frequencies from the INFO column
data2$EUR_AF = sapply(strsplit(data2$INFO, ";"), function(x) x[9])
data2$EUR_AF = sapply(strsplit(data2$EUR_AF, "="), function(x) x[2])

# extract the African allele frequencies from the INFO column
data2$AFR_AF = sapply(strsplit(data2$INFO, ";"), function(x) x[8])
data2$AFR_AF = sapply(strsplit(data2$AFR_AF, "="), function(x) x[2])

data2 = data2 %>% relocate(EUR_AF:AFR_AF, .after=ALT)
datatable(data2[,c(1:7, 12:15)], rownames=F)

Data Cleaning

Merge Data

Next, the 1000 Genomes data needs to be merged with the PGS and 23 and Me data. Remember, the 1000 Genomes allele counts refer to the number of alternate alleles, so we want the effect_allele to correspond with the ALT allele. We also want the allelefrequency_effect and EUR_AF columns to be similar since the PGS data is derived from individuals with European ancestry. But, looking at the merged data below, that is not the case for some of the SNPs. However, we can switch the reference and alternate alleles for a particular SNP as long as we also update the allele counts (two minus the alternate allele count) and allele frequencies (one minus the alternate allele frequency).

You may also notice that the effect allele for one SNP (rs3825807) doesn’t match either the reference or alternate allele from 1000 Genomes. Additionally, the effect allele for another SNP (rs216172) does correspond with the alternate allele but the allelefrequency_effect and EUR_AF columns don’t match. For these two SNPs, we also need to do a strand flip. Remember that DNA is double-stranded and that A alleles on one strand pair with T alleles on the other and vice versa while C alleles on one strand pair with G alleles on the other and vice versa. We know that the 1000 Genomes and 23 and Me data refer to the positive strand, but we are unsure which strand the effect alleles from the PGS data are on. So, we will flip the strands for the two previously mentioned SNPs and assume that the effect alleles of the remaining SNPs are on the positive strand. However, this would have to be verified in real clinical or research applications.

# merge the 1000 genomes data with the pgs data
merged = left_join(pgs2, data2, by=c("rsID"="ID")) %>% 
  select(-QUAL, -FILTER, -INFO, -FORMAT) %>%
  relocate(CHROM:ALT, .after=rsID) %>%
  relocate(EUR_AF:AFR_AF, .after=allelefrequency_effect) %>%
  mutate(EUR_AF=as.numeric(EUR_AF), AFR_AF=as.numeric(AFR_AF))
datatable(merged[,1:20], rownames=F)

Looking at the updated data below, you can see that the effect_allele and ALT alleles now match and that the allelefrequency_effect and EUR_AF columns are more similar.

# update the 1000 genomes data based on the effect alleles in the pgs data
merged2 = update(merged)
datatable(merged2[,1:20], rownames=F)

Subset Ancestry

Now, we can look at the metadata for the 1000 Genomes data, which shows information about each individual. We can use this information to subset the merged data by ancestry population. In this exercise, we will calculate and compare the polygenic scores for the 403 Non-Finnish European individuals (EUR Superpopulation.code, excluding Finnish Population.name) with the scores for the 661 African individuals (AFR Superpopulation.code).

# read in the 1000 genomes metadata for the samples
info = read.table("1000G.phase3.metadata.tsv", sep="\t", header=T) %>% filter(Sample.name %in% samples) %>% arrange(Sample.name)
datatable(info[,-9], rownames=F)

# subset the sample names by population
NFE.samples = info %>% filter(Superpopulation.code=="EUR", Population.code!="FIN") %>% select(Sample.name)
AFR.samples = info %>% filter(Superpopulation.code=="AFR") %>% select(Sample.name)

# subset the 1000 genomes data by population
NFE = merged2 %>% select(rsID:AC2, NFE.samples$Sample.name) %>% select(-AFR_AF)
AFR = merged2 %>% select(rsID:AC2, AFR.samples$Sample.name) %>% select(-EUR_AF)

PGS

Calculations

To calculate the polygenic scores for the 1000 Genomes populations, we will multiply the effect_weight column by the allele counts in each HG column of the NFE and AFR data. We will also calculate the scores for the 23 and Me examples by multiplying the effect_weight column by the AC1 and AC2 columns. So, we will have 403 scores for the NFE population and 661 scores for the AFR population as well as two scores for the example data.

Note, prior to calculating the scores, we will remove the two SNPs missing from the 23 and Me data (rs3825807 and rs9411489) so all of the scores will be based on the same 25 SNPs (instead of the original 27 in the PGS data).

# remove NA values
NFE2 = NFE %>% filter(geno1!="No Data")
AFR2 = AFR %>% filter(geno1!="No Data")

# calculate the pgs for each population
NFE.pgs = NFE2$effect_weight %*% data.matrix(NFE2[,16:ncol(NFE2)])
AFR.pgs = AFR2$effect_weight %*% data.matrix(AFR2[,16:ncol(AFR2)])

# calculate the pgs for the examples
ex.pgs1 = NFE2$effect_weight %*% NFE2$AC1
ex.pgs2 = NFE2$effect_weight %*% NFE2$AC2

Plots

We can then plot the distributions of the scores from the two populations as well as the two 23 and Me example scores.

# re-format the data for plotting
pgs.data = data.frame(Population=c(rep("NFE", length(NFE.pgs)), rep("AFR", length(AFR.pgs))),
                      Score=c(NFE.pgs, AFR.pgs))
ex.data = data.frame(Example=c("1", "2"), Score=c(ex.pgs1, ex.pgs2))

# plot the pgs distributions
ggplot() +
  geom_density(data=pgs.data, aes(x=Score, fill=Population), alpha=0.75) +
  geom_vline(data=ex.data, aes(xintercept=Score, linetype=Example)) +
  scale_linetype_manual(values=c("1"="dashed", "2"="dotted")) +
  theme_bw() + ylab("Density")

The plot shows that the distributions look very different for the two populations, with the average score for the NFE population 3.09 and the average score for the AFR population 2.41. The plot also shows that the score for Example 2 is about one more than the score for Example 1. Comparing the example scores to the population distributions, Example 1 would have a low score compared to the NFE population, but a pretty high score compared to the AFR population. That is why it is SO important to account for ancestry when interpreting polygenic scores. Since the PGS data is derived from individuals with European ancestry, the example scores should only be compared to a European population in practice.

Percentages

We can also calculate the exact percentiles of the example scores relative to each population.

# estimate the distributions based on the data (empirical cdfs)
NFE.cdf = ecdf(NFE.pgs)
AFR.cdf = ecdf(AFR.pgs)

# calculate the percentiles for the first example
ex1.NFE = NFE.cdf(ex.pgs1)
ex1.AFR = AFR.cdf(ex.pgs1)

# calculate the percentiles for the second example
ex2.NFE = NFE.cdf(ex.pgs2)
ex2.AFR = AFR.cdf(ex.pgs2)

Based on the distributions shown in the plot above, the score for Example 1 is in the 29.03 percentile (low risk) for the NFE population, but is in the 91.53 percentile (high risk) for the AFR population. Additionally, the score for Example 2 is in the 95.53 percentile (high risk) for the NFE population, but is in the 99.85 percentile (very high risk) for the AFR population.

Note that the individuals from the 1000 Genomes data are considered “healthy”, so the distributions and percentages might be different if we instead compared the 23 and Me examples to a population of people with coronary artery disease or a family history of coronary artery disease.