HLA divergence calculation might not be correct

[steletvinicius]

Hi again.

Could you please clarify how the HLA divergence score is calculated?
What is the math behind it?

I have got this comment from the hla_divergence.R file: https://github.com/slowkow/hlabud/blob/main/R/hla_divergence.R

The divergence is the sum of the distances between each pair of amino acids at each position, divided by the total sequence length.

Take into consideration the position argument You added in response to my request about calculating HLA divergence for HLA allele protein domains (substring of the alignment / subset of columns from the hla_alignment matrix).

Taking as an example the allele pair B44:02+B44:03, I have run some tests and It seems to me that something might be not working like expected.

To summarize:

1. HLA protein sequence length is overestimated as it is considering the number of columns from the hla_alignments object which do not correspond to the correct number of amino acids forming the full HLA protein
   For HLA B, the matrix has 380 columns; however, the HLA-B protein has 362 residues

2. When applying the protein segmentation strategy, a new feature implemented as discussed on the issue HLA divergence restricted to the peptide binding groove #4 ,
   the HLA divergence calculation is not assuming the new potein segment length informed on the positions argument.
   Given that, in the end, we will not control/normalize the sequence divergence accordingly as we expect to have a higher divergence index (normalized by area) on the peptide binding groove (wher the amount of differences is more concentrated) when compared to what is found for complete HLA protein.

Please, check my tests on the gist mentioned below.
You will also find there an IMGT protein alignment example.

https://gist.github.com/steletvinicius/7789538387ac7d94c417b747db2be655

[slowkow]

@steletvinicius I've been unable to work this past week but I'll get back to you soon. Thanks for reporting the issue.

[slowkow]

@steletvinicius Hi Vinicius, I tried to address your very helpful comments. Thank you for taking the time to test this so thoroughly. I really appreciate it! 🤟

Below, you can see my re-run of your code from your gist. I used reprex::reprex(input = 'issues/issue5.R') to generate this output.

As far as I can tell, it seems like I have successfully addressed both of the issues you raised:

• 362 length instead of 380 for HLA-B
• shorter length for the subsequence of HLA-B

Please let me know if you can test the new branch:
main...issue5

If you can think of additional examples that I should test, please let me know. I think I should probably add a few more tests to the test suite to check that the divergence calculations are working as expected, so your contributions are very welcome here. Thank you!

---

library(tidyverse)

# devtools::install_github("slowkow/hlabud@issue5")

library(hlabud)

b_align <- hla_alignments("B", type = 'prot', verbose = TRUE)
#> Reading /Users/ks38/Library/Application Support/hlabud/3.59.0/alignments/B_prot.txt

Calculate overall divergence (whole protein sequence)

dosage(b_align$onehot, c("B*44:02", "B*44:03"))
#>               D156 L156
#> B*44:02:01:01    1    0
#> B*44:03:01:01    0    1

overall_divergence <- hlabud::hla_divergence("B*44:02:01:01,B*44:03:01:01",method = "grantham")
overall_divergence
#> B*44:02:01:01,B*44:03:01:01 
#>                   0.4751381

Distance matrix Grantham

matrix   <- amino_distance_matrix(method = "grantham")
distance <- matrix["D","L"]
distance
#> [1] 172

protein_length <- distance / overall_divergence
protein_length
#> B*44:02:01:01,B*44:03:01:01 
#>                         362

dim(b_align$alleles[c("B*44:02:01:01", "B*44:03:01:01"),])
#> [1]   2 380

However, HLA full length proteins (signal peptide + mature protein) have 362 amino acids

HLA protein domains coordinates

positions_leader_b <- sprintf("n%s", 30:7)
positions_exon2_b  <-c(sprintf("n%s", 6:1), 1:84 )
positions_exon3_b  <- as.character(c(85:176))
positions_pbg_b    <- c(positions_exon2_b, positions_exon3_b)

