Calculating PRS from a single WGS sample

[mfasold]

Hi! I would like to calculate PGS using DNA variants from a single WGS sample. I processed my data according to this discussion thread linked in the docs.

Is there a possibility to match homref variants when no ALT is given in the input VCF file?

Here is an example entry from my GVCF file that I obtained using the suggested processing steps with GATK:
chr3 12351626 . C . . . DP=40 GT:AD:DP:RGQ 0/0:40:40:99

My scorefile contains the following entry for this location:
3 12351626 C G 0.0453 additive False PGS000033_hmPOS_GRCh38 3

So the effect allele corresponds to the allele in my GVCF. Contrary to my expectation, the variant is not matched. It is however matched, if I manually put in G as ALT allele in the GVCF (replacing the "."). Since the genotype is 0/0, shouldn’t it be matched no matter what allele is in the ALT column?

[smlmbrt]

Hi @mfasold, you'll need to find a way to add in the reference base to the VCF in order to get the matching working. You may be able to do this using the reference dataset we provide (because all variants will inevitably have to intersect that if you would like to use --run_ancestry). It is important that the variants that are homozygous REF are tallied in the score, as they represent variation that needs to be compared with other individuals.

[mfasold]

Thank you for your reply! Please allow me to give a bit more detail on my toy example.

We are looking here at SNP rs1801282 which has the REF allele C, and (according to dbSNP) the ALT alleles G and T. In the scorefiles for my example trait PGS000033 (PGS000033_hmPOS_GRCh38.txt.gz and scorefiles.txt.gz) the effect_allele is C and the other_allele is G.

Since my original VCF did not contain hom-ref entries, I did variant calling with the GATK haplotype caller together with the positions from the scorefile to create a special GVCF that contains the explicit hom-ref call for all positions in the scorefile. For our SNP rs1801282, that gave the entry
chr3 12351626 . C . . . DP=40 GT:AD:DP:RGQ 0/0:40:40:99
(REF is C and ALT is empty “.”, with a 0/0 genotype.)

Note that I here focus on the “difficult” case, the homozygous reference SNPs, as other variants with a non-empty ALT allele will match, even with the original VCF. So I am already aware that the homref variants need to be included in the scoring and a special GVCF input file is needed for pgsc_calc. As mentioned, I did this and received an output where the two homref variants still did not match - see the *log.csv.gz file:

So when the reference base is already contained in the (G)VCF, do you mean I should add the other_allele to the (G)VCF? (And yes, I would like to use –run_ancestry later.) Because this is what I did by adding the G manually in the ALT column of the GVCF. Now the example SNP is not unmatched anymore (but also not matched like the others):

This seems not the way to go.

Where should I add the reference base? And what do refer to with the reference dataset?

[smlmbrt]

Hi @mfasold - adding the G (other_allele when the locus is HOMREF for the effect_allele) did fix it. However, the variant matched but ends up getting excluded from the scoring file because it is an ambiguous variant (C/G or A/T variants that can be strand ambiguous, these can be included with the --keep_ambiguous flag in matching options). Basically the . allele (e.g. C/.) will fail matching to a genotype or scoring file variant that is C/G as it is ambiguous - only scoring file variants can have a single effect allele, but the genotypes must always have a REF and ALT.

[Fiwx]

To my knowledge, it doesn't matter what the alternate allele is for PRS calculation as long as the sample doesn't have that allele. [Edit: this is wrong.] Here is something I have tried before. I would like to know if this is ok to do for pgsc_calc:

1. Extract the genotypes from the last column of the VCF line.
2. Check if the alternate allele is missing and the genotype is homozygous reference.
3. If the conditions in step 2 are met, replace the missing alternate allele with a complementary nucleotide based on the reference allele. (Or any other nucleotide.)
4. Returns the modified row list and replace the original genotypes with the new ones.

[PCyph]

@smlmbrt Okay, great. Would a solution be to simply add all possible alt alleles when the VCF position is homozygous reference?

For example (AA genotype):
Scorefile: A (effect allele) C (other allele)

Original VCF position: 0/0 A .
(Fails matching.)

Subpar solution: add an arbitrary alt allele: 0/0 A G
(Fails matching, as the scorefile is A C not A G.)

Current proposed solution: 0/0 A G,C,T
(Desired behavior: able to match as the scorefile, A C, is included in the possible alleles. Adds 2*effect allele beta to the score.)

[smlmbrt]

@PCyph that would probably allow the pipeline to run with our errors, but there is no guarantee it's correct (also I think you may have to recode the genotype fields?). The pgsc_calc pipeline is currently optimised to calculate PGS on VCFs that are output from common imputation pipelines. Our long-term goal is to support sequencing VCFs (specifically gVCFs where you know which blocks of genome were sequenced are are homref) - this will involve reworking a lot of the variant matching and scoring utilities though so it's a major piece of work. I think there may be other ways to get the scoring file positions and explicitly call those in your VCF with knowledge of which ALT alleles you're comparing to?

