Quantifying tumor-infiltrating lymphocytes (TILs)

Introduction¶

This tutorial demonstrates a WSI-level analysis workflow focused on tumor-infiltrating lymphocytes (TILs). We will use cervix biopsy WSI segmentation data as an example. Specifically, in this tutorial, we will show:

Spatial Querying of TILs: Spatial querying allows you to precisely locate TILs relative to other tissue compartments. We will show how to:

1. Query for intratumoral TILs: By performing a spatial containment query, you can identify all lymphocytes located entirely within the boundaries of a segmented tumor region.

2. Identify TILs at the tumor-stroma interface: Using tissue intersection partitioning, you can select all lymphocytes that are located within a specified distance of the tumor boundary.

Quantifying Lymphocyte "Hotness" and "Coldness": Beyond a simple count, the spatial distribution of TILs provides valuable prognostic information. A "hot" tumor is characterized by a high density of TILs, while a "cold" tumor has a low or absent immune infiltrate. We will demonstrate how to quantify this using grid based aggregation., where we partition the stromal regions into small, uniform grid cells and aggregating the immune cell densities within those.

Quantifying Neoplastic-Lymphocyte Links: We will create a spatial graph where nodes represent individual nuclei (both tumor cells and lymphocytes). Links (edges) will be established based on a proximity rule. We'll then show how to to identify and count tumor-lymphocyte links. This gives you a direct measure of close proximity between tumor cells and lymphocytes.

In [1]:

import matplotlib.pyplot as plt
from histolytics.data import cervix_nuclei, cervix_tissue

# Let's load example data
tis = cervix_tissue()
nuc = cervix_nuclei()

fig, ax = plt.subplots(figsize=(10, 10))
tis.plot(ax=ax, column="class_name", aspect=1, alpha=0.5, legend=True)
nuc.plot(ax=ax, column="class_name", aspect=1, legend=False)
ax.set_axis_off()

import matplotlib.pyplot as plt from histolytics.data import cervix_nuclei, cervix_tissue # Let's load example data tis = cervix_tissue() nuc = cervix_nuclei() fig, ax = plt.subplots(figsize=(10, 10)) tis.plot(ax=ax, column="class_name", aspect=1, alpha=0.5, legend=True) nuc.plot(ax=ax, column="class_name", aspect=1, legend=False) ax.set_axis_off()

Intratumoral TILs¶

Querying intratumoral TILs is straightforward. By performing a spatial containment query, you can identify all lymphocytes located entirely within the boundaries of a segmented tumor region. In this example, we are actually talking about a cervical lesion and pre-cancerous cells but the principles remain the same.

In [2]:

from histolytics.spatial_ops import get_objs

# get the CIN (lesion) tissue
cin_tissue = tis[tis["class_name"] == "cin"]

# select all the nuclei contained within CIN tissue
nuc_within_cin = get_objs(cin_tissue, nuc, predicate="contains")

ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False)
nuc_within_cin.plot(ax=ax, column="class_name", aspect=1, legend=True)
ax.set_axis_off()

from histolytics.spatial_ops import get_objs # get the CIN (lesion) tissue cin_tissue = tis[tis["class_name"] == "cin"] # select all the nuclei contained within CIN tissue nuc_within_cin = get_objs(cin_tissue, nuc, predicate="contains") ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False) nuc_within_cin.plot(ax=ax, column="class_name", aspect=1, legend=True) ax.set_axis_off()

Next we can count the intratumoral TILs which can be a clinically relevant feature to assess.

In [3]:

nuc_within_cin.value_counts("class_name")

nuc_within_cin.value_counts("class_name")

Out[3]:

class_name
neoplastic           3501
inflammatory          540
connective              3
glandular_epithel       1
squamous_epithel        1
Name: count, dtype: int64

Visualizing the TILs:

In [4]:

TILS = nuc_within_cin.loc[nuc_within_cin["class_name"] == "inflammatory"]

ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False)
TILS.plot(
    ax=ax,
    column="class_name",
    aspect=1,
    legend=True,
)
ax.set_axis_off()

TILS = nuc_within_cin.loc[nuc_within_cin["class_name"] == "inflammatory"] ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False) TILS.plot( ax=ax, column="class_name", aspect=1, legend=True, ) ax.set_axis_off()

TILs at the Tumor-Stroma Interface¶

Now we'll identify TILs that are located at the tumor-stroma interface (actually lesion-stroma-interface). This can be done by first partitioning the stromal tissue into the tumor-stroma-interface and then performing a spatial query to find all lymphocytes that are within the interface.

In [5]:

