In 2024, we put out a preprint that made public an initial version of the CMV ECOcluster. In 2026, having rebuilt ECOclusters with more data and improved algorithms, we updated the ECOclusters preprint significantly and put a new version of the CMV ECOcluster in the supplemental data. The two versions both define sensitive, specific classifiers for CMV, but they’re very different from each other. Here, I’ll characterize the differences and try to explain them a bit.
Since anyone familiar with the 2024 CMV ECOcluster is likely to be a bit of an ECOclusters enthusiast, already, this post will get fairly far into some weeds.
Basic imports.
# useful imports
import pandas as pd
import seaborn as sns
from cmvividly.data import access as cmv_access
from cmvividly.data.access import load_cmv_ecocluster_2026, load_cmv_ecocluster_2024
from matplotlib import pyplot as plt
from matplotlib_venn import venn2
import itables
from cmvividly.plots.style import set_style
set_style()
#stats
from scipy.stats import mannwhitneyuLoad both versions of the CMV ECOcluster and count the rows.
# load the 2024 CMV ECOcluster.
pdf_cmv_ecocluster_2024 = load_cmv_ecocluster_2024()
print(f"Rows in 2024 ECOcluster: {len(pdf_cmv_ecocluster_2024)}")
# load the 2026 CMV ECOcluster.
pdf_cmv_ecocluster_2026 = load_cmv_ecocluster_2026()
print(f"Rows in 2026 ECOcluster: {len(pdf_cmv_ecocluster_2026)}")Rows in 2024 ECOcluster: 26139
Rows in 2026 ECOcluster: 52447First, a note: the 2024 CMV ECOcluster has some TCRs that occur multiple times, while the 2026 CMV ECOcluster does not:
print(f"Max occurrences of a TCR in 2026:", pdf_cmv_ecocluster_2026["tcr"].value_counts().max())
print(f"Occurrence count breakdown in 2024:")
pdf_cmv_ecocluster_2024.tcr.value_counts().reset_index()["count"].value_counts()Max occurrences of a TCR in 2026: 1
Occurrence count breakdown in 2024:count
1 26073
2 33
Name: count, dtype: int64That’s due to methodological differences in TCR-HLA allele association between 2024 and 2026. It’s pretty much a rounding error, but it’ll throw our accounting off slightly, so let’s deduplicate the 2024 version by TCR (taking an HLA association at random).
pdf_cmv_ecocluster_2024_dedup = pdf_cmv_ecocluster_2024.drop_duplicates(subset=["tcr"])
print(f"Rows in deduplicated 2024 dataframe: {len(pdf_cmv_ecocluster_2024_dedup)}")Rows in deduplicated 2024 dataframe: 26106tcr_set_2024 = set(pdf_cmv_ecocluster_2024_dedup["tcr"].unique())
tcr_set_2026 = set(pdf_cmv_ecocluster_2026["tcr"].unique())
all_tcrs_set = tcr_set_2024 | tcr_set_2026
proportion_all_tcrs_shared = len(tcr_set_2024 & tcr_set_2026) / len(all_tcrs_set)
f, ax = plt.subplots()
venn2([tcr_set_2024, tcr_set_2026], set_labels=("2024", "2026"), ax=ax, subset_label_formatter=lambda n: f"{n:,}")
ax.set_title(f"TCRs ({proportion_all_tcrs_shared*100:.0f}% of total shared)", y=.95)
Wow! The 2026 ECOcluster is much bigger, but it loses half the TCRs from the 2024 version.
# build a 2x2 table: # of TCRs in 2024 where hla == "ci" and where hla == "cii", and same in 2026
n_ci_2024 = sum(pdf_cmv_ecocluster_2024_dedup["hla_class"] == "ci")
n_cii_2024 = sum(pdf_cmv_ecocluster_2024_dedup["hla_class"] == "cii")
n_ci_2026 = sum(pdf_cmv_ecocluster_2026["hla_class"] == "ci")
n_cii_2026 = sum(pdf_cmv_ecocluster_2026["hla_class"] == "cii")
# Create a 2x2 table
pdf_hlaclass_by_version = pd.DataFrame({
"Class I": [n_ci_2024, n_ci_2026],
"Class II": [n_cii_2024, n_cii_2026],
"% Class II": [n_cii_2024 / (n_ci_2024 + n_cii_2024) * 100, n_cii_2026 / (n_ci_2026 + n_cii_2026) * 100]
}, index=["2024", "2026"])
# make a stacked barplot
f, ax = plt.subplots(figsize=(5, 4))
pdf_hlaclass_by_version[["Class I", "Class II"]].plot(kind="bar", stacked=True, ax=ax)
ax.set_xlabel("Version")
ax.set_ylabel("Number of TCRs")
ax.set_title("TCRs by HLA Class and Version")
ax.legend(title="HLA Class")
plt.xticks(rotation=0)
plt.tight_layout()
pdf_hlaclass_by_version| Class I | Class II | % Class II | |
|---|---|---|---|
| 2024 | 9242 | 16864 | 64.598177 |
| 2026 | 5323 | 47124 | 89.850706 |
The 2026 ECOcluster is 90% Class-II-associated TCRs, vs. 65% in 2024. That’s a big shift! And it’s consistent with the reported changes in HLA models between the two versions: as described in the preprint, Class II alleles tended to have dramatically more associated TCRs than Class II.
But, also, in absolute terms, the 2026 version actually has many fewer Class I-associated TCRs than the 2024 version did. That’s a bit of a surprise; I’ll look into that some more, further down.
# first, add an hla_gene column so we can break things down by gene later
pdf_cmv_ecocluster_2024_dedup["hla_gene"] = pdf_cmv_ecocluster_2024_dedup.hla.apply(
lambda x: x.split("*")[0])
pdf_cmv_ecocluster_2026["hla_gene"] = pdf_cmv_ecocluster_2026.hla.apply(
lambda x: x.split("*")[0])
# set up an outer join on tcr, to build a dataframe that contains all TCRs in both ECOcluster versions
merge_cols = ["tcr", "cdr3", "vgene", "jgene"]
pdf_2024_merge = pdf_cmv_ecocluster_2024_dedup.copy()
pdf_2024_merge["source_2024"] = True
pdf_2026_merge = pdf_cmv_ecocluster_2026.copy()
pdf_2026_merge["source_2026"] = True
pdf_merged_tcr_outer = pdf_2024_merge.merge(
pdf_2026_merge,
on=merge_cols,
how="outer",
suffixes=("_2024", "_2026"),
indicator=True,
)
# annotate each row with HLA class agreement and where the TCR came from
pdf_merged_tcr_outer["hla_class_agree"] = (
pdf_merged_tcr_outer["hla_class_2024"] == pdf_merged_tcr_outer["hla_class_2026"]
)
pdf_merged_tcr_outer["source"] = pdf_merged_tcr_outer["_merge"].map(
{"left_only": "2024 only", "right_only": "2026 only", "both": "shared"})
pdf_merged_tcr_outer["hla_class_sourced"] = pdf_merged_tcr_outer.apply(
lambda row: row["hla_class_2024"] if row["source"] == "2024 only"
else row["hla_class_2026"] if row["source"] == "2026 only"
else row["hla_class_2024"] if row["hla_class_agree"]
else "disagree",
axis=1
)
pdf_merged_tcr_outer["source_2026"] = pdf_merged_tcr_outer.source_2026.fillna(False).astype(bool)
pdf_merged_tcr_outer["source_2024"] = pdf_merged_tcr_outer.source_2024.fillna(False).astype(bool)
print(f"TCRs in merged table: {len(pdf_merged_tcr_outer)}")TCRs in merged table: 66024I don’t have convenient access to a pGen estimate for the 2024 version, but I can look at how it differs in the 2026 version by HLA class and by whether the TCR was in the 2024 version.
pdf_2026_annot_2024 = pdf_merged_tcr_outer[pdf_merged_tcr_outer.source_2026 == True]
pdf_2026_annot_2024_droplowp = pdf_2026_annot_2024[pdf_2026_annot_2024.tcr_pgen >= 1e-25]
f, ax = plt.subplots()
sns.boxplot(data=pdf_2026_annot_2024_droplowp,
x="log10_tcr_pgen_eps",
y="hla_class_2026", hue="source_2024", ax=ax, orient="horiz")
ax.set_xlabel("log10(pGen + eps)")
ax.set_ylabel("HLA Class")
ax.set_title("pGen of 2026 TCRs by HLA class and 2024 presence")
plt.tight_layout()
pdf_2026_annot_2024_cii = pdf_2026_annot_2024[pdf_2026_annot_2024["hla_class_2026"] == "cii"]
pdf_2026_annot_2024_ci = pdf_2026_annot_2024[pdf_2026_annot_2024["hla_class_2026"] == "ci"]
higher_pgen_hlaclass = "Class II" if pdf_2026_annot_2024_cii["tcr_pgen"].mean() > pdf_2026_annot_2024_ci["tcr_pgen"].mean() else "Class I"
mwu_p_hlaclass = mannwhitneyu(
pdf_2026_annot_2024_cii["tcr_pgen"],
pdf_2026_annot_2024_ci["tcr_pgen"],
).pvalue
mwu_p_cii = mannwhitneyu(
pdf_2026_annot_2024_cii[pdf_2026_annot_2024_cii["source_2024"] == True]["tcr_pgen"],
pdf_2026_annot_2024_cii[pdf_2026_annot_2024_cii["source_2024"] == False]["tcr_pgen"],
).pvalue
mwu_p_ci = mannwhitneyu(
pdf_2026_annot_2024_ci[pdf_2026_annot_2024_ci["source_2024"] == True]["tcr_pgen"],
pdf_2026_annot_2024_ci[pdf_2026_annot_2024_ci["source_2024"] == False]["tcr_pgen"],
).pvalue
print(f"Two-sided Mann-Whitney U test p-values:")
print(f"Between HLA classes: {mwu_p_hlaclass:.2e} ({higher_pgen_hlaclass} higher)")
print(f"Within HLA class, 2024 vs. 2026: cii = {mwu_p_cii:.2e}, ci = {mwu_p_ci:.2e}")
print("")
print("(dropping very-low-pGen outliers for visibility)")Two-sided Mann-Whitney U test p-values:
Between HLA classes: 1.22e-03 (Class I higher)
Within HLA class, 2024 vs. 2026: cii = 0.00e+00, ci = 1.84e-78
(dropping very-low-pGen outliers for visibility)
Overall, the Class II TCRs have slightly (and significantly) lower pGen. But the much bigger effect is that, in both HLA classes, the TCRs that are unique to the 2026 version have strongly lower pGen than the ones the two versions share. One way to interpret that is that the much larger dataset that built the new ECOclusters provides more statistical power to associate TCRs that are observed in fewer people.
That leaves the question: why are more than half of the 2024 TCRs missing from the 2026 version? I have a hypothesis that it’s not because they’re not actually CMV-specific. I’ll explore that further down in this post, and likely return to that question in a later post.
First let’s visualize how TCRs are either present or absent in each version, and associated with HLA Class I or Class II alleles in each version, with a Sankey plot.
from cmvividly.plots.sankey import plot_sankey
plot_sankey(pdf_merged_tcr_outer.fillna("MISSING"), source_col="hla_class_2024",
target_col="hla_class_2026")
passThat gives us a feel for what’s going on: HLA class is ~100% consistent when present in both versions, but there are a lot of TCRs just missing from one version. Almost all the new ones in 2026 are Class II, while the ones unique to 2024 are an even split.
Let’s look at that a little more precisely.
source_hlaclass_breakdown = pd.crosstab(pdf_merged_tcr_outer.hla_class_sourced, pdf_merged_tcr_outer.source, margins=True)
source_hlaclass_breakdown.loc["% Class I"] = (source_hlaclass_breakdown.loc["ci"] / source_hlaclass_breakdown.loc["All"] * 100).map("{:.0f}%".format)
source_hlaclass_breakdown.loc["% Class II"] = (source_hlaclass_breakdown.loc["cii"] / source_hlaclass_breakdown.loc["All"] * 100).map("{:.0f}%".format)
source_hlaclass_breakdown| source | 2024 only | 2026 only | shared | All |
|---|---|---|---|---|
| hla_class_sourced | ||||
| ci | 6940 | 3016 | 2301 | 12257 |
| cii | 6637 | 36902 | 10221 | 53760 |
| disagree | 0 | 0 | 7 | 7 |
| All | 13577 | 39918 | 12529 | 66024 |
| % Class I | 51% | 8% | 18% | 19% |
| % Class II | 49% | 92% | 82% | 81% |
There’s a lot in that table:
hla_set_2024 = set(pdf_cmv_ecocluster_2024["hla"].unique())
hla_set_2026 = set(pdf_cmv_ecocluster_2026["hla"].unique())
all_hlas_set = hla_set_2024 | hla_set_2026
proportion_all_hlas_shared = len(hla_set_2024 & hla_set_2026) / len(all_hlas_set)
f, ax = plt.subplots()
venn2([hla_set_2024, hla_set_2026], set_labels=("2024", "2026"), ax=ax, subset_label_formatter=lambda n: f"{n:,}")
ax.set_title(f"Alleles ({proportion_all_hlas_shared*100:.0f}% of total shared)", y=.95)
pass
Let’s look at the specific alleles that drop out or show up in 2026 vs. 2024.
def get_hlalist_by_gene(hlas) -> pd.DataFrame:
genes = set([hla.split("*")[0] for hla in hlas])
gene_list_map = {gene: sorted([hla for hla in hlas if hla.startswith(gene)]) for gene in set([hla.split("*")[0] for hla in hlas])}
return pd.DataFrame([{"gene": gene, "hlas": ", ".join(gene_list_map[gene])} for gene in sorted(genes)])
itables.init_notebook_mode()
hlas_2024_only = hla_set_2024.difference(hla_set_2026)
hlas_2026_only = hla_set_2026.difference(hla_set_2024)
print("Alleles in 2024 only:")
display(get_hlalist_by_gene(hlas_2024_only))
print("Alleles in 2026 only:")
display(get_hlalist_by_gene(hlas_2026_only))
print(f"2026-only 'alleles' that are P groups: {sum([hla.endswith('P') for hla in hlas_2026_only])} of {len(hlas_2026_only)}")Alleles in 2024 only:
Loading ITables v2.9.1 from the init_notebook_mode cell...
(need help?)
|
| ⓘgene | hlas |
|---|---|
| A | A*02:06, A*25:01, A*26:01, A*32:01, A*34:01 |
| B | B*13:02, B*14:02, B*35:02, B*35:03, B*39:01, B*42:01, B*45:01, B*52:01, B*57:03, B*58:02 |
| C | C*04:01, C*07:01 |
| DPA1 | DPA1*02:01+DPB1*11:01, DPA1*02:12+DPB1*85:01 |
| DQA1 | DQA1*01:03+DQB1*06:02, DQA1*01:03+DQB1*06:09, DQA1*02:01+DQB1*02:02, DQA1*03:03+DQB1*03:01 |
| DRB1 | DRB1*14:54 |
| DRB4 | DRB4*01:03 |
| DRB5 | DRB5*02:02 |
Alleles in 2026 only:
Loading ITables v2.9.1 from the init_notebook_mode cell...
(need help?)
|
| ⓘgene | hlas |
|---|---|
| A | A*23:01P, A*31:01 |
| B | B*40:02, B*49:01 |
| DPA1 | DPA1*01:03+DPB1*04:02, DPA1*01:03+DPB1*104:01, DPA1*02:01+DPB1*131:01 |
| DQA1 | DQA1*01:01+DQB1*06:02, DQA1*01:02+DQB1*05:02, DQA1*01:02+DQB1*06:04, DQA1*01:02+DQB1*06:09, DQA1*01:03+DQB1*06:01, DQA1*04:01+DQB1*04:02, DQA1*05:09+DQB1*03:01 |
| DRB1 | DRB1*01:03, DRB1*04:03, DRB1*04:05, DRB1*04:07, DRB1*08:01, DRB1*08:03, DRB1*08:04, DRB1*11:01, DRB1*11:03, DRB1*14:01P |
| DRB3 | DRB3*01:01P |
| DRB4 | DRB4*01:01P |
2026-only 'alleles' that are P groups: 4 of 26In the 2026 version, we combined several groups of HLA alleles in the same p group (alleles with the same binding domains) into a single model for each p group, after observing that they tended to have the same associated TCRs.
P groups gave us more statistical power to model alleles in 2026, but in the CMV ECOcluster they only replaced two 2024 alleles:
DRB1*14:54 -> DRB1*14:01PDRB4*01:03 -> DRB4*01:01Ppdf_merged_tcr_inner = pdf_merged_tcr_outer[pdf_merged_tcr_outer["source"] == "shared"]
pdf_merged_tcr_inner["hla_agree"] = pdf_merged_tcr_inner["hla_2024"] == pdf_merged_tcr_inner["hla_2026"]
proportion_disagree = (pdf_merged_tcr_inner.hla_agree == False).sum() / len(pdf_merged_tcr_inner)
print(f"{proportion_disagree*100:.1f}% of TCR HLA associations disagree between versions")
pdf_merged_tcr_inner.hla_agree.value_counts()22.9% of TCR HLA associations disagree between versionsThat’s a lot of disagreement! Are there many shifts between HLA genes?
pdf_merged_tcr_inner["hla_gene_agree"] = pdf_merged_tcr_inner["hla_gene_2024"] == pdf_merged_tcr_inner["hla_gene_2026"]
proportion_disagree = (pdf_merged_tcr_inner.hla_gene_agree == False).sum() / len(pdf_merged_tcr_inner)
print(f"{proportion_disagree*100:.1f}% of TCR HLA gene associations disagree between versions")
pdf_merged_tcr_inner.hla_gene_agree.value_counts()14.1% of TCR HLA gene associations disagree between versionsThat’s still a lot of disagreement at the HLA gene level! Let’s visualize it.
from cmvividly.plots.sankey import plot_sankey
plot_sankey(pdf_merged_tcr_inner, source_col="hla_gene_2024",
target_col="hla_gene_2026")
passSo, most of that movement is between DQA1 and DRB1, in both directions. That could be due to more LD between some of those alleles. Let’s get more granular and look at individual gene pairs.
pdf_tcrs_disagree_hla = pdf_merged_tcr_inner[pdf_merged_tcr_inner.hla_agree == False]
pdf_disagreeing_allelepair_counts = pdf_tcrs_disagree_hla[["hla_2024", "hla_2026"]].value_counts().reset_index()
f, ax = plt.subplots()
sns.histplot(data=pdf_disagreeing_allelepair_counts, x="count", ax=ax
)
ax.set_xlabel("Number of TCRs")
ax.set_ylabel("Frequency")
ax.set_title("Distribution of disagreeing HLA allele pair counts")
plt.tight_layout()
print("Top pairs of disagreeing HLA alleles:")
pdf_disagreeing_allelepair_counts[:8]Top pairs of disagreeing HLA alleles:
Loading ITables v2.9.1 from the init_notebook_mode cell...
(need help?)
|
| ⓘhla_2024 | hla_2026 | count |
|---|---|---|
| DPA1*01:03+DPB1*03:01 | DPA1*01:03+DPB1*104:01 | 579 |
| DRB1*03:01 | DQA1*05:01+DQB1*02:01 | 383 |
| DRB1*14:54 | DRB1*14:01P | 220 |
| DQA1*01:04+DQB1*05:03 | DRB1*14:01P | 190 |
| DQA1*02:01+DQB1*02:02 | DRB1*07:01 | 149 |
| DRB4*01:03 | DRB4*01:01P | 135 |
| DQA1*01:02+DQB1*06:02 | DRB5*01:01 | 106 |
| DQA1*01:03+DQB1*06:03 | DRB1*13:01 | 103 |
So, a handful of outlier allele pairs account for most of the allele association disagreement.
DRB1*14:54 -> DRB1*14:01P and DRB4*01:03 -> DRB4*01:01P, as I noted above.DPB1*03:01 -> DPB1*104:01
import numpy as np
pdf_hla_count_2024 = pdf_cmv_ecocluster_2024_dedup.groupby(["hla", "hla_class"]).size().reset_index()
pdf_hla_count_2024.rename(columns={0: "tcrs_2024"}, inplace=True)
pdf_hla_count_2026 = pdf_cmv_ecocluster_2026.groupby(["hla", "hla_class"]).size().reset_index()
pdf_hla_count_2026.rename(columns={0: "tcrs_2026"}, inplace=True)
pdf_hla_counts_bothversions = pdf_hla_count_2024.merge(
pdf_hla_count_2026,
on=["hla_class", "hla"],
how="outer",
suffixes=("_2024", "_2026")
).fillna(0)
pdf_hla_counts_bothversions["log10_tcrs_2024"] = np.log10(pdf_hla_counts_bothversions["tcrs_2024"] + 1)
pdf_hla_counts_bothversions["log10_tcrs_2026"] = np.log10(pdf_hla_counts_bothversions["tcrs_2026"] + 1)
f, ax = plt.subplots()
sns.scatterplot(data=pdf_hla_counts_bothversions, x="log10_tcrs_2024", y="log10_tcrs_2026", ax=ax,
hue="hla_class", s=10)
# add abline
sns.lineplot(x=np.linspace(0, 3), y=np.linspace(0, 3), color="black", linewidth=1, linestyle="--", alpha=0.5, ax=ax)
ax.set_xlabel("${log_{10}}$ 2024 count")
ax.set_ylabel("${log_{10}}$ 2026 count")
ax.set_title("TCR counts per HLA allele in 2026 vs 2024\n(log scale, missing = 0.0, line is 1:1)")
plt.tight_layout()
f.set_size_inches(7, 6)
Takeaways:
Those top alleles in 2026 are:
pdf_hla_counts_bothversions.sort_values("tcrs_2026", ascending=False)[:8]
Loading ITables v2.9.1 from the init_notebook_mode cell...
(need help?)
|
| ⓘ | hla | hla_class | tcrs_2024 | tcrs_2026 | log10_tcrs_2024 | log10_tcrs_2026 |
|---|---|---|---|---|---|---|
| 83 | DRB1*07:01 | cii | 871.0 | 6562.0 | 2.940516 | 3.817102 |
| 73 | DRB1*01:01 | cii | 606.0 | 3565.0 | 2.783189 | 3.552181 |
| 105 | DRB5*01:01 | cii | 151.0 | 2929.0 | 2.181844 | 3.466868 |
| 77 | DRB1*04:01 | cii | 778.0 | 2918.0 | 2.891537 | 3.465234 |
| 70 | DQA1*05:01+DQB1*02:01 | cii | 151.0 | 2191.0 | 2.181844 | 3.340841 |
| 98 | DRB1*15:01 | cii | 221.0 | 1897.0 | 2.346353 | 3.278296 |
| 89 | DRB1*11:01 | cii | 0.0 | 1832.0 | 0.000000 | 3.263162 |
| 57 | DQA1*01:02+DQB1*06:02 | cii | 386.0 | 1595.0 | 2.587711 | 3.203033 |
We’ve summarized the biggest differences between the two versions of the CMV ECOcluster. But given that the 2024 ECOcluster has been available for two years and has been independently established as a sensitive, specific biomarker of CMV infection, it’s concerning that half those TCRs just go away in the new version. To put a point on it: the old one was good, the new one is good… but is some kind of hybrid better (or at least more comprehensive) than either?
I won’t answer that question completely in this post. In this post, I’ll just catalog what we’re “losing” in the 2026 ECOcluster and speculate about why. I plan to dig further in another post, using CMV-labeled repertoires.
Let’s look at the HLA alleles that lost the most TCRs. Here, I don’t mean the TCR shows up in the 2026 ECOcluster associated with a different allele. I mean actually lost. We already learned above that those are about 50/50 Class I vs. Class II. How does that break down by allele?
pdf_cmv_ecocluster_2024_dedup["in_2026"] = pdf_cmv_ecocluster_2024_dedup["tcr"].isin(tcr_set_2026)
pdf_for_agg = pdf_cmv_ecocluster_2024_dedup.copy()
pdf_for_agg["not_in_2026"] = ~pdf_for_agg.in_2026
pdf_for_agg["n_lost"] = pdf_for_agg["not_in_2026"]
pdf_for_agg["n_original"] = 1
pdf_loss_summary = pdf_for_agg.groupby(["hla", "hla_class", "hla_gene"]).agg({
"in_2026": "mean",
"n_original": "sum",
"n_lost": "sum"
}).reset_index()
pdf_loss_summary["proportion_lost"] = pdf_loss_summary["n_lost"] / pdf_loss_summary.n_original
pdf_loss_summary.sort_values("n_lost", ascending=False)[:10]
Loading ITables v2.9.1 from the init_notebook_mode cell...
(need help?)
|
| ⓘ | hla | hla_class | hla_gene | in_2026 | n_original | n_lost | proportion_lost |
|---|---|---|---|---|---|---|---|
| 39 | DPA1*01:03+DPB1*03:01 | cii | DPA1 | 0.508317 | 1503 | 739 | 0.491683 |
| 24 | B*42:01 | ci | B | 0.000000 | 599 | 599 | 1.000000 |
| 38 | DPA1*01:03+DPB1*02:01 | cii | DPA1 | 0.520513 | 1170 | 561 | 0.479487 |
| 80 | DRB5*02:02 | cii | DRB5 | 0.091727 | 556 | 505 | 0.908273 |
| 18 | B*35:01 | ci | B | 0.269058 | 669 | 489 | 0.730942 |
| 15 | B*15:01 | ci | B | 0.370803 | 685 | 431 | 0.629197 |
| 37 | C*08:02 | ci | C | 0.150743 | 471 | 400 | 0.849257 |
| 21 | B*38:01 | ci | B | 0.061321 | 424 | 398 | 0.938679 |
| 33 | C*04:01 | ci | C | 0.007792 | 385 | 382 | 0.992208 |
| 64 | DRB1*04:04 | cii | DRB1 | 0.025140 | 358 | 349 | 0.974860 |
n_alleles_totally_lost = sum(pdf_loss_summary.proportion_lost == 1)
print(f"We lost all the TCRs from {n_alleles_totally_lost} alleles.")We lost all the TCRs from 18 alleles.Let’s look at the ones that had the most TCRs in 2024.
pdf_loss_summary.sort_values(["proportion_lost", "n_lost"], ascending=False)
Loading ITables v2.9.1 from the init_notebook_mode cell...
(need help?)
|
| ⓘ | hla | hla_class | hla_gene | in_2026 | n_original | n_lost | proportion_lost |
|---|---|---|---|---|---|---|---|
| 24 | B*42:01 | ci | B | 0.000000 | 599 | 599 | 1.000000 |
| 8 | A*32:01 | ci | A | 0.000000 | 250 | 250 | 1.000000 |
| 20 | B*35:03 | ci | B | 0.000000 | 192 | 192 | 1.000000 |
| 22 | B*39:01 | ci | B | 0.000000 | 160 | 160 | 1.000000 |
| 13 | B*13:02 | ci | B | 0.000000 | 149 | 149 | 1.000000 |
| 29 | B*52:01 | ci | B | 0.000000 | 143 | 143 | 1.000000 |
| 6 | A*25:01 | ci | A | 0.000000 | 136 | 136 | 1.000000 |
| 19 | B*35:02 | ci | B | 0.000000 | 130 | 130 | 1.000000 |
| 2 | A*02:06 | ci | A | 0.000000 | 115 | 115 | 1.000000 |
| 7 | A*26:01 | ci | A | 0.000000 | 115 | 115 | 1.000000 |
| (71 more rows not shown) | |||||||
Plot the proportion TCRs lost vs. the original TCR count.
f, ax = plt.subplots()
sns.scatterplot(x="n_original", y="proportion_lost", data=pdf_loss_summary,
hue="hla_class")
ax.set_title("Proportion TCRs 'lost' per allele vs. original count")
So, the biggest losers, in terms of proportion of TCRs lost, seem to be Class I alleles. That’s as expected, given that the 2026 ECOcluster is far more Class II-dominated and that, even though it’s much bigger overall, it actually has many fewer Class I-associated TCRs.
Highly sequence-similar TCRs tend to bind the same antigen. So, if many of those “lost” 2024 TCRs are sequence-similar to a 2026 TCR (I’ll define that here as same V gene, same J gene, same CDR3 length, one CDR3 edit), that’s strong evidence they shouldn’t have been lost. 1-Hamming doesn’t always mean same-antigen binding, but the additional evidence that it was CMV-associated in the 2024 version makes it pretty compelling.
from cmvividly.data.hamming1_pairs import find_hamming1_pairs_same_vj
# find all Hamming-1 pairs with the same v and j gene in the outer join between 2024 and 2026.
# That's overkill, but that's the tooling I've built so far.
pdf_ham1_pairs_all = find_hamming1_pairs_same_vj(
pdf_merged_tcr_outer)
pdf_ham1_pairs_all = (
pdf_ham1_pairs_all
.drop(columns=["row_i", "row_j"])
.rename(columns={
"tcr_i": "tcr_2024",
"tcr_j": "tcr_2026"
})
)
# filter to just the ones present in 2024...
pdf_ham1_pairs_in2024 = pdf_ham1_pairs_all[pdf_ham1_pairs_all.tcr_2024.isin(set(pdf_cmv_ecocluster_2024.tcr))]
# ...with a match in 2026.
pdf_ham1_pairs_in2024_pairin2026 = pdf_ham1_pairs_in2024[pdf_ham1_pairs_in2024.tcr_2026.isin(set(pdf_cmv_ecocluster_2026.tcr))]
tcrs_2024_ham1_match_2026 = set(pdf_ham1_pairs_in2024_pairin2026.tcr_2024)
lost_tcrs_2024 = set(pdf_cmv_ecocluster_2024_dedup[~pdf_cmv_ecocluster_2024_dedup.in_2026].tcr)
lost_tcrs_2024_ham1_match_2026 = tcrs_2024_ham1_match_2026.intersection(lost_tcrs_2024)
n_lost_2024 = len(lost_tcrs_2024)
n_lost_with2026match = len(lost_tcrs_2024_ham1_match_2026)
proportion_oflost_withmatch = n_lost_with2026match / n_lost_2024
print(f"{n_lost_with2026match} of {n_lost_2024} 'lost' TCRs ({proportion_oflost_withmatch*100:.1f}%) have a Hamming-1 match in 2026")2132 of 13577 'lost' TCRs (15.7%) have a Hamming-1 match in 2026A lot of those “lost” TCRs have a match! Do they have just one? Lots of matches would be more evidence that they shouldn’t have been lost.
pdf_losttcr_matchcounts = pdf_ham1_pairs_in2024_pairin2026.tcr_2024.value_counts().reset_index()
f, ax = plt.subplots()
sns.histplot(pdf_losttcr_matchcounts["count"])
f.set_size_inches(8, 3)
ax.set_title("Hamming-1 matches in 2026 per matched 'lost' 2024 TCR")
Almost half of the “lost” 2024 TCRs have at least 1 Hamming-1 match in 2026… and more than a quarter have 2+ matches… and some have as many as 20 Hamming-1 matches! That’s very strong evidence that some of those TCRs, at least, shouldn’t have been lost.
Why were they lost, then? Some ideas:
In any case, this analysis is a strong suggestion that many of those “lost” 2024 TCRs are actually CMV-specific.
The 2026 CMV ECOcluster is bigger and possibly “better”, overall. But the 2024 version has many TCRs and HLA alleles missing from 2026. If you’re applying the 2026 ECOcluster to repertoires, it’s worth trying the same analysis with the 2024 ECOcluster, and also combining them… carefully.
At the very least, we need to deal with the combining of groups of HLA alleles into p-group models that was done in the 2026 version but not in 2024. That allele mapping is in a big table in the supplementary materials for the 2026. The practical impact for the CMV ECOcluster is as described above:
DRB1*14:54 -> DRB1*14:01PDRB4*01:03 -> DRB4*01:01PIn the code block below, I’ll take a quick stab at combining the two and dealing with those two p groups. This isn’t the final word on the best way to combine the two ECOclusters, so I won’t put this in library code yet.
def combine_2024_2026_cmv_ecoclusters(pdf_cmv_ecocluster_2024, pdf_cmv_ecocluster_2026):
"""This is a quick and dirty stab at combining the two ECOclusters.
I'l button it up later, after seeing how the 2024 and 2026 TCRs behave in repertoires"""
pdf_cmv_ecocluster_2024_formerge = pdf_cmv_ecocluster_2024_dedup.copy()
pdf_cmv_ecocluster_2024_formerge["hla"] = (
pdf_cmv_ecocluster_2024_formerge.hla
.replace("DRB1*14:54", "DRB1*14:01P")
.replace("DRB4*01:03", "DRB4*01:01P")
)
tcr_columns = ["cdr3", "vgene", "jgene"]
pdf_combined_ecocluster = pdf_cmv_ecocluster_2024_formerge.merge(
pdf_cmv_ecocluster_2026, on=tcr_columns, how="outer",
suffixes=["_2024", "_2026"])
pdf_combined_ecocluster["hla"] = [
row["hla_2026"] if isinstance(row["hla_2026"], str)
else row["hla_2024"]
for _, row in pdf_combined_ecocluster.iterrows()
]
return pdf_combined_ecocluster
pdf_combined_ecocluster = combine_2024_2026_cmv_ecoclusters(pdf_cmv_ecocluster_2024, pdf_cmv_ecocluster_2026)Let’s see how many TCRs and HLA alleles that combined ECOcluster comprises.
n_tcrs = len(pdf_combined_ecocluster)
n_hlas = len(set(pdf_combined_ecocluster.hla))
print(f"Combined ECOcluster has {n_tcrs} TCRs associated with {n_hlas} HLA alleles")Combined ECOcluster has 66024 TCRs associated with 105 HLA allelesIn a future post, I plan to explore the ideas that the 2024-only TCRs are CMV-specific, and that the best way to use the two ECOclusters is to combine them, by seeing how those TCRs behave in CMV-labeled repertoires. I’ll also take a more formal stab at combining them in a way that deals with the differences in how they were built.