Effect of anti-PD-L1 therapy on CD4+ T cells in chronic infection

Background

The PD1:PD-L1 immune checkpoint is a crucial target for immunotherapies. However, the effect of PD-L1 blockade on CD4+ T cells is not well understood. Snell et al. (2021) set out to characterize CD4+ T cell heterogeneity during chronic viral infections, and how these cells were affected by anti-PD-L1 therapy. They found that anti-PD-L1 blockade specifically increased the proportion of Th1 cells, and restored their cytotoxic capacities.

In their study, the also generated scRNA-seq data of virus-specific CD4+ T cells isolated from mice chronically infected with LMCV (33 days p.i.), after treatment with an anti-PD-L1 antibody or with an isotype antibody (control). We will apply ProjecTILs to re-analyse these samples in the context of a reference atlas of CD4+ T cells in viral infection. This atlas was described by Andreatta et al. (2022) eLife and can also be explored interactively through the SPICA portal at: https://spica.unil.ch/refs/viral-CD4-T

R Environment

Load the required packages. If you clone the repository, you may reconstitute the R environment using renv::restore()

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

Download scRNA-seq data

We can download the single-cell data for this study directly from Gene Expression Omnibus (GEO).

Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS = "true")
options(timeout = max(300, getOption("timeout")))

library(GEOquery)
geo_acc <- "GSE163345"
datadir <- "input/Snell"

gse <- getGEO(geo_acc)

system(paste0("mkdir -p ", datadir))
getGEOSuppFiles(geo_acc, baseDir = datadir)

Read in the count matrices in 10x format

system(sprintf("tar -xvf %s/%s/GSE163345_RAW.tar -C %s", datadir, geo_acc, datadir))

ids <- c("GSM4977292", "GSM4977293")

iso <- ReadMtx(mtx = sprintf("%s/%s_CD4_ISO_matrix.mtx.gz", datadir, ids[1]), cells = sprintf("%s/%s_CD4_ISO_barcodes.tsv.gz",
    datadir, ids[1]), features = sprintf("%s/%s_CD4_ISO_genes.tsv.gz", datadir, ids[1]))

pdl <- ReadMtx(mtx = sprintf("%s/%s_CD4_PDL_matrix.mtx.gz", datadir, ids[2]), cells = sprintf("%s/%s_CD4_PDL_barcodes.tsv.gz",
    datadir, ids[2]), features = sprintf("%s/%s_CD4_PDL_genes.tsv.gz", datadir, ids[2]))

scRNA-seq preprocessing and QC

Create Seurat objects

# Unique barcodes
colnames(iso) <- paste("Sne1", colnames(iso), sep = "_")
colnames(pdl) <- paste("Sne2", colnames(pdl), sep = "_")

data_iso <- CreateSeuratObject(counts = iso, assay = "RNA")
data_iso$Condition <- "isotype"

data_pdl <- CreateSeuratObject(counts = pdl, assay = "RNA")
data_pdl$Condition <- "PDL1_treated"

# Merge in single object
data_seurat <- merge(data_iso, data_pdl)
table(data_seurat$Condition)


     isotype PDL1_treated 
         965         1743

Basic quality control

data_seurat <- NormalizeData(data_seurat)

patterns <- c("^Rp", "^mt-")
data_seurat <- AddMetaData(data_seurat, metadata = PercentageFeatureSet(data_seurat,
    pattern = patterns[1]), col.name = "percent.ribo")
data_seurat <- AddMetaData(data_seurat, metadata = PercentageFeatureSet(data_seurat,
    pattern = patterns[2]), col.name = "percent.mito")

Idents(data_seurat) <- "Condition"
VlnPlot(data_seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.ribo", "percent.mito"),
    ncol = 2, pt.size = 0.01)

cutoffs <- list()
cutoffs[["percent.ribo"]] <- c(min = 0, max = 50)
cutoffs[["percent.mito"]] <- c(min = 0, max = 6)
cutoffs[["nFeature_RNA"]] <- c(min = 500, max = 3000)
cutoffs[["nCount_RNA"]] <- c(min = 1000, max = 10000)

data_seurat <- subset(data_seurat, subset = nFeature_RNA > cutoffs[["nFeature_RNA"]]["min"] &
    nFeature_RNA < cutoffs[["nFeature_RNA"]]["max"] & nCount_RNA > cutoffs[["nCount_RNA"]]["min"] &
    nCount_RNA < cutoffs[["nCount_RNA"]]["max"] & percent.ribo > cutoffs[["percent.ribo"]]["min"] &
    percent.ribo < cutoffs[["percent.ribo"]]["max"] & percent.mito > cutoffs[["percent.mito"]]["min"] &
    percent.mito < cutoffs[["percent.mito"]]["max"])

Calculate low dimensional embeddings for the two samples

sample.list <- SplitObject(data_seurat, split.by = "Condition")
names(sample.list)

[1] "isotype"      "PDL1_treated"

sample.list <- lapply(sample.list, function(x) {
    ScaleData(x) |>
        FindVariableFeatures(nfeatures = 500) |>
        RunPCA(npcs = 15) |>
        RunUMAP(dims = 1:15)
})
DimPlot(sample.list$isotype) | DimPlot(sample.list$PDL1_treated)

Load reference map

The CD4+ T cell map is described in this paper, and can be downloaded from Figshare at: doi.org/10.6084/m9.figshare.16592693.v2