b_align$alleles[c("B*44:02:01:01", "B*44:03:01:01"),positions_pbg_b]
#>               n6  n5  n4  n3  n2  n1  1   2   3   4   5   6   7   8   9   10 
#> B*44:02:01:01 "L" "T" "E" "T" "W" "A" "G" "S" "H" "S" "M" "R" "Y" "F" "Y" "T"
#> B*44:03:01:01 "L" "T" "E" "T" "W" "A" "G" "S" "H" "S" "M" "R" "Y" "F" "Y" "T"
#>               11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26 
#> B*44:02:01:01 "A" "M" "S" "R" "P" "G" "R" "G" "E" "P" "R" "F" "I" "T" "V" "G"
#> B*44:03:01:01 "A" "M" "S" "R" "P" "G" "R" "G" "E" "P" "R" "F" "I" "T" "V" "G"
#>               27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42 
#> B*44:02:01:01 "Y" "V" "D" "D" "T" "L" "F" "V" "R" "F" "D" "S" "D" "A" "T" "S"
#> B*44:03:01:01 "Y" "V" "D" "D" "T" "L" "F" "V" "R" "F" "D" "S" "D" "A" "T" "S"
#>               43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58 
#> B*44:02:01:01 "P" "R" "K" "E" "P" "R" "A" "P" "W" "I" "E" "Q" "E" "G" "P" "E"
#> B*44:03:01:01 "P" "R" "K" "E" "P" "R" "A" "P" "W" "I" "E" "Q" "E" "G" "P" "E"
#>               59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74 
#> B*44:02:01:01 "Y" "W" "D" "R" "E" "T" "Q" "I" "S" "K" "T" "N" "T" "Q" "T" "Y"
#> B*44:03:01:01 "Y" "W" "D" "R" "E" "T" "Q" "I" "S" "K" "T" "N" "T" "Q" "T" "Y"
#>               75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90 
#> B*44:02:01:01 "R" "E" "N" "L" "R" "T" "A" "L" "R" "Y" "Y" "N" "Q" "S" "E" "A"
#> B*44:03:01:01 "R" "E" "N" "L" "R" "T" "A" "L" "R" "Y" "Y" "N" "Q" "S" "E" "A"
#>               91  92  93  94  95  96  97  98  99  100 101 102 103 104 105 106
#> B*44:02:01:01 "G" "S" "H" "I" "I" "Q" "R" "M" "Y" "G" "C" "D" "V" "G" "P" "D"
#> B*44:03:01:01 "G" "S" "H" "I" "I" "Q" "R" "M" "Y" "G" "C" "D" "V" "G" "P" "D"
#>               107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
#> B*44:02:01:01 "G" "R" "L" "L" "R" "G" "Y" "D" "Q" "D" "A" "Y" "D" "G" "K" "D"
#> B*44:03:01:01 "G" "R" "L" "L" "R" "G" "Y" "D" "Q" "D" "A" "Y" "D" "G" "K" "D"
#>               123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
#> B*44:02:01:01 "Y" "I" "A" "L" "N" "E" "D" "L" "S" "S" "W" "T" "A" "A" "D" "T"
#> B*44:03:01:01 "Y" "I" "A" "L" "N" "E" "D" "L" "S" "S" "W" "T" "A" "A" "D" "T"
#>               139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
#> B*44:02:01:01 "A" "A" "Q" "I" "T" "Q" "R" "K" "W" "E" "A" "A" "R" "V" "A" "E"
#> B*44:03:01:01 "A" "A" "Q" "I" "T" "Q" "R" "K" "W" "E" "A" "A" "R" "V" "A" "E"
#>               155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
#> B*44:02:01:01 "Q" "D" "R" "A" "Y" "L" "E" "G" "L" "C" "V" "E" "S" "L" "R" "R"
#> B*44:03:01:01 "Q" "L" "R" "A" "Y" "L" "E" "G" "L" "C" "V" "E" "S" "L" "R" "R"
#>               171 172 173 174 175 176
#> B*44:02:01:01 "Y" "L" "E" "N" "G" "K"
#> B*44:03:01:01 "Y" "L" "E" "N" "G" "K"

peptide binding groove domain length

segment_length <- dim(b_align$alleles[c("B*44:02:01:01", "B*44:03:01:01"),positions_pbg_b])
segment_length
#> [1]   2 182
length(positions_pbg_b)
#> [1] 182

Calculating pbg divergence

pbg_divergence <- hlabud::hla_divergence("B*44:02:01:01,B*44:03:01:01",method = "grantham",positions = positions_pbg_b)
pbg_divergence
#> B*44:02:01:01,B*44:03:01:01 
#>                   0.9450549

This should be equal to the length of the pbg_b sequence

distance / pbg_divergence
#> B*44:02:01:01,B*44:03:01:01 
#>                         182
length(positions_pbg_b)
#> [1] 182

Should return true as there is only one amino acid difference,
it is located inside the pbg region and the peptide length is shorter

overall_divergence < pbg_divergence
#> B*44:02:01:01,B*44:03:01:01 
#>                        TRUE

Adjusted by peptide binding groove length

pbg_divergence_adjusted <- distance / length(positions_pbg_b)
pbg_divergence_adjusted
#> [1] 0.9450549

pbg_difference <- overall_divergence - pbg_divergence_adjusted
pbg_difference
#> B*44:02:01:01,B*44:03:01:01 
#>                  -0.4699168

^{Created on 2025-04-08 with reprex v2.0.2}