[PCyph]

Okay, thank you. That makes sense. I could get the scoring file and add the corresponding alt allele to the VCF for each homozygous reference position.

By not creating new multiallelic sites, and instead explicitly calling VCF alt alleles based on which ALT alleles are being compared:

• If the 1000G panel and the scorefile have two different alt alleles at a position, wouldn't this make it impossible for this variant to be matched correctly?

• A similar issue might arise if running two different scorefiles: if a position in the first scorefile is (effect=A, other=G), but a position in the second scorefile is (effect=A, other=T), then would this cause failure to match for at least one of these scores, regardless of the genotype?

[smlmbrt]

If the 1000G panel and the scorefile have two different alt alleles at a position, wouldn't this make it impossible for this variant to be matched correctly?

Currently, yes.

[...] if a position in the first scorefile is (effect=A, other=G), but a position in the second scorefile is (effect=A, other=T), then would this cause failure to match for at least one of these scores, regardless of the genotype?

Yes, but variants don't need to be in all scoring files for them to be calculated so it will still be included in the scoring file where it does match. You can find information about the variants that matched in the log file.

[mfasold]

I am dealing with a lot of different cases... Could you please help clear my mind and clarify if the following should work? We are looking at a hom-refs in the sample, genotyped for all locations occurring in the selected scorefiles. I want to modify those VCF entries so that they score correctly for all matching effect_allele / other-allele pairs in the selected score
files. For simplicity, we only consider SNPs.

When the location matches AND

• REF == effect_allele AND other_allele available -> put other_allele in ALT
• REF == effect_allele AND other_allele unavailable -> put any other base in in ALT
• other_allele available AND REF == other_allele -> put in effect_allele in ALT
• otherwise: accept that no matching is possible

[mfasold]

Do you have a suggestion on how to best obtain a Super-Scoring-File: it should contain, for my genome build hg38, all the variants (with chr, pos, effect allele, other allele) for all possible scores in the catalog? The idea, of course, would be to run the variant calling and VCF preparation once, so that afterward I will be able to run pgsc_calc on any list of pgs_id's that I'm interested in.

I saw that there is the download_scorefiles tool which could be used to loop over the pgs_scores_list.txt and then merge the results. But maybe there is an easier way.

[smlmbrt]

You can use the pgscatalog_utils python package (also on pypi) by first using the download_scorefiles command and merging those files with combine_scorefiles afterwards. This is what is happening within the pipeline currently.

Alternatively, you can do this within the pipeline by adding the --only_input flag. A combined long-form scoring file (with variants in your build) will be in the work directory where the COMBINE_SCOREFILES process has completed. This probably has higher overhead though.

[mfasold]

I went ahead with adding artificial ALT alleles to my VCF file which correspond to the effect alleles in the scoring file.

I find several cases where the variant does not match, opposed to my expectation. Here are two examples - first the VCF entry and then the output from the match log:

chr20   32759707        .       A       C       .       .       DP=62   GT:AD:DP:RGQ    0/0:62:62:99
241,PGS000662_hmPOS_GRCh38,20,32759707,C,CA,0.046670938,additive,,,,,,,,,,,,,,unmatched,my-wgs

chrX    51502425        .       A       GA      .       .       DP=19   GT:AD:DP:RGQ    0/0:0:19:0
264,PGS000662_hmPOS_GRCh38,X,51502425,GA,G,0.099701455,additive,,,,,,,,,,,,,,unmatched,my-wgs

Any ideas what could be the reason?

[smlmbrt]

For the first variant the scoring file is for a C/CA allele combination (INDEL) but your VCF is an A/C (SNP). Same for your second example.

[mfasold]

What would be the correct way to score these InDels when I have a hom-ref in the sample? Is it known if it is a deletion or insertion (compared to the reference) based on the effect_allele or other_allele entries in the scorefile? As far as I can see, the InDels will never be matched on hom-ref VCF entries.

[smlmbrt]

That's probably true. It's important to note that the pipeline has been setup and optimised for the scenario where you have array genotypes and submitted the sample to an imputation server so all dosages are recovered and used for GWAS (and subsequent PGS calculation). We do intend to make tools/optimise the pipeline for WGS [g]VCFs, but this is a long-term goal and requires significant re-factoring of the scoring and matching.

[PCyph]

@mfasold Would you mind sharing the code that adds artificial ALT alleles to the VCF file which correspond to the effect alleles in the scoring file? I was just about to implement this myself.

[mfasold]

Absolutely. Please give me 1-2 weeks to do some additional checks. I still see some quite low matched % with some scores. It's only one-liners with a lot of cases not considered. After this, I will post them here