Exhausted CD8 T cells are associated with ICB responsiveness in breast cancer

Background

The single-cell dataset generated by Bassez et al. (2021) Nature Medicine contains single-cell data from breast carcinoma biopsies from 29 patients taken immediately before (pre-treatment) and 9 +/- 2 days after (on-treatment) receiving anti-PD1 immunotherapy. Lacking an objective measure of clinical benefit following therapy, in this study the authors use T cell expansion following treatment as a proxy for response.

Full Data Original Source

Case Study Dataset Internal link

Goal

To compare the T cell subtype composition between responders and non-responders using ProjecTILs and reference maps for CD8 T cells and CD4 T cells.

R Environment

library(renv)
renv::restore()

library(patchwork)
library(ggplot2)
library(reshape2)
library(Seurat)
library(ProjecTILs)

scRNA-seq data preparation

options(timeout = 5000)
ddir <- "input/Bassez_breast_data"

if (!file.exists(ddir)) {
    dir.create(ddir)
    dataUrl <- "https://figshare.com/ndownloader/files/47900725"
    download.file(dataUrl, paste0(ddir, "/tmp.zip"))
    unzip(paste0(ddir, "/tmp.zip"), exdir = ddir)
    file.remove(paste0(ddir, "/tmp.zip"))
}
# Count matrices
f1 <- sprintf("%s/1864-counts_tcell_cohort1.rds", ddir)

cohort1 <- readRDS(f1)
dim(cohort1)

[1] 25288 53382

# Metadata
meta1 <- read.csv(sprintf("%s/1870-BIOKEY_metaData_tcells_cohort1_web.csv", ddir))
rownames(meta1) <- meta1$Cell

head(meta1)

                                                             Cell nCount_RNA
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1 BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1       1252
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1 BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1       1869
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1 BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1       1000
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1 BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1       1288
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1 BIOKEY_13_Pre_AAAGATGGTTTGACAC-1       2056
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1 BIOKEY_13_Pre_AAAGATGTCAACGCTA-1       1224
                                 nFeature_RNA patient_id timepoint expansion
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1          700  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1          995  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1          627  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1          681  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1          789  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1          634  BIOKEY_13       Pre       n/a
                                 BC_type cellSubType          cohort
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1   HER2+         gdT treatment_naive
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1   HER2+     NK_REST treatment_naive
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1   HER2+       CD4_N treatment_naive
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1   HER2+      CD8_EM treatment_naive
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1   HER2+       CD8_N treatment_naive
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1   HER2+      CD8_RM treatment_naive

data.seurat <- CreateSeuratObject(cohort1, project = "Cohort1_IT", meta.data = meta1)
data.seurat <- NormalizeData(data.seurat)

Subset pre-treatment biopsies

table(data.seurat$timepoint)


   On   Pre 
27528 25854

data.seurat <- subset(data.seurat, subset = timepoint == "Pre")

Subset T cells according to annotation by the authors

table(data.seurat$cellSubType)


               CD4_EM                CD4_EX  CD4_EX_Proliferating 
                 4810                  1707                   107 
                CD4_N               CD4_REG CD4_REG_Proliferating 
                 3007                  2861                    69 
               CD8_EM              CD8_EMRA                CD8_EX 
                 5443                   290                  2286 
 CD8_EX_Proliferating                 CD8_N                CD8_RM 
                  340                   354                  1963 
                  gdT               NK_CYTO               NK_REST 
                 1000                   230                   938 
           Vg9Vd2_gdT 
                  449

data.seurat <- subset(data.seurat, subset = cellSubType %in% c("NK_REST", "Vg9Vd2_gdT",
    "gdT", "NK_CYTO"), invert = T)

We will remove samples that are too small and downsample large ones, in order to have similar contribution from all patients/samples. Downsampling large samples also speeds up downstream computations.

ds <- 1000
min.cells <- 200

tab <- table(data.seurat$patient_id)
keep <- names(tab)[tab > min.cells]
data.seurat <- subset(data.seurat, patient_id %in% keep)


