Protocol for defining locus coordinates and genotyping definitions

[hdashnow]

I'm starting a conversation about general strategies to define locus coordinates after useful discussions with @laurelhiatt @avvaruakshay @penhoorn @egor-dolzhenko @bw2. I wanted to give you all a chance to discuss this before we implement anything.

All genomic coordinates in STRchive should be 0-based (0-start, half-open). See here for more details.
Coordinates will need to be adjusted for linking with UCSC, see #149.

Defining the "core" STRchive locus coordinates:

• Start with the coordinates or region defined in the literature
• Use TRF annotations or similar to refine the coordinates if needed
• Examine the reference genome at these coordinates and find the longest exact stretch of the motif (allowing for "N" bases if relevant, e.g. polyalanine loci)
• The coordinates should include partial motifs i.e. the total length does not need to be a multiple of the motif length
• The coordinates may additionally include interrupted motifs if they are described in the literature
• The coordinates should first be defined in one reference genome (typically hg38), and then lifted over to the others. After leftover, confirm that the coordinates still meet the above criteria and adjust them as needed.

The reasoning:
This is designed specifically to provide the coordinates of the pathogenic part of the sequence.

Defining coordinates for TRGT genotyping

• Start with the "core" STRchive locus coordinates
• Add any flanking repeat or known variable regions
• Include all known reference, benign, and pathogenic motifs in the motif field (how should these be ordered?)

Defining coordinates for ExpansionHunter genotyping

• Start with the "core" STRchive locus coordinates
• Limit the region to the longest exact stretch of the motif (allowing for "N" bases if relevant, e.g. polyalanine loci)
• Exclude 3' (downstream) partial motifs (unsure if this is needed?)

The reasoning:
TRGT performs better when flanking variable regions are included. ExpansionHunter conversely performs better when only pure regions are included.

[egor-dolzhenko]

This sounds great to me! I'd only maybe consider changing "The coordinates should include partial motifs" to "The coordinates could include partial motifs".

[hdashnow]

Does the presence of partial motifs affect either the TRGT genotyping or the new motif counting algorithm? I'd like to choose a strategy and use it uniformly if possible.

[laurelhiatt]

@egor-dolzhenko if you don't feel that the presence of partial motifs is an issue, I can correct all the STRchive coordinates based on these principles!

[egor-dolzhenko]

@hdashnow: Sorry for the late reply. The new motif counting algorithm will most likely be included in TRGT v2.0.0 (which will hopefully be ready for release in the coming months). With this new algorithm, partial STR motif matches will be excluded from motif count (but will still appear in the allele sequence).

@laurelhiatt: It would be great to do this, thank you! Maybe the updated repeat catalog could be released at the same time as TRGT v2?

[hdashnow]

@egor-dolzhenko that's great!
I think we will include partial motifs for all loci then. I think it best reflects the underlying biology of instability.

We will be working on updating the STRchive loci. This particular fix is pretty easy but the motif stuff might be a bit more involved. Keep us posted on the TRGT v2 timing and we can try to match up the timeline if needed :)

[egor-dolzhenko]

This sounds great @hdashnow!

[penhoorn]

I just want to echo @egor-dolzhenko comment. The plan of action looks great. I would be happy to "test drive" the new coordinates for TRGT once they're available.

[bw2]

At first glance, this all sounds reasonable to me also.
Unfortunately, at least for the existing known loci, I think there are remaining issues.
In addition to published literature, there are a lot of established resources that have already used and will likely continue to use the alternative locus definitions from the original ExpansionHunter repo (including DRAGEN ExpansionHunter calls now available to all researchers in the UKBB RAP, etc.).

I just ran a 4-way comparison of the locus definitions in the

1. gnomAD catalog (after adding the latest set of loci recently included in STRchive) vs.
2. the definitions in the STRchive json vs.
3. STRipy definitions vs.
4. the original ExpansionHunter catalog

Here are the full results of the comparison.

More than half (38 out of 70) of the loci present in both the gnomAD and STRchive jsons currently have different start/end coordinates(!) Many loci differ by only 1 or 2 bases, but a few differ by a lot.

For example, our PABPN1 definitions differ by 4 repeats, with STRchive choosing the longer definition while the other sources use the shorter one (we discuss these different competing definitions in the Tandem Repeat Catalogs preprint):

