CCI network map for global LR interactions between cell types

[Rafael-Silva-Oliveira]

Is your feature request related to a problem? Please describe.

Would be super cool to expand some of the plotting functionality of LIANA+ and make it even more of a strong contestant for the go-to tool for CCC/LR analysis for both scRNA and Spatial transcriptomics, so this request would be to create a plot like this (adapted from the ccinet_plot from stLearn, which unfortunately isn't quite compatible with HD spatial data nor well maintained - https://stlearn.readthedocs.io/en/stable/tutorials/stLearn-CCI.html)

Describe the solution you'd like From their plotting functionality, something like this:

Contribution

I've made an attempt at adapting the code using a 10X Visium HD dataset after processing with Bin2Cell to extract actual cell locations and using CellTypist for cell type inference

from matplotlib import pyplot as plt
import pandas as pd
import numpy as np
import networkx as nx
import matplotlib.patches as patches

def CCIMap(
    interaction_data,
    lr: str = None,
    pos: dict = None,
    return_pos: bool = False,
    cmap: str = "tab10",
    font_size: int = 12,
    node_size_exp: int = 1,
    node_size_scaler: int = 1,
    min_counts: int = 0,
    sig_interactions: bool = True,
    fig: plt.Figure = None,
    ax: plt.Axes = None,
    pad=0.25,
    title: str = None,
    figsize: tuple = (10, 10),
):
    """Circular celltype-celltype interaction network based on LR-CCI analysis using LIANA results adapted from stLearn (https://stlearn.readthedocs.io/en/stable/tutorials/stLearn-CCI.html)."""
    # Filter interactions based on the LR pair if specified
    if sig_interactions:
        interaction_data = interaction_data[
            (interaction_data["cellphone_pvals"] < 0.05)
            & (interaction_data["lr_logfc"] > 0.5)
        ]

    if lr is not None:
        ligand, receptor = lr.split("^")
        interaction_data = interaction_data[
            (interaction_data["ligand_complex"] == ligand)
            & (interaction_data["receptor_complex"] == receptor)
        ]

    # Creating the interaction graph #
    graph = nx.MultiDiGraph()
    cell_types = np.unique(
        np.concatenate([interaction_data["source"], interaction_data["target"]])
    )

    for cell in cell_types:
        graph.add_node(cell)

    for _, row in interaction_data.iterrows():
        ligand = row["source"]
        receptor = row["target"]
        count = row["lr_means"]

        if count > min_counts:
            graph.add_edge(ligand, receptor, weight=count)

    if pos is None:
        pos = nx.circular_layout(graph)

    total = sum([d["weight"] for (u, v, d) in graph.edges(data=True)])

    if total == 0:
        print("No interactions found with the specified criteria.")
        return

    node_sizes = np.array(
        [
            (
                (sum([d["weight"] for _, _, d in graph.edges(cell, data=True)]) / total)
                * 10000
                * node_size_scaler
            )
            ** node_size_exp
            for cell in graph.nodes()
        ]
    )
    node_sizes[node_sizes == 0] = 0.1  # pseudocount

    edge_weights = [d["weight"] / total for (_, _, d) in graph.edges(data=True)]

    #### Drawing the graph #####
    if fig is None or ax is None:
        fig, ax = plt.subplots(figsize=figsize, facecolor=[0.7, 0.7, 0.7, 0.4])

    color_map = plt.get_cmap(cmap)
    node_colors = [color_map(i % color_map.N) for i, _ in enumerate(graph.nodes())]

    # Adding in the self-loops #
    z = 55
    for u, v, d in graph.edges(data=True):
        if u == v:
            x, y = pos[u]
            angle = np.degrees(np.arctan(y / x))
            if x > 0:
                angle += 180
            arc = patches.Arc(
                xy=(x, y),
                width=0.3,
                height=0.025,
                lw=5,
                ec=plt.cm.get_cmap("Blues")(d["weight"] / total),
                angle=angle,
                theta1=z,
                theta2=360 - z,
            )
            ax.add_patch(arc)

    nx.draw_networkx(
        graph,
        pos,
        node_size=node_sizes,
        node_color=node_colors,
        arrowstyle="->",
        arrowsize=50,
        width=5,
        font_size=font_size,
        font_weight="bold",
        edge_color=edge_weights,
        edge_cmap=plt.cm.Blues,
        ax=ax,
    )

    if title is None:
        title = f"CCI Interaction Network for ({lr})"
    fig.suptitle(title, fontsize=30)
    plt.tight_layout()

    xlims = ax.get_xlim()
    ax.set_xlim(xlims[0] - pad, xlims[1] + pad)
    ylims = ax.get_ylim()
    ax.set_ylim(ylims[0] - pad, ylims[1] + pad)

    if return_pos:
        return pos

The input would be the adata.uns["liana_res"] as follows:

ligrec_pairs = li.mt.cellphonedb(
    adata,
    groupby="cell_type",
    resource_name="consensus",
    expr_prop=0.1,
    verbose=True,
    use_raw=False,
    key_added="cpdb_res",
)

liana_results = li.mt.rank_aggregate(
    adata,
    groupby="cell_type",
    use_raw=False,
    groupby_pairs=adata.uns["cpdb_res"],
)

lr_pair = "TIMP1^CD63"
CCIMap(interaction_data=adata.uns["liana_res"], lr=lr_pair)

The output is:

Which seems to be very "uni-directional". It does seem that (as expected) a particular receptor is associated with a particular "source" cell type and thus we see one main cell type with the arrows pointing to all other cell types, instead of something like we see in the stLearn plot. Perhaps there's another function in Liana that would allow this visualization to be more insightful?

Here's the dataframe in case you just want to load it and test it out inter_data.xlsx

Happy to discuss this further!

Thanks

Other plots form stLearn that could be re-adapted and used in LIANA

Another possible implementation of the heatmap/dotplot plots, using lr_mean instead of # interactions as seen in the stLearn plot:

df = adata.uns["liana_res"]

# Create the new 'interaction' and 'cci' columns
df["interaction"] = df["ligand_complex"] + "-" + df["receptor_complex"]
df["cci"] = df["source"] + "-" + df["target"]

# Pivot the data for heatmap-style plotting
df_pivot = df.pivot_table(
    index="interaction", columns="cci", values="lr_means", aggfunc="mean"
)
df_pivot.fillna(0, inplace=True)
df_pivot
import plotly.express as px

fig = px.scatter(df, y="interaction", x="cci", color="lr_means")
fig.update_traces(marker_size=10)
fig.show()
plt.savefig("./dotplot_test.png")

EDIT: Ignore the last feature of the dotplot, completely forgot it already has one dotplot function which already shows the source and target as well as the interactions :)