Idents(data.seurat) <- "patient_id"
data.seurat <- subset(data.seurat, cells = WhichCells(data.seurat, downsample = ds))
table(data.seurat$patient_id)


 BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16 
     1000      1000       555      1000      1000       939      1000      1000 
BIOKEY_19  BIOKEY_2 BIOKEY_21 BIOKEY_24 BIOKEY_27 BIOKEY_28  BIOKEY_3 BIOKEY_31 
     1000       960       203       347       601       594       266       349 
 BIOKEY_4  BIOKEY_5  BIOKEY_6 
     1000      1000       616

ProjecTILs analysis

Load reference maps

Here we will project the query data onto human reference TIL atlas for CD4 and CD8 T cells, with the goal to interpret T cell diversity in the context of reference T cell subtypes. First, load the reference maps:

download.file("https://figshare.com/ndownloader/files/41414556", destfile = "CD8T_human_ref_v1.rds")
ref.cd8 <- load.reference.map("CD8T_human_ref_v1.rds")

[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map Human CD8 TILs"

download.file("https://figshare.com/ndownloader/files/39012395", destfile = "CD4T_human_ref_v1.rds")
ref.cd4 <- load.reference.map("CD4T_human_ref_v1.rds")

[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map custom_reference"

a <- DimPlot(ref.cd8, cols = ref.cd8@misc$atlas.palette, label = T) + theme(aspect.ratio = 1) +
    ggtitle("CD8 T reference") + NoLegend()

b <- DimPlot(ref.cd4, cols = ref.cd4@misc$atlas.palette, label = T) + theme(aspect.ratio = 1) +
    ggtitle("CD4 T reference") + NoLegend()

a | b

Show key marker genes for the reference subtypes

DefaultAssay(ref.cd8) <- "RNA"
Idents(ref.cd8) <- "functional.cluster"

genes <- c("SELL", "TCF7", "LEF1", "CCR7", "S1PR1", "LMNA", "IL7R", "GZMK", "FGFBP2",
    "FCGR3A", "XCL1", "XCL2", "CD200", "CRTAM", "GNG4", "TOX", "PDCD1", "HAVCR2",
    "GZMB", "PRF1", "LAG3", "KLRB1", "TRAV1-2")

DotPlot(ref.cd8, features = genes, cols = "RdBu", scale = T, col.max = 1.5) + theme(axis.text.x = element_text(angle = 45,
    hjust = 1)) + ggtitle("CD8 T cell reference map")

DefaultAssay(ref.cd4) <- "RNA"
Idents(ref.cd4) <- "functional.cluster"

genes <- c("SELL", "TCF7", "S1PR1", "KLF2", "IL7R", "CCL5", "GZMK", "EOMES", "TBX21",
    "CXCR6", "GZMA", "GZMB", "PRF1", "GNLY", "CXCL13", "PDCD1", "TOX", "LAG3", "HAVCR2",
    "KLRB1", "IL17A", "FOXP3", "IL2RA")

DotPlot(ref.cd4, features = genes, cols = "RdBu", scale = T, col.max = 1.5) + theme(axis.text.x = element_text(angle = 45,
    hjust = 1)) + ggtitle("CD4 T cell reference map")

Reference-based analysis

To classify CD4 and CD8 T cell subtypes, we will project CD4 and CD8 T cells separately into the two reference maps. Running ProjecTILs.classifier sequentially on the two maps allows combining the subtype predictions for a complete annotation of T cells subtypes.

For optimal batch-effect correction, we recommend projecting each patient/batch separately (set split.by="patient_id"). For large samples, Run.ProjecTILs can take several minutes. Set the number of parallel jobs with ncores according to your computational resources.

# Classify CD8 T subtypes
ncores = 8
DefaultAssay(ref.cd8) <- "integrated"
data.seurat <- ProjecTILs.classifier(data.seurat, ref.cd8, ncores = ncores, split.by = "patient_id")

After running the classifier with the CD8T reference, we obtain CD8 T cell subtype annotations:

table(data.seurat$functional.cluster, useNA = "ifany")


       CD8.CM        CD8.EM      CD8.MAIT CD8.NaiveLike     CD8.TEMRA 
         1557          1859            16           388           102 
      CD8.TEX      CD8.TPEX          <NA> 
         1169            49          9290

We can now run the classifier on the same object with the CD4 T reference, specifying the parameter overwrite=FALSE. This adds labels for the CD4 T cells without replacing the CD8 T cell subtypes. Note that only cells that can be unequivocally assigned to one cell type are labeled; the remaining cells (including those that receive a label by more than one reference) are assigned to NA.

# Classify CD4 T subtypes
DefaultAssay(ref.cd4) <- "integrated"
data.seurat <- ProjecTILs.classifier(data.seurat, ref.cd4, ncores = ncores, split.by = "patient_id",
    overwrite = FALSE)

The labels now contain both CD4T and CD8T subtypes:

table(data.seurat$functional.cluster, useNA = "ifany")


CD4.CTL_EOMES   CD4.CTL_Exh  CD4.CTL_GNLY CD4.NaiveLike       CD4.Tfh 
          910           472           106          4318           628 
     CD4.Th17      CD4.Treg        CD8.CM        CD8.EM      CD8.MAIT 
           24          1522          1554          1848            16 
CD8.NaiveLike     CD8.TEMRA       CD8.TEX      CD8.TPEX          <NA> 
          387           100          1166            49          1330

We can verify the expression profile using a panel of marker genes.

genes4radar <- c("CD4", "CD8A", "TCF7", "CCR7", "IL7R", "LMNA", "GZMA", "GZMK", "FGFBP2",
    "XCL1", "CD200", "CRTAM", "TOX", "PDCD1", "HAVCR2", "PRF1", "GLNY", "GZMB", "TRAV1-2",
    "KLRB1", "FOXP3")

plot.states.radar(ref = ref.cd8, query = data.seurat, genes4radar = genes4radar,
    min.cells = 20)

plot.states.radar(ref = ref.cd4, query = data.seurat, genes4radar = genes4radar,
    min.cells = 20)

Compare subtype composition of responders vs. non-responders (based on clonal expansion) on pre-therapy samples

a <- ref.cd8@misc$atlas.palette
b <- ref.cd4@misc$atlas.palette
cols <- c(a, b)

by.exp <- SplitObject(data.seurat, split.by = "expansion")

tb <- table(factor(by.exp$E$functional.cluster, levels = names(cols)))
tb <- tb * 100/sum(tb)
tb.m <- reshape2::melt(tb)
colnames(tb.m) <- c("Cell_state", "Perc_cells")

p1 <- ggplot(tb.m, aes(x = Cell_state, y = Perc_cells, fill = Cell_state)) + geom_bar(stat = "identity") +
    theme_bw() + scale_fill_manual(values = cols) + theme(axis.text.x = element_blank(),
    legend.position = "none") + ggtitle("Responders") + ylim(0, 60)


tb <- table(factor(by.exp$NE$functional.cluster, levels = names(cols)))
tb <- tb * 100/sum(tb)
tb.m <- reshape2::melt(tb)
colnames(tb.m) <- c("Cell_state", "Perc_cells")

p2 <- ggplot(tb.m, aes(x = Cell_state, y = Perc_cells, fill = Cell_state)) + geom_bar(stat = "identity") +
    theme_bw() + scale_fill_manual(values = cols) + theme(axis.text.x = element_blank(),
    legend.position = "right") + ggtitle("Non-Responders") + ylim(0, 60)

p1 | p2

The analysis shows that ICB-responding tumors have a higher proportion of CD8 exhausted (CD8.TEX), CD8 precursor exhausted (CD8.TPEX) and CD4 cytotoxic exhausted (CD4.CTL_Exh) T cells; on the other hand, CD4 Naive-like and CD8 TEMRA are more abundant in non-responders.

A more compact way of displaying enrichment/depletion of specific subtypes is to calculate and visualize fold-changes of TIL composition in responders vs non-responders

which.types <- table(data.seurat$functional.cluster) > 50  # disregard small populations
which.types <- names(which.types)[which.types]

states_all <- names(cols)

cols_use <- cols[which.types]

norm.c <- table(by.exp$NE$functional.cluster)/sum(table(by.exp$NE$functional.cluster))
norm.q <- table(by.exp$E$functional.cluster)/sum(table(by.exp$E$functional.cluster))

foldchange <- norm.q[which.types]/norm.c[which.types]
foldchange <- sort(foldchange, decreasing = T)

tb.m <- reshape2::melt(foldchange)
colnames(tb.m) <- c("Cell_state", "Fold_change")
ggplot(tb.m, aes(x = Cell_state, y = Fold_change, fill = Cell_state)) + geom_bar(stat = "identity") +
    scale_fill_manual(values = cols_use) + theme_bw() + geom_hline(yintercept = 1) +
    scale_y_continuous(trans = "log2") + theme(axis.text.x = element_blank(), legend.position = "left") +
    ggtitle("Responders vs. Non-responders")

Tumors in ICB-responders have a >10 fold-change increase in CD4 exhausted CTL and CD8 exhausted T cells.

Correspondence analysis of subtype composition

Do patients cluster based on their subtype composition? are subtype composition “classes” related to response to ICB? we can perform correspondence analysis (CA) to try to answer these questions, and group patients based on their composition of T cell subtypes.

library("FactoMineR")
min.cells <- 30

# Only use samples with sufficient amount of cells
tab <- table(data.seurat$patient_id)
samples.use <- names(tab)[tab > min.cells]
data.seurat <- subset(data.seurat, subset = patient_id %in% samples.use)

subtype_comp_table <- prop.table(table(data.seurat$functional.cluster, data.seurat$patient_id),
    margin = 2)
res.CA <- CA(t(subtype_comp_table), graph = FALSE)
a <- plot(res.CA) + theme(aspect.ratio = 1)

# Color by response
coords <- as.data.frame(res.CA$row$coord)
coords.ct <- as.data.frame(res.CA$col$coord)
coords.all <- rbind(coords[, 1:5], coords.ct[, 1:5])

is.R <- unique(data.seurat$patient_id[data.seurat$expansion == "E"])
is.NR <- unique(data.seurat$patient_id[data.seurat$expansion == "NE"])
coords[["Responder"]] <- "n/a"
coords$Responder[rownames(coords) %in% is.R] <- "Responder"
coords$Responder[rownames(coords) %in% is.NR] <- "Non-responder"
coords[["Patient"]] <- rownames(coords)

b <- ggplot(coords, aes(x = `Dim 1`, y = `Dim 2`, color = Responder)) + geom_point(size = 4) +
    theme_bw() + theme(aspect.ratio = 1) + scale_color_manual(values = c("#E9E9E9",
    "#fb0505", "#08e579")) + xlim(min(coords.all[, "Dim 1"]), max(coords.all[, "Dim 1"])) +
    ylim(min(coords.all[, "Dim 2"]), max(coords.all[, "Dim 2"]))

a | b

Discussion

Our automated ProjecTILs analysis of this patient cohort indicate that compared to non-responders, early responders to PD-1 blockade (in terms of T cell expansion) have higher relative amounts of exhausted T cells (both CD8+ and CD4+) infiltrating their tumors, in agreement with previous studies (Daud et al. 2016; Thommen et al. 2018; Kumagai et al. 2020; Bassez et al. 2020). This supports the notion that patients with a pre-existing tumor-specific T cell response in their tumors are more likely to respond to PD-1 blockade therapy.

Because tumoricidal CD8 Tex are derived from Tpex (Siddiqui et al. 2019; Miller et al. 2019), the presence of CD8 Tex seems to be a good proxy for the presence of smaller populations of quiescent Tpex, that have the ability to proliferate and differentiate into Tex cells following PD-1 blockade (Siddiqui et al. 2019; Miller et al. 2019). Indeed we observed a small population of Tpex-like cells, with increased levels R vs RN. Because there are very few Tpex their gene profile might be less robust than others, but compared to Tex they display higher levels of Tpex-associated genes, such as XCL1 and CD200, and lower levels of terminal-exhaustion markers such as HAVCR2 and ENTPD1.

References

Bassez, A., Vos, H., Van Dyck, L. et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 27, 820–832 (2021). https://doi.org/10.1038/s41591-021-01323-8
Daud, A.I., Loo, K., Pauli, M.L., Sanchez-Rodriguez, R., Sandoval, P.M., Taravati, K., Tsai, K., Nosrati, A., Nardo, L., Alvarado, M.D., Algazi, A.P., Pampaloni, M.H., Lobach, I.V., Hwang, J., Pierce, R.H., Gratz, I.K., Krummel, M.F., Rosenblum, M.D., 2016. Tumor immune profiling predicts response to anti-PD-1 therapy in human melanoma. J. Clin. Invest. 126, 3447–3452. https://doi.org/10.1172/JCI87324
Gros, A., Robbins, P.F., Yao, X., Li, Y.F., Turcotte, S., Tran, E., Wunderlich, J.R., Mixon, A., Farid, S., Dudley, M.E., Hanada, K., Almeida, J.R., Darko, S., Douek, D.C., Yang, J.C., Rosenberg, S. a, 2014. PD-1 identifies the patient-specific in filtrating human tumors. J. Clin. Invest. 124, 2246–59. https://doi.org/10.1172/JCI73639.2246
Miller, B.C., Sen, D.R., Al Abosy, R., Bi, K., Virkud, Y.V., LaFleur, M.W., Yates, K.B., Lako, A., Felt, K., Naik, G.S., Manos, M., Gjini, E., Kuchroo, J.R., Ishizuka, J.J., Collier, J.L., Griffin, G.K., Maleri, S., Comstock, D.E., Weiss, S.A., Brown, F.D., Panda, A., Zimmer, M.D., Manguso, R.T., Hodi, F.S., Rodig, S.J., Sharpe, A.H., Haining, W.N., 2019. Subsets of exhausted CD8+ T cells differentially mediate tumor control and respond to checkpoint blockade. Nat. Immunol. 20, 326–336. https://doi.org/10.1038/s41590-019-0312-6
Siddiqui, I., Schaeuble, K., Chennupati, V., Fuertes Marraco, S.A., Calderon-Copete, S., Pais Ferreira, D., Carmona, S.J., Scarpellino, L., Gfeller, D., Pradervand, S., Luther, S.A., Speiser, D.E., Held, W., 2019. Intratumoral Tcf1+PD-1+CD8+ T Cells with Stem-like Properties Promote Tumor Control in Response to Vaccination and Checkpoint Blockade Immunotherapy. Immunity 50, 195–211.e10. https://doi.org/10.1016/J.IMMUNI.2018.12.021
Thommen, D.S., Koelzer, V.H., Herzig, P., Roller, A., Trefny, M., Dimeloe, S., Kiialainen, A., Hanhart, J., Schill, C., Hess, C., Prince, S.S., Wiese, M., Lardinois, D., Ho, P.C., Klein, C., Karanikas, V., Mertz, K.D., Schumacher, T.N., Zippelius, A., 2018. A transcriptionally and functionally distinct pd-1 + cd8 + t cell pool with predictive potential in non-small-cell lung cancer treated with pd-1 blockade. Nat. Med. 24, 994. https://doi.org/10.1038/s41591-018-0057-z
Kumagai, S., Togashi, Y., Kamada, T. et al. The PD-1 expression balance between effector and regulatory T cells predicts the clinical efficacy of PD-1 blockade therapies. Nat Immunol 21, 1346–1358 (2020). https://doi.org/10.1038/s41590-020-0769-3

Exhausted CD8 T cells are associated with ICB responsiveness in breast cancer