PABPN1 hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chr14:23321472-23321490  ( 6.0 x GCG     18bp) official catalog from EH repo
    chr14:23321472-23321490  ( 6.0 x GCG     18bp) gnomAD    7 = pathogenic.min.  hg38 seq:             GCGGCGGCGGCGGCGGCG
    chr14:23321472-23321502  (10.0 x GCN     30bp) STRchive  12 = pathogenic.min.  hg38 seq:            GCGGCGGCGGCGGCGGCGGCAGCAGCAGCG
    chr14:23321472-23321490  ( 6.0 x GCG     18bp) STRipy    7 = pathogenic.min.  hg38 seq:             GCGGCGGCGGCGGCGGCG

Large differences also exist for C9orf72, AR, FMR1, JPH3, MARCHF6, AFF2, DAB1, BEAN1, FGF14, etc.

It would be great to resolve these differences as much as possible.
As part of that, what do you think about extending the

The coordinates may additionally include interrupted motifs if they are described in the literature

rule to say that existing definitions in the original Illumina ExpansionHunter catalog should take precedence (given the number of resources that rely on them)?

[hdashnow]

Hi Ben,

We recently updated the coordinates using the protocol discussed above. But were waiting to merge any changes until we'd had these kinds of community discussions.
If you want to give feedback on the proposed new coordinates, they are here (although may still change as this discussion evolves): #155
All the off-by-one errors should be fixed, but there are likely many loci where the inclusion of partial motifs will result in differences.

I'm going to have to give some more thought to the question of consistency. Our goal here was to produce the "best" set of loci (by some definition), and then adapt it to the needs of specific variant callers. I'm not sure how to reconcile that with consistency. I will keep thinking.

Warm regards,
Harriet

[bw2]

I just looked at the updated STRchive catalog and see the improvements.

I think PABPN1 is a good example for this discussion..
AFAIK the narrower locus definition is no better or worse than the wider definition both for short read and long read calling. Some resources and papers use one definition, while ~ the same number use the other. I take your point about not tieing STRchive to one specific STR caller. On the other hand, for situations where the alternative definitions are ~ equally good, and one of them is present in the original ExpansionHunter catalog, I think it provides an appropriate tie breaker because of the number of studies that were/will be based on it.

Redoing the cross-catalog comparison vs. the PR version of STRchive-loci.json, and adding in the TRGT catalog, these are loci where the STRchive and gnomAD catalogs still differ by more than one motif length:

------
AFF2 hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chrX:148500631-148500691  (20.0 x GCC     60bp) official catalog from EH repo
    chrX:148500631-148500691  (20.0 x GCC     60bp) gnomAD    200 = pathogenic.min.  hg38 seq:                                     GCCGCTGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCTGCCGCC
    chrX:148500604-148500753  (49.7 x GCC    149bp) STRchive  201 = pathogenic.min.  hg38 seq:          GCCGCGGCCGCCGCCGCCGCCTGTGCAGCCGCTGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCTGCCGCCCCGGCTGCCGCGCCGCGCCGCTGCCTCTGCCCCGGCCGCCCCCGCCGCCGCTGCCGCCGCCG
    chrX:148500638-148500683  (15.0 x CCG     45bp) STRipy    200 = pathogenic.min.  hg38 seq:                                            CCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCG
    chrX:148500631-148500691  (20.0 x CCG     60bp) TRGT                        hg38 seq:                                          GCCGCTGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCCGCTGCCGCC
------
AR hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chrX:67545316-67545385  (23.0 x GCA     69bp) official catalog from EH repo
    chrX:67545316-67545385  (23.0 x GCA     69bp) gnomAD    40 = pathogenic.min.  hg38 seq:             GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA
    chrX:67545316-67545419  (34.3 x GCA    103bp) STRchive  38 = pathogenic.min.  hg38 seq:             GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAAGAGACTAGCCCCAGGCAGCAGCAGCAGCAGCAG
    chrX:67545316-67545385  (23.0 x CAG     69bp) STRipy    40 = pathogenic.min.  hg38 seq:             GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA
    chrX:67545316-67545385  (23.0 x AGC     69bp) TRGT                        hg38 seq:                 GCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCA
