HLA divergence restricted to the peptide binding groove

[steletvinicius]

Once again, thank you for developing and maintaining such a great toolbox.

I am working on calculating HLA divergence from a set of allele pairs.

From the tests I have made, I am assuming that the divergence calculation includes ANY amino acid differences among the allele pair (predicted protein pair), right?

> library(hlabud)
> #Same P Group; different proteins
> hlabud::hla_divergence("A*01:02:01:01,A*01:412")
A*01:02:01:01,A*01:412 
            0.05511811 
> 
> #Same P Group; same protein
> hlabud::hla_divergence("A*01:02:01:01,A*01:02:02")
A*01:02:01:01,A*01:02:02 
                       0 
> 
> #Different P Groups
> hlabud::hla_divergence("A*01:02:01:01,A*01:315")
A*01:02:01:01,A*01:315 
             0.8057743 

Is there any possible way of restricting this calculation to amino acid differences located ONLY at the regions forming the peptide binding groove (exon 2 +3, class I genes; exon 2, class II genes)?

[slowkow]

Hey Vinicius, thanks for opening this interesting issue!

the divergence calculation includes ANY amino acid differences among the allele pair (predicted protein pair), right?

That's right, the entire length of the sequence is used for the calculation, just like the original code by Pierini & Lenz 2018 MolBiolEvol.

And we can see a bit more in ?hla_divergence:

     When comparing two amino acid sequences, only characters that are
     one of the 20 amino acids are considered in the divergence
     calculation, so gaps (and any other characters) do not count.

Is there any possible way of restricting this calculation to amino acid differences located ONLY at the regions forming the peptide binding groove (exon 2 +3, class I genes; exon 2, class II genes)?

Right now, there isn't any straightforward way to quickly change the hlabud code to support calculation over just the peptide binding groove amino acid positions. But this is definitely something that should be possible. So, I'll work on it and see how it goes. Thanks for the request.

Could I please ask if you might be able to share a link to your favorite publication that clearly defines which portions of the sequence are part of the peptide binding groove (or other sequence features)? I'm always looking for more references.

To implement this feature:

• First, I'll have to try to figure out how to retain the exon information from the original input files, so we can reference it later.
• Second, I'll try to modify the divergence function so we can have some control over what positions are part of the calculation.

This will take some time, so I'll try to post an update here when I eventually get around to this. Pull requests are welcome, if anyone wants to contribute.

[steletvinicius]

There is a whole book available online describing many aspects about HLA molecules

The HLA FactsBook
https://www.sciencedirect.com/book/9780125450256/the-hla-factsbook

Chapters 04 and 05 describe the structure molecule structure in details and show the correspondence exon number - protein domain

This is also a classic review about the HLA genes function and structure:
https://pubmed.ncbi.nlm.nih.gov/10974135/

I could not find any reference to the ordinary genomic coordinates for the exon frontiers.

But, checking the information from the protein sequence alignments on IPD-IMGT/HLA database website,
if we select only the signal peptide - derived from exon 1 sequence - we will find out that:

For classical HLA class I alleles (HLA-A, -B, -C):

1. signal peptide is formed by the first 24 amino acids (coordinates -24 to -1 );
2. peptide binding groove extends from amino acid residues +1 to 182, composed by two domains:
   2.1 alpha 1 domain - derived from exon 2 (formed by the first 90 amino acids [coordinates 1 to 90])
   2.2 alpha 2 domain - derived from exon 3 (formed by the next 92 amino acids [coordinates 91 to 182])

What if when calculating peptide binding groove divergence, the algorithm could take the partial protein sequence based on those coordinates (amino acid coordinates 1 to 182 like shown on the image below ) and them submit them to the function responsible for the divergence calculation?

[slowkow]

The HLA FactsBook

Thank you for the links! I didn't know about the HLA FactsBook. I will eventually try to read this.

This is also a classic review

And thank you so much for the link to the review! The figures are absolutely gorgeous, and there are so many incredibly interesting references to follow. It's such a pleasure to explore, thank you. However, I'm not seeing an explicit mention of which amino acid positions are the most important for defining the floor and walls of the peptide binding groove in each gene.