from histolytics.spatial_ops import get_interfaces

stroma = tis[tis["class_name"] == "stroma"]
cin_tissue = tis[tis["class_name"] == "cin"]

# Partition interface
interface = get_interfaces(cin_tissue, stroma, buffer_dist=300)
interface = interface.assign(class_name="lesion-stroma-interface")

# get the nuclei within the interface
nuc_within_interface = get_objs(interface, nuc, predicate="contains")

# plot
ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False)
interface.plot(ax=ax, column="class_name", aspect=1, alpha=0.3, legend=True)
nuc_within_interface.plot(
    ax=ax,
    column="class_name",
    aspect=1,
    legend=False,
)
ax.set_axis_off()

# Get the absolute inflammatory counts or the fractions by setting `normalize=True`
nuc_within_interface.value_counts("class_name", normalize=False)

from histolytics.spatial_ops import get_interfaces stroma = tis[tis["class_name"] == "stroma"] cin_tissue = tis[tis["class_name"] == "cin"] # Partition interface interface = get_interfaces(cin_tissue, stroma, buffer_dist=300) interface = interface.assign(class_name="lesion-stroma-interface") # get the nuclei within the interface nuc_within_interface = get_objs(interface, nuc, predicate="contains") # plot ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False) interface.plot(ax=ax, column="class_name", aspect=1, alpha=0.3, legend=True) nuc_within_interface.plot( ax=ax, column="class_name", aspect=1, legend=False, ) ax.set_axis_off() # Get the absolute inflammatory counts or the fractions by setting `normalize=True` nuc_within_interface.value_counts("class_name", normalize=False)

Out[5]:

class_name
inflammatory         2820
connective           1141
glandular_epithel      38
neoplastic              7
squamous_epithel        1
Name: count, dtype: int64

Now select again the immune nuclei within the interface and plot them.

In [6]:

TILS = nuc_within_interface.loc[nuc_within_interface["class_name"] == "inflammatory"]

ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False)
interface.plot(ax=ax, column="class_name", aspect=1, alpha=0.3, legend=True)
TILS.plot(
    ax=ax,
    column="class_name",
    aspect=1,
    legend=True,
)
ax.set_axis_off()

TILS = nuc_within_interface.loc[nuc_within_interface["class_name"] == "inflammatory"] ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False) interface.plot(ax=ax, column="class_name", aspect=1, alpha=0.3, legend=True) TILS.plot( ax=ax, column="class_name", aspect=1, legend=True, ) ax.set_axis_off()

As we see, there is a heavy concentration of TILs at the interface region. This can be a sign of an activated immune response to the CIN lesion.

Immune Hotness/Coldness¶

Next, we'll quantify the immune "hotness" and "coldness" of the tumor microenvironment by analyzing the spatial density distribution of the TILs. We will use the h3 spatial indexing to create a hexagonal grid over the stromal region and compute the density of TILs within each hexagon.

In [7]:

from histolytics.spatial_ops.h3 import h3_grid

h3 = h3_grid(stroma, resolution=10)
ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False)

# Let's plot the grid to see the hexagonal partitioning
h3.plot(
    ax=ax, aspect=1, legend=True, facecolor="none", edgecolor="blue", lw=1, alpha=0.5
)
ax.set_axis_off()

h3.head(5)

from histolytics.spatial_ops.h3 import h3_grid h3 = h3_grid(stroma, resolution=10) ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False) # Let's plot the grid to see the hexagonal partitioning h3.plot( ax=ax, aspect=1, legend=True, facecolor="none", edgecolor="blue", lw=1, alpha=0.5 ) ax.set_axis_off() h3.head(5)

Out[7]:

	geometry