You can also use the commands below to download the atlas directly within R, and load it into memory.

cd4.atlas.file <- "ref_LCMV_CD4_mouse_release_v1.rds"
if (!file.exists(cd4.atlas.file)) {
    dataUrl <- "https://figshare.com/ndownloader/files/31057081"
    download.file(dataUrl, cd4.atlas.file)
}
ref <- load.reference.map(cd4.atlas.file)

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

palette <- ref@misc$atlas.palette

DimPlot(ref, cols = palette) + theme(aspect.ratio = 1)

Run ProjecTILs

ProjecTILs can be used in two modes: 1) as a classifier, to annotate query cells in their original space; 2) as a method to embed the query dataset into the same space of the reference map. We will apply both below.

1) ProjecTILs classifier

The ProjecTILs.classifier function integrates the query with the reference and transfers cell type labels from reference to query. The original low-dimensional spaces of the query are not modified:

sample.class <- list()
sample.class$isotype <- ProjecTILs.classifier(query = sample.list$isotype, ref = ref,
    filter.cells = FALSE)

[1] "Using assay RNA for query"
[1] "Aligning query to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space

Projecting corrected query onto Reference UMAP space

sample.class$PDL1_treated <- ProjecTILs.classifier(query = sample.list$PDL1_treated,
    ref = ref, filter.cells = FALSE)

[1] "Using assay RNA for query"
[1] "Aligning query to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space

Projecting corrected query onto Reference UMAP space

a <- DimPlot(sample.class$isotype, group.by = "functional.cluster", cols = palette) +
    theme(aspect.ratio = 1) + ggtitle("Isotype")
b <- DimPlot(sample.class$PDL1_treated, group.by = "functional.cluster", cols = palette) +
    theme(aspect.ratio = 1) + ggtitle("PDL1-treated")
a | b

We can compare the subtype composition between anti-PD-L1 treated and isotype-treated CD4+ T cells. We see a nearly three-fold increase in Th1 subtypes upon anti-PD-L1 blockade compared to control.

a <- plot.statepred.composition(ref = ref, query = sample.class$isotype, metric = "Percent") +
    ggtitle("Isotype") + ylim(0, 40)
b <- plot.statepred.composition(ref = ref, query = sample.class$PDL1_treated, metric = "Percent") +
    ggtitle("PDL1_treated") + ylim(0, 40)
a | b

2) ProjecTILs reference embedding

The Run.ProjecTILs function integrates the query with the reference, transfers cell type labels from reference to query, and embeds the cells from the query in the same low-dimensional embeddings (PCA and UMAP) of the reference. This allows comparing different conditions in the same system of coordinates:

obj.projected <- Run.ProjecTILs(query = sample.list, ref = ref, filter.cells = FALSE)

[1] "Using assay RNA for isotype"
[1] "Aligning isotype to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space

Projecting corrected query onto Reference UMAP space
[1] "Using assay RNA for PDL1_treated"
[1] "Aligning PDL1_treated to reference map for batch-correction..."

Projecting corrected query onto Reference PCA space

Projecting corrected query onto Reference UMAP space

a <- plot.projection(ref = ref, query = obj.projected$isotype, linesize = 0.5) +
    theme(aspect.ratio = 1) + ggtitle("Isotype")
b <- plot.projection(ref = ref, query = obj.projected$PDL1_treated, linesize = 0.5) +
    theme(aspect.ratio = 1) + ggtitle("PDL1-treated")
a | b

As in the classifier mode, we can examine subtype composition

a <- plot.statepred.composition(ref = ref, query = obj.projected$isotype, metric = "Percent") +
    ggtitle("Isotype") + ylim(0, 40)


b <- plot.statepred.composition(ref = ref, query = obj.projected$PDL1_treated, metric = "Percent") +
    ggtitle("PDL1-treated") + ylim(0, 40)
a | b

Expression profiles for a panel of marker genes shows good agreement of the projected data with the reference

genes4radar <- c("Cxcr6", "Gzmb", "Ccl5", "Nkg7", "Ly6c2", "Cxcr5", "Tox", "Izumo1r",
    "Tnfsf8", "Ccr7", "Il7r", "Tcf7", "Eomes", "Gzmk", "Crtam", "Ifit1", "Foxp3",
    "Ikzf2", "Il2ra")

plot.states.radar(ref = ref, query = obj.projected, genes4radar = genes4radar, min.cells = 20)

Are there any differences between the control Th1 cells and the anti-PD-L1 treated Th1 cells?

library(EnhancedVolcano)
deg <- find.discriminant.genes(ref, query = obj.projected$PDL1_treated, query.control = obj.projected$isotype,
    state = "Th1_Effector")

EnhancedVolcano(deg, lab = rownames(deg), x = "avg_log2FC", y = "p_val", pCutoff = 0.01,
    FCcutoff = 0.2, labSize = 5, legendPosition = "none", drawConnectors = F, title = "Anti-PD-L1 vs. Isotype (Th1_Effector)")

We observe that Th1 effector cells after anti-PD-L1 treatment upregulated a Th1-associated gene module that includes Klrd1, Plac8, Ctla2a and Ly6c2. This observation, together with the amplification of Th1 effectors upon treatment, confirm the findings of the original study by Snell et al. (2021).

Effect of anti-PD-L1 therapy on CD4+ T cells in chronic infection