checking the information from the protein sequence alignments on IPD-IMGT/HLA database website

Could I please ask if you can share the link to the page you're viewing? I browsed around a bit, but I am not seeing a page with the coordinates you mentioned. Maybe I'm looking in the wrong place.

I can go to PDB and manually figure out the correspondence of which positions in the IMGT files correspond to which positions in the PDB structure, and go from there. But I haven't yet spent enough time on this.

What if when calculating peptide binding groove divergence, the algorithm could take the partial protein sequence based on those coordinates (amino acid coordinates 1 to 182 like shown on the image below ) and them submit them to the function responsible for the divergence calculation?

Sure, I made a new branch called divergence_positions that implements your suggestion. Let me know what you think, and we can keep going from there.

Here's an example that seems to work:

# Install the new branch
# devtools::install_github("slowkow/hlabud@divergence_positions")

hla_divergence("A*01:02:01:01,A*01:315")
#> A*01:02:01:01,A*01:315 
#>              0.8057743
hla_divergence("A*01:02:01:01,A*01:315", positions = 1:182)
#> A*01:02:01:01,A*01:315 
#>              0.7217848
hla_divergence("A*01:02:01:01,A*01:315", positions = 1:90)
#> A*01:02:01:01,A*01:315 
#>              0.6955381
hla_divergence("A*01:02:01:01,A*01:315", positions = 91:182)
#> A*01:02:01:01,A*01:315 
#>             0.02624672

[steletvinicius]

Could I please ask if you can share the link to the page you're viewing? I browsed around a bit, but I am not seeing a page with the coordinates you mentioned. Maybe I'm looking in the wrong place.

This is the link for the IMGT alignment tool https://www.ebi.ac.uk/ipd/imgt/hla/alignment/

After accessing, if you choose the "Protein - Signal Petide" option and run the alignment, you will take as a result for each protein, the first 24 peptides from the IMGT alignment files available on github.

This is the leader peptide and it does not take part on the final protein. (Although there are many evidences suggesting its relevance on the HLA-E protein modulatory role over NK-cells)

However, I'm not seeing an explicit mention of which amino acid positions are the most important for defining the floor and walls of the peptide binding groove in each gene.

Unfortunatelly, I could not find any cannonical reference to the coordinates I mentioned as boundaries for the alpha 1 and alpha 2 domains. I took them from a NGS read alignment alignment software that follows the IMGT official coordinates.

These coordinates corresponds to the amino acids encoded by exons 2 and 3, which corresponds to the alpha 1 and 2 domains of the classical class I genes.