8a82a9395007fff	POLYGON ((6754.29642 740.90756, 6795.0466 681....
8a82a939a0cffff	POLYGON ((3777.49042 7514.3929, 3818.20685 745...
8a82a939a227fff	POLYGON ((3270.67389 7313.52131, 3344.20357 73...
8a82a939e217fff	POLYGON ((4353.07966 5392.93377, 4385.89665 54...
8a82a939a787fff	POLYGON ((4415.58687 7992.8303, 4448.40275 806...

Next we'll compute the immune densities. We will use the fraction of the inflammatory nuclei as the density metric but you could also use other metrics such as the absolute count or absolute count divided by the area of the grid cells.

In [8]:

from histolytics.spatial_agg import grid_aggregate
from histolytics.utils.plot import legendgram


# This function will compute the percentage of immune nuclei within each grid cell
# In general, any function that takes a GeoDataFrame and returns a scalar
# can be used here. Typically, this will be a function that calculates
# a count, sum, mean, or other statistic of interest out of the nuclei.
def immune_density(nuclei):
    """Calculate the immune cell count in a grid cell."""
    if "inflammatory" in nuclei.value_counts("class_name"):
        cnt = nuclei.value_counts("class_name", normalize=False)["inflammatory"]
    else:
        cnt = 0
    return float(cnt)


h3 = grid_aggregate(
    objs=nuc,
    grid=h3,
    metric_func=immune_density,
    new_col_names=["immune_density"],
    predicate="contains",
    num_processes=2,
)

ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False)
nuc.plot(ax=ax, column="class_name", aspect=1, legend=True)
h3.plot(
    ax=ax,
    column="immune_density",
    cmap="Reds",
    legend=True,
    aspect=1,
    facecolor="none",
)
ax.set_axis_off()

# Add a legendgram to visualize the distribution of immune cell density
ax = legendgram(
    gdf=h3,
    column="immune_density",
    n_bins=30,
    cmap="Reds",
    ax=ax,
)

from histolytics.spatial_agg import grid_aggregate from histolytics.utils.plot import legendgram # This function will compute the percentage of immune nuclei within each grid cell # In general, any function that takes a GeoDataFrame and returns a scalar # can be used here. Typically, this will be a function that calculates # a count, sum, mean, or other statistic of interest out of the nuclei. def immune_density(nuclei): """Calculate the immune cell count in a grid cell.""" if "inflammatory" in nuclei.value_counts("class_name"): cnt = nuclei.value_counts("class_name", normalize=False)["inflammatory"] else: cnt = 0 return float(cnt) h3 = grid_aggregate( objs=nuc, grid=h3, metric_func=immune_density, new_col_names=["immune_density"], predicate="contains", num_processes=2, ) ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False) nuc.plot(ax=ax, column="class_name", aspect=1, legend=True) h3.plot( ax=ax, column="immune_density", cmap="Reds", legend=True, aspect=1, facecolor="none", ) ax.set_axis_off() # Add a legendgram to visualize the distribution of immune cell density ax = legendgram( gdf=h3, column="immune_density", n_bins=30, cmap="Reds", ax=ax, )

Here we see that the immune dense regions are located at the interface between the CIN lesion and the stroma, indicative of immune activation. Next we'll compute the distances of the different density grid cells to the lesion to actually quantify that the immune dense regions are indeed closer to the lesion.

In [9]:

import pandas as pd
import mapclassify

lesion = tis[tis["class_name"] == "cin"]

distances = {}
for i, lesion in lesion.reset_index().iterrows():
    dist = h3.distance(lesion.geometry)
    distances[i] = dist

min_dists = pd.DataFrame(distances).min(axis=1)
min_dists.name = "min_dist"

# join the distances to the grid
h3 = h3.join(other=min_dists, how="left")

# Let's bin the immune density for visualization, we'll use quantile binning
bins = mapclassify.Quantiles(h3["immune_density"], k=4)
h3["immune_density_level"] = bins.yb

h3

import pandas as pd import mapclassify lesion = tis[tis["class_name"] == "cin"] distances = {} for i, lesion in lesion.reset_index().iterrows(): dist = h3.distance(lesion.geometry) distances[i] = dist min_dists = pd.DataFrame(distances).min(axis=1) min_dists.name = "min_dist" # join the distances to the grid h3 = h3.join(other=min_dists, how="left") # Let's bin the immune density for visualization, we'll use quantile binning bins = mapclassify.Quantiles(h3["immune_density"], k=4) h3["immune_density_level"] = bins.yb h3

Out[9]:

	geometry	immune_density	min_dist	immune_density_level
8a82a9395007fff	POLYGON ((6754.29642 740.90756, 6795.0466 681....	0.0	3210.127874	0
8a82a939a0cffff	POLYGON ((3777.49042 7514.3929, 3818.20685 745...	0.0	1539.781970	0
8a82a939a227fff	POLYGON ((3270.67389 7313.52131, 3344.20357 73...	0.0	1000.778838	0
8a82a939e217fff	POLYGON ((4353.07966 5392.93377, 4385.89665 54...	0.0	1663.882983	0
8a82a939a787fff	POLYGON ((4415.58687 7992.8303, 4448.40275 806...	0.0	1943.510700	0
...	...	...	...	...
8a82a92b42dffff	POLYGON ((2321.10708 10094.14135, 2280.40587 1...	2.0	0.000000	2
8a82a92b0927fff	POLYGON ((2435.32504 10045.4735, 2508.84185 10...	0.0	0.000000	0
8a82a92b429ffff	POLYGON ((2615.17352 10135.55733, 2541.65625 1...	3.0	0.000000	3
8a82a92b465ffff	POLYGON ((2746.43098 10413.05959, 2713.61664 1...	10.0	0.000000	3
8a82a92b09b7fff	POLYGON ((2705.72793 10472.07979, 2746.43098 1...	1.0	0.000000	1

2313 rows × 4 columns

Let's now plot the distributions of the distances to the lesion for each immune density level.

In [10]:

# !pip install seaborn # Uncomment this line to install seaborn for plotting

# !pip install seaborn # Uncomment this line to install seaborn for plotting

In [ ]:

import seaborn as sns

tidy = h3.reset_index().set_index("immune_density_level")
tidy = tidy[["min_dist"]]
tidy = tidy.stack()
tidy = tidy.reset_index()
tidy = tidy.rename(
    columns={
        "immune_density_level": "Immune Density Level",
        "level_1": "Attribute",
        0: "Distance to Lesion",
    }
)

fig, ax = plt.subplots(1, 2, figsize=(10, 5))

ax[0] = sns.kdeplot(
    ax=ax[0],
    data=tidy,
    x="Distance to Lesion",
    hue="Immune Density Level",
    fill=True,
    alpha=0.5,
    palette="viridis",
)

ax[1] = sns.swarmplot(
    ax=ax[1],
    data=tidy,
    y="Distance to Lesion",
    x="Immune Density Level",
    hue="Immune Density Level",
    size=1.5,
    orient="v",
    legend=False,
    warn_thresh=0.5,
    palette="viridis",
)

ax[1] = sns.boxplot(
    ax=ax[1],
    data=tidy,
    x="Immune Density Level",
    y="Distance to Lesion",
    showfliers=False,
    color="black",
    linewidth=1.5,
    width=0.2,
    fill=False,
    whis=1.0,
    showcaps=False,
)

import seaborn as sns tidy = h3.reset_index().set_index("immune_density_level") tidy = tidy[["min_dist"]] tidy = tidy.stack() tidy = tidy.reset_index() tidy = tidy.rename( columns={ "immune_density_level": "Immune Density Level", "level_1": "Attribute", 0: "Distance to Lesion", } ) fig, ax = plt.subplots(1, 2, figsize=(10, 5)) ax[0] = sns.kdeplot( ax=ax[0], data=tidy, x="Distance to Lesion", hue="Immune Density Level", fill=True, alpha=0.5, palette="viridis", ) ax[1] = sns.swarmplot( ax=ax[1], data=tidy, y="Distance to Lesion", x="Immune Density Level", hue="Immune Density Level", size=1.5, orient="v", legend=False, warn_thresh=0.5, palette="viridis", ) ax[1] = sns.boxplot( ax=ax[1], data=tidy, x="Immune Density Level", y="Distance to Lesion", showfliers=False, color="black", linewidth=1.5, width=0.2, fill=False, whis=1.0, showcaps=False, )

As seen in the KDE and swarm-plots, the immune dense regions are clearly closer to the lesion compared to the less dense regions (which is also clearly visible in the spatial plots). This is a nice way to quantify the spatial distribution of immune cells in relation to the tumor (lesion) tissue in segmented histology images.

Neoplastic-Lymphocyte Links¶

Next we'll quantify direct proximity based links between neoplastic cells and lymphocytes. This will involve creating a spatial graph where nodes represent individual nuclei (both tumor cells and lymphocytes) and edges represent proximity relationships. This way we can have a more detailed view of possible interactions between neoplastic and immune cells.

In [ ]:

from histolytics.spatial_graph.graph import fit_graph
from histolytics.utils.gdf import set_uid

nuc = set_uid(nuc)  # ensure unique IDs for nuclei
w, w_gdf = fit_graph(
    nuc, "delaunay", id_col="uid", threshold=100
)  # Let's use Delaunay graph

# visualize the different links on the tissue
ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.3)
w_gdf.plot(ax=ax, linewidth=1, column="class_name", legend=True, aspect=1, lw=0.5)
ax.set_axis_off()
w_gdf.value_counts("class_name")

from histolytics.spatial_graph.graph import fit_graph from histolytics.utils.gdf import set_uid nuc = set_uid(nuc) # ensure unique IDs for nuclei w, w_gdf = fit_graph( nuc, "delaunay", id_col="uid", threshold=100 ) # Let's use Delaunay graph # visualize the different links on the tissue ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.3) w_gdf.plot(ax=ax, linewidth=1, column="class_name", legend=True, aspect=1, lw=0.5) ax.set_axis_off() w_gdf.value_counts("class_name")

Out[ ]:

class_name
connective-connective                  15216
connective-inflammatory                12399
inflammatory-inflammatory              11953
neoplastic-neoplastic                   8470
inflammatory-neoplastic                 2526
glandular_epithel-glandular_epithel     1142
connective-neoplastic                    780
glandular_epithel-inflammatory           439
connective-glandular_epithel             382
glandular_epithel-neoplastic              80
connective-squamous_epithel                5
inflammatory-squamous_epithel              2
glandular_epithel-squamous_epithel         1
neoplastic-squamous_epithel                1
squamous_epithel-squamous_epithel          1
Name: count, dtype: int64

In total, there are more than 2500 links between neoplastic and inflammatory cells. This link count can be used as a metric for immune infiltration as it is or you can derive additional features such as link density at different regions etc. Let's now select these links and visualize them.

In [14]:

imm_neo_links = w_gdf[w_gdf["class_name"] == "inflammatory-neoplastic"]
ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.3)
imm_neo_links.plot(ax=ax, linewidth=0.5, column="class_name", legend=True, aspect=1)
ax.set_axis_off()

imm_neo_links = w_gdf[w_gdf["class_name"] == "inflammatory-neoplastic"] ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.3) imm_neo_links.plot(ax=ax, linewidth=0.5, column="class_name", legend=True, aspect=1) ax.set_axis_off()

As we see, these links between neoplastic and inflammatory cells concentrate on specific regions on the lesion. We can check which parts of the lesion are the most enriched for these links by overlaying a grid on the lesion and computing the link density. It should be noted that one inflammatory cell can be linked to multiple neoplastic cells and vice versa, so the link density will give us a sense of the overall interaction landscape.

In [ ]:

def immune_neo_link_density(nuclei):
    """Calculate the immune cell count in a grid cell."""
    if "inflammatory-neoplastic" in nuclei.value_counts("class_name"):
        cnt = nuclei.value_counts("class_name", normalize=False)[
            "inflammatory-neoplastic"
        ]
    else:
        cnt = 0
    return float(cnt)


lesion = tis[tis["class_name"] == "cin"]
h3_lesion = h3_grid(lesion, resolution=10)

h3_lesion = grid_aggregate(
    objs=w_gdf,
    grid=h3_lesion,
    metric_func=immune_neo_link_density,
    new_col_names=["immune_link_density"],
    predicate="contains",
    num_processes=2,
)

ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False)
imm_neo_links.plot(ax=ax, column="class_name", aspect=1, legend=True, lw=0.3)
h3_lesion.plot(
    ax=ax,
    column="immune_link_density",
    cmap="Reds",
    legend=True,
    aspect=1,
    facecolor="none",
)
ax.set_axis_off()

# Add a legendgram to visualize the distribution of immune cell density
ax = legendgram(
    gdf=h3_lesion,
    column="immune_link_density",
    n_bins=15,
    cmap="Reds",
    ax=ax,
)

def immune_neo_link_density(nuclei): """Calculate the immune cell count in a grid cell.""" if "inflammatory-neoplastic" in nuclei.value_counts("class_name"): cnt = nuclei.value_counts("class_name", normalize=False)[ "inflammatory-neoplastic" ] else: cnt = 0 return float(cnt) lesion = tis[tis["class_name"] == "cin"] h3_lesion = h3_grid(lesion, resolution=10) h3_lesion = grid_aggregate( objs=w_gdf, grid=h3_lesion, metric_func=immune_neo_link_density, new_col_names=["immune_link_density"], predicate="contains", num_processes=2, ) ax = tis.plot(figsize=(10, 10), column="class_name", aspect=1, alpha=0.5, legend=False) imm_neo_links.plot(ax=ax, column="class_name", aspect=1, legend=True, lw=0.3) h3_lesion.plot( ax=ax, column="immune_link_density", cmap="Reds", legend=True, aspect=1, facecolor="none", ) ax.set_axis_off() # Add a legendgram to visualize the distribution of immune cell density ax = legendgram( gdf=h3_lesion, column="immune_link_density", n_bins=15, cmap="Reds", ax=ax, )

Conclusions¶

In this workflow tutorial, we showed how to analyze the immuno oncological relationships between inflammatory cells and tumor cells in a segmented histology WSI using Histolytics. We demonstrated how to quantify intratumoral TILs, tumor-stroma-interface TILs, immune cell densities in the stroma and immune cell link densities within the lesion. The flexible tools in Histolytics enable easy exploration, quantification and visualization of these complex spatial relationships.