------
ATXN3 hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chr14:92071009-92071042  (11.0 x GCT     33bp) official catalog from EH repo
    chr14:92071010-92071040  (10.0 x CTG     30bp) gnomAD    56 = pathogenic.min.  hg38 seq:            CTGCTGCTGCTGCTGCTGCTGCTGTTGCTG
    chr14:92071010-92071052  (14.0 x CTG     42bp) STRchive  60 = pathogenic.min.  hg38 seq:            CTGCTGCTGCTGCTGCTGCTGCTGTTGCTGCTTTTGCTGCTG
    chr14:92071010-92071040  (10.0 x CAG     30bp) STRipy    56 = pathogenic.min.  hg38 seq:            CTGCTGCTGCTGCTGCTGCTGCTGTTGCTG
    chr14:92071009-92071042  (11.0 x AGC     33bp) TRGT                        hg38 seq:                CCTGCTGCTGCTGCTGCTGCTGCTGTTGCTGCT
------
C9ORF72 hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chr9:27573528-27573546  ( 3.0 x GGCCCC     18bp) official catalog from EH repo
    chr9:27573528-27573546  ( 3.0 x GGCCCC     18bp) gnomAD    31 = pathogenic.min.  hg38 seq:                                                      GGCCCCGGCCCCGGCCCC
    chr9:27573484-27573546  (10.3 x GGCCCC     62bp) STRchive  251 = pathogenic.min.  hg38 seq:         GCCCCGCCCCGGGCCCGCCCCCGGGCCCGCCCCGACCACGCCCCGGCCCCGGCCCCGGCCCC
    chr9:27573528-27573546  ( 3.0 x GGGGCC     18bp) STRipy    24 = pathogenic.min.  hg38 seq:                                                      GGCCCCGGCCCCGGCCCC
    chr9:27573528-27573546  ( 3.0 x CCCCGG     18bp) TRGT                        hg38 seq:                                                          GGCCCCGGCCCCGGCCCC
------
JPH3 hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chr16:87604287-87604329  (14.0 x CTG     42bp) official catalog from EH repo
    chr16:87604287-87604329  (14.0 x CTG     42bp) gnomAD    41 = pathogenic.min.  hg38 seq:                 CTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTG
    chr16:87604282-87604329  (15.7 x CTG     47bp) STRchive  40 = pathogenic.min.  hg38 seq:            GCTGCCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTG
    chr16:87604287-87604329  (14.0 x CTG     42bp) STRipy    41 = pathogenic.min.  hg38 seq:                 CTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTG
    chr16:87604287-87604329  (14.0 x AGC     42bp) TRGT                        hg38 seq:                     CTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTGCTG
------
PABPN1 hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chr14:23321472-23321490  ( 6.0 x GCG     18bp) official catalog from EH repo
    chr14:23321472-23321490  ( 6.0 x GCG     18bp) gnomAD    7 = pathogenic.min.  hg38 seq:             GCGGCGGCGGCGGCGGCG
    chr14:23321472-23321503  (10.3 x GCN     31bp) STRchive  12 = pathogenic.min.  hg38 seq:            GCGGCGGCGGCGGCGGCGGCAGCAGCAGCGG
    chr14:23321472-23321490  ( 6.0 x GCG     18bp) STRipy    7 = pathogenic.min.  hg38 seq:             GCGGCGGCGGCGGCGGCG
    chr14:23321472-23321490  ( 6.0 x CCG     18bp) TRGT                        hg38 seq:                GCGGCGGCGGCGGCGGCG
------
RUNX2 hg38 ReferenceRegion isn't the same between gnomAD and STRchive
    chr6:45422750-45422792  (14.0 x GCN     42bp) gnomAD    27 = pathogenic.min.  hg38 seq:             GCGGCGGCGGCGGCTGCGGCGGCGGCGGCGGCTGCGGCGGCG
    chr6:45422750-45422801  (17.0 x GCN     51bp) STRchive  20 = pathogenic.min.  hg38 seq:             GCGGCGGCGGCGGCTGCGGCGGCGGCGGCGGCTGCGGCGGCGGCAGCTGCA
    chr6:45422750-45422801  (17.0 x GCN     51bp) STRipy    20 = pathogenic.min.  hg38 seq:             GCGGCGGCGGCGGCTGCGGCGGCGGCGGCGGCTGCGGCGGCGGCAGCTGCA
    chr6:45422750-45422802  (17.3 x CNG     52bp) TRGT                        hg38 seq:                 GCGGCGGCGGCGGCTGCGGCGGCGGCGGCGGCTGCGGCGGCGGCAGCTGCAG

The full comparison results are here.

[bw2]

I was syncing the latest pathogenic thresholds from STRchive to gnomAD and the same question as for PABPN1 applies to the AR locus. Does STRchive suggest counting the "AGAGACTAGCCCCAGGCAGCAGCAGCAGCAGCAG" sequence when comparing to the pathogenic threshold of 38 x CAG?