We have also this interesting web tool that organizes 3D structures from HLA molecules.
Sharing if it might be of some help.
[https://www.phla3d.com.br/alleles/index ]

Sure, I made a new branch called divergence_positions that implements your suggestion. Let me know what you think, and we can keep going from there.

Thank you. I am going to perform some experiments and share the results.

Here's an example that seems to work:

These examples are taking substrings from the fasta protein sequence and them performing the divergence calculation, right?
So, if I want do remove the leader peptide - first 24 amino acids - and include the following 182 amino acids (peptide binding groove coordinates, if I am correct), I should inform:

hla_divergence("A*01:02:01:01,A*01:315", positions = 25:206)

To be sure if my calcs and the boudnaries are correct, could you add an option to this new feature?
If it could return not only the divergence calculation but also the peptide string selected by the user on the interval he could check by its own if the boundaries are correct. Not sure if it would be good for production, but for development and validation, it would help me.

[slowkow]

https://www.phla3d.com.br/alleles/index

Thanks for sharing this resource, wonderful!

I found this in the HLA Factsbook:

I might go ahead and try to work with PDB files in the future, unless I can find some other published resource with more explicit descriptions like this one.

---

These examples are taking substrings from the fasta protein sequence and them performing the divergence calculation, right?

The new positions argument accepts column names from the a$alleles matrix:

> a <- hla_alignments("A")

> a$alleles[1:5,1:35]
               n26 n25 n24 n23 n22 n21 n20 n19 n18 n17 n16 n15 n14 n13 n12 n11 n10 n9  n8  n8_n7
A*01:01:01:01  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."
A*01:01:01:02N "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."
A*01:01:01:03  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."
A*01:01:01:04  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."
A*01:01:01:05  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."
               n7  n6  n5  n4  n3  n2  n2_n1 n1  1   2   3   4   5   6   7
A*01:01:01:01  "T" "Q" "T" "W" "A" "G" "."   "S" "H" "S" "M" "R" "Y" "F" "F"
A*01:01:01:02N "T" "Q" "T" "W" "A" "G" "."   "S" "H" "S" "M" "R" "Y" "F" "F"
A*01:01:01:03  "T" "Q" "T" "W" "A" "G" "."   "S" "H" "S" "M" "R" "Y" "F" "F"
A*01:01:01:04  "T" "Q" "T" "W" "A" "G" "."   "S" "H" "S" "M" "R" "Y" "F" "F"
A*01:01:01:05  "T" "Q" "T" "W" "A" "G" "."   "S" "H" "S" "M" "R" "Y" "F" "F"

So, for example, we can run it like this:

> my_positions <- sprintf("n%s", 26:1)

> my_positions
 [1] "n26" "n25" "n24" "n23" "n22" "n21" "n20" "n19" "n18" "n17" "n16" "n15" "n14" "n13" "n12" "n11"
[17] "n10" "n9"  "n8"  "n7"  "n6"  "n5"  "n4"  "n3"  "n2"  "n1"

> hla_divergence("A*01:02:01:01,A*01:315", positions = my_positions)
A*01:02:01:01,A*01:315
             0.0839895

I am not sure if I like the column names like "1", "2", etc. I think there is too much confusion between integer indexing with literal numbers like 1 and character numbers like "1". So, I think maybe the convention should be to name the positions like "nX" or "pX" where X is an integer. There are some awkward position names like "n2_n1" when there are gaps in the alignment, and I hope that is clear enough. But I'm open to suggestions.

If you pass invalid positions, it should stop with an error like this:

> bad_positions <- c("a", "b", "-1", "wrong")

> hla_divergence("A*01:02:01:01,A*01:315", positions = bad_positions)
Error: positions must be a subset of column names from the hla_alignments()$alleles matrix:
n26,n25,n24,n23,n22,n21,n20,n19,n18,n17 ... 339,340,341,342,343,344,345,346,347,348

The function is automatically converting the values given in positions with as.character(), so these two calls are identical:

> hla_divergence("A*01:02:01:01,A*01:315", positions = 1:90)
A*01:02:01:01,A*01:315
             0.6955381

> hla_divergence("A*01:02:01:01,A*01:315", positions = as.character(1:90))
A*01:02:01:01,A*01:315
             0.6955381

[steletvinicius]

On my last message I was planning to share this table you foudn but I thought It could bring more confusion on what we are planning here. The option to investigate divergence only on the amino acid forming the pocket residues could be a new inprovenment for the hla_divergence() function.

I think the greatest challenge here is to match the hlabud alleles matrix columns numbers and the IPD-IMGT/HLA coordinates like I mentioned on my last closed issue.

Despite the fact you implemented the matrix column naming using the prefix n and p (negative and positeive, respectivelly),
based on my experimetns, there are some small modifications to be made:

If you compare
the image I have shared showing the signal peptide alignment and
the alignment matrix generated by hlabud that you shared, you will notice that
the IMGT alignment coordinates ranges from -24 to -1,
while the same region on the hlabud matrix is n26 to n3.

AND

If I mange to take the first 24 nucleotides, using the interval approach It won`t be feasible beacause sometimes you might have INDEL columns like in this example. So, If I run:

a_align$alleles[1:5, 1:24]

 a_align$alleles[1:5, 1:24]
               n26 n25 n24 n23 n22 n21 n20 n19 n18 n17 n16 n15 n14 n13 n12 n11 n10 n9  n8  n8_n7 n7  n6  n5  n4 
A*01:01:01:01  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."  "T" "Q" "T" "W"
A*01:01:01:02N "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."  "T" "Q" "T" "W"
A*01:01:01:03  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."  "T" "Q" "T" "W"
A*01:01:01:04  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."  "T" "Q" "T" "W"
A*01:01:01:05  "M" "A" "V" "M" "A" "P" "R" "T" "L" "L" "L" "L" "L" "S" "G" "A" "L" "A" "L" ".."  "T" "Q" "T" "W "

I will miss the last peptide from this region (column n3, that should be named n1, according to IMGT coordinates)

One suggestion: for the indel columns, instead of naming them like n8_n7, You can use the nomenclature n8_1, which would be the corresponding for -8:1 that is not appropriate for matrix column names.

After that, we can use the argument position as reference for slicing the peptide and feed the divergence calculation function.
However, if we take the position input as a range of columns, it will always have problems selecting correct columns when we have INDEL columns.

In this case, when the user informs
position = 1 : 25
the function should take the amino acids from columns p1 until p25 (it does not matter how many columns are in between) and them use these two substrings (peptides) for divergence calculation.

It could be an option for this position argument: instead of the range, the user could inform the coordinates from IMGT.
(But HLAbud has to stick to IMGT coordinates, as I mentioned)

===========================================================================

Take a look at these images that make the correspondence HLA gene region - HLA molecule domain

For HLA class I genes, the peptide binding groove is formed by information from exons 2 and 3, enconding the alpha 2 and 3 domains

If you access that same alignment tool from IMGT, choosing only exon 2 and 3 for HLA-A and a few alleles,
you will get this result that shows:

• in green the exon 2 region that encodes the first 90 codons enconding the 90 amino acids of the alpha 1 domain;
• in red the exon 3 region that encodes the next 92 codons that will inform the 92 amino acids of the alpha 2 domain;

I believe it might be enough to stablish the exact coordinates for these regions, right?

[slowkow]

I think the greatest challenge here is to match the hlabud alleles matrix columns numbers and the IPD-IMGT/HLA coordinates like I mentioned on my last closed issue.

I thought I had solved the problem in the previous issue, but I was wrong. I didn't run enough tests, so there are still some lingering issues with numbering for some genes. Thank you again for catching this!

Going forward, I will need to be a lot more careful about this. The files are surprisingly difficult to parse correctly...

One suggestion: for the indel columns, instead of naming them like n8_n7, You can use the nomenclature n8_1, which would be the corresponding for -8:1 that is not appropriate for matrix column names.

What does the :1 mean in this name? Is there any reference that supports this naming convention?

In this case, when the user informs position = 1 : 25 the function should take the amino acids from columns p1 until p25 (it does not matter how many columns are in between) and them use these two substrings (peptides) for divergence calculation.
It could be an option for this position argument: instead of the range, the user could inform the coordinates from IMGT.

In the new branch, I designed the function such that these three calls are identical:

# 3 ways to do the same thing
hla_divergence(..., position = 1:5)
hla_divergence(..., position = as.character(1:5))
hla_divergence(..., position = c("1", "2", "3", "4", "5"))

The function takes the values in positions and uses them to select the columns in a$alleles with those names.

IMAGE with green and red boxes

Do you know why there is a | in the string G|GT? What does this indicate? I don't see any mention of the | character here: https://www.ebi.ac.uk/ipd/mhc/alignment/help/

[steletvinicius]

What does the :1 mean in this name? Is there any reference that supports this naming convention?

I am not aware of any reference supporting this nomenclature. I have taken this from a HLA read alignment software.
For each HLA gene, there is one reference sequence and all of the other sequences for the same gene are aligned to this reference
and represented in terms of the diferences among them and the reference. Like we can found on the alignment files on the IMGT github repo.

On the image below, you can see what I mean with an example. In this case, A*01:01:01:01 is the reference allele.

The function takes the values in positions and uses them to select the columns in a$alleles with those names.

I got it. As a workaround, I have made my own research comparing the IPD protein alignment and the matrix generated by the hla_alignments( ) function. After that, I believe I was able to find and determine the specific coordinates (matrix columns) for each relevante region of the HLA genes (protein), like the peptide binding groove and each anchor pocket.

Could you please check if it make sense? My code is extremelly verbose by I believe it is working fine from the tests I made

https://gist.github.com/steletvinicius/4913c02676a97cd82df353e8cbb156b9

Do you know why there is a | in the string G|GT? What does this indicate? I don't see any mention of the | character here:

The | (pipe) is the genomic region delimiter. In this example, it is showing the nucleotides on the borders among exon 2 | exon 3

I have found this link on the IPD website that might be of some help.
https://www.ebi.ac.uk/ipd/imgt/hla/about/help/genomics/

[slowkow]

I have taken this from a HLA read alignment software.

Thanks for the screenshot! It is very helpful to see other HLA software.

It seems that they have numbered each residue inside the indel with 56:1 for the first residue and 56:2 for the second one, and so on. This is an interesting way to do it, and I'll consider adopting the same convention in the future.

It seems a bit odd to me, because there is a chance that a future undiscovered allele could disrupt this naming convention, resulting in a hypothetical situation where the position we label 56:1 today might be labeled 56:2 tomorrow. I wonder if this is something that we should worry about?

For now, I'm using the convention 56_57 to indicate that there are one or more residues at this indel. I suppose this might be inconvenient if we want easy access to each residue inside an indel, because the user would need to write more code to retrieve each position inside the indel.

The | (pipe) is the genomic region delimiter. In this example, it is showing the nucleotides on the borders among exon 2 | exon 3

I would have expected the pipe | to separate the exons, but in your screenshot you drew the red box including one nucleotide to the left of the pipe. Indeed, I checked that exon 3 does start with codon 91 G|GT, encoding the amino acid Gly (G). I would have expected the pipe to be on the left of the first G like |GGT instead of G|GT. I don't understand why the pipe is positioned where it is. Maybe there is some text about this somewhere? Any tips are appreciated!

Could you please check if it make sense?

Please try it again with the latest code from commit 28ae7b7. This is my latest attempt to number the positions correctly. Before this commit, some genes were numbered correctly and others were incorrect. After this commit, I am confident that the positions are named correctly, because I manually verified each gene and also wrote some tests to automatically verify that the positions are named correctly.

However, there is always a chance that I might have missed something! So, I would be grateful if you would report any issues that you find. Thank you in advance for your patience! 🙏

[steletvinicius]

It seems a bit odd to me, because there is a chance that a future undiscovered allele could disrupt this naming convention, resulting in a hypothetical situation where the position we label 56:1 today might be labeled 56:2 tomorrow. I wonder if this is something that we should worry about?

The only change you will have to consider is extended this insertion region length if another allele is reported with 7 or more insertion bases. In this case, you will have the coordinates 56:1 56:2 56:3 56:4 56:5 56:6 56:7 - or more.

I would have expected the pipe | to separate the exons, but in your screenshot you drew the red box including one nucleotide to the left of the pipe. Indeed, I checked that exon 3 does start with codon 91 G|GT, encoding the amino acid Gly (G). I would have expected the pipe to be on the left of the first G like |GGT instead of G|GT. I don't understand why the pipe is positioned where it is. Maybe there is some text about this somewhere? Any tips are appreciated!

I am not sure if it happens on other than HLA genes but the first codon of each exon is formed by:
the last base of the previous exon + the first two bases of the current exon.

Please try it again with the latest code from commit 28ae7b7. This is my latest attempt to number the positions correctly. Before this commit, some genes were numbered correctly and others were incorrect. After this commit, I am confident that the positions are named correctly, because I manually verified each gene and also wrote some tests to automatically verify that the positions are named correctly.

I will give it a try latter cause I am busy running some statistics on the divergence data I generated using that customized solution described on the gist that I came up with.

Thanks in advance for all support.