[dbdimitrov]

Hi @Rafael-Silva-Oliveira,

Thanks a lot for the suggestion and the code. I hope that soon we'll get someone to work a bit on the visualisations of liana.

In the meantime, there is also this thread with code to make some really nice plots by users: https://github.com/saezlab/liana-py/issues/85

[Rafael-Silva-Oliveira]

Hi @Rafael-Silva-Oliveira,

Thanks a lot for the suggestion and the code. I hope that soon we'll get someone to work a bit on the visualisations of liana.

In the meantime, there is also this thread with code to make some really nice plots by users: #85

Super, thank you for mentioning that, I'll have a look! :)

[Rafael-Silva-Oliveira]

Also, @dbdimitrov what version of plotnine are you using? I'm getting an error here:

li.multi.nmf(lrdata, n_components=None, inplace=True, random_state=0, max_iter=200, verbose=True)

when n_components = None since this if statement prints an elbow plot and I'm getting an error which I believe is version issue; I have the latest plotnine version, so I'd like to check which one was yours just to re-run the code and see how many factors I should use :)

[dbdimitrov]

Hi @Rafael-Silva-Oliveira,

I should definitely change this - i.e. the plot should not be part of this function.

Re version: plotnine version = "0.12.4"

[Rafael-Silva-Oliveira]

Hi @Rafael-Silva-Oliveira,

I should definitely change this - i.e. the plot should not be part of this function.

Re version: plotnine version = "0.12.4"

I think the elbow plot makes a lot of sense there, just perhaps another if statement and a parameter associated with the function like "plot_elbow=True" per default and then the user can just change to False if they don't want to plot it

Thanks, works with that version :)