Individual / translational use case.

[bgulko]

First, thank you for this profound piece of work. It holds great promise for work like mine, but I have struggled a bit to configure the pgs_calc reproducible pipeline for my use case (individual WGS, conditioned on models). It seems like the current pipeline focus is on larger samples of highly confident (hard) calls, while own use case (and I expect a growing interest) will be on risk analysis of individual genomes.

My case has a small number of samples (less than 10, and generally just one) and may have a softer (more probabilistic) assessment of each variant. Such analyses would be conditioned on each score (model) and ultimately seek risk projections and variant risk attributions rather than broader demographic distributions of scores.

In particular, with a very small sample size, there seem a number of issues that make this pipeline more difficult to use. I'd be happy to discuss (and work with) this further if there is interest in developing more support for this focus. I don't mean to diminish your efforts here at all.

Below are some issues I (and perhaps others) have encountered – any feedback or approaches I have missed (quite likely) are welcome and would be appreciated:

• The likely number variants from the reference will be small which may induce issues in the "combine matches" as the pipeline may report "Error: no target variants match any variants" (confident reference calls at a position are common but seem excluded from the count). I understand this can be addressed by reducing the --min_overap argument. However, doing this may disable checking for reference build mismatches, which is also important.
• As a related concern, accepting the gVCF format as an input would be most helpful as it compactly covers an entire genome by combining long stretches of confident reference calls into a single record (common for individual genomes). This is more important when the sample size small (or singular). Right now I create a filter from the union of all input to all models and force a hard call at each position in the filter. Perhaps PGS Catalog could offer such a filter?
• There are a number of variants that are not lifted-over during the harmonization process, that might be assigned a harmonized coordinate (for example in GRCh38, and eventually T2T's CHM13). This happens particularly in variants with a valid rsid tags but an "Unknown" source reference tag or “NR” genome_build tag, though it also occurs in some novel variants. For example PGS003753_hmPOS_GRCh38.txt.gz rs147937400 (T>C), which seems to correspond to an original reference of GRCh37 (chr6 29068745), but appears unmapped in the PGS Catalog’s GRCh38 harmonized score file (UCSC reports it mapped to chr6 29100968 in GRCh38)
• Scores (models) may have identifiable source references (via rsid positions or reference nucleotides, required by –liftover option) but these references may not be annotated in the harmonized files. For example PGS003753_hmPOS_GRCh38.txt.gz seems to report #genome_build=NR , but based on the locations assigned to rsIDs, seems to originate from GRCh37.

There are a few other difficulties if they are of interest, but primarily Thank you for your work here! The PGS catalog seems a splendid resource for supporting this growing area.

[smlmbrt]

Hi @bgulko, thanks for your comments on your experiences using the pipeline and suggestions. We're currently busy finalising an imminent software release so will address your queries after that has been finished.

[bgulko]

Bump! I am writing code to handle "NR" builds in models and also unmapped variants in harmonized files, but this is already done somewhere, I'd love to avoid a redundant effort!

[smlmbrt]

Hi @bgulko, what do you mean handle NR? If the PGS had an unknown build and the variants have rsIDs we use a harmoniser (HmPos command from this tool: https://github.com/PGScatalog/pgscatalog_utils) which also does liftover if the build is known. We do plan to use the variant matching utilities of our pgscatalog_utils to do this more efficiently than what we did in the past.

[smlmbrt]

For PGS003753_hmPOS_GRCh38.txt.gz, the NR doesn't seem to be a problem because it uses the rsIDs to fill positions in 38 (only 104/35445 seem to not be mapped, and we'll look into this) but the pipeline then uses hm_chr/hm_pos for variant matching.

[bgulko]

I appreciate your attention and the use of rsid tags (or HLA identifiers) where positions are missing is what I expected from harmonized files. I was surprised to find otherwise (my apologies for the mistake in the example).

In some harmonized PGRS files, the covariates have no label or a non-rs SNP label (suggesting a custom variant included by the author) AND no specified build (the model header carries the genome_build tag value of "NR") and there is no position value in the harmonized file (hm_chr/hm_pos fields are simply blank). It seems like harmonizer software is failing here.

While the genome_build build might be inferred by the positions of other rsid tags in the same model, it seems contrary to the goal of a harmonization to a specific build (in this case grch38). Here is a better example:

File: PGS002746_hmPOS_GRCh38.txt.gz
#genome_build=NR
rsID chr_name chr_position effect_allele other_allele effect_weight OR hm_source hm_rsID hm_chr hm_pos hm_inferOtherAllele
6 150268228 A G -0.17800289564993818 0.83694 Unknown

or

rsID	chr_name	chr_position	effect_allele	other_allele	effect_weight	OR	hm_source	hm_rsID	hm_chr	hm_pos	hm_inferOtherAllele
	6	150268228	A	G	-0.178002896	0.83694	Unknown

From other (mapped) rsid positions in this file it can be inferred that the genome_build tag or hm_source could be identified as grch37/hg19 (assuming a single build for the model). According to UCSC liftover , that position in the harmonized model for grch38 file would be chr6 149947092 .

Perhaps it is possible to offer the harmonizer a default build to use if the "genome_tag" in the header is "NR" and the "hm_source" field is "Unknown". While this sort of missing data seems unusual, there are several harmonized model files that seem to carry an absence like this.

[smlmbrt]

You're right that we can sometimes infer the build from the scoring file, but it is our practice not to as the genome_build is technically curated directly from the author (e.g. whatever is mentioned in the paper). From the perspective of a resource/Catalog we are reporting exactly what we get and not making assumptions, which faithfully reproduces the paper but does result in some missingness (like the variant you highlight). The harmoniser does have an option to override/specify, if you did that and said it should be considered hg19 it would perform the liftover (however our default is to not do this) - but this is also why we allow users to provide their own scores.

Re: rs147937400: we are currently tracking down where this happened (and believe it might be a small DB synchronisation issue) that we're looking to publish correct files for.

There is also a known issue of some positions that don't liftover with the UCSC chain but do liftover with other tools, but this happens to very few variants.