In [1]:
import warnings
warnings.filterwarnings("ignore")
# Import necessary libraries
import squidpy as sq
import scanpy as sc
import anndata as ad
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats # For statistical tests
import omnipath as op # Squidpy uses omnipath for LR pairs
print(f"Squidpy version: {sq.__version__}")
print(f"Scanpy version: {sc.__version__}")
print(f"Omnipath version: {op.__version__}")
adata = sc.read_h5ad("./102_Scvi_visium_results_Explanted1/processed_visium_adata_scvi.h5ad")
adata
import squidpy as sq
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# --- Configuration ---
# Set plotting parameters for better visualization
sc.set_figure_params(figsize=(8, 8), facecolor="white")
# --- Load your AnnData Object ---
# Replace this with your actual loading code if adata isn't already in memory
# For example:
# adata = sc.read_h5ad("path/to/your/processed_adata.h5ad")
# --- Verify Input Data ---
print("AnnData object info:")
print(adata)
# Check if the required inputs exist
if 'counts' not in adata.layers:
raise KeyError("The 'counts' layer is missing from your AnnData object.")
if 'leiden_scvi' not in adata.obs:
raise KeyError("The 'leiden_scvi' column is missing from adata.obs.")
if 'spatial' not in adata.obsm:
raise KeyError("The 'spatial' key is missing from adata.obsm (needed for some downstream visualizations).")
# Ensure cluster labels are categorical (often helpful for plotting)
if not pd.api.types.is_categorical_dtype(adata.obs['leiden_scvi']):
print("Converting 'leiden_scvi' to categorical.")
adata.obs['leiden_scvi'] = adata.obs['leiden_scvi'].astype('category')
# Define the cluster key we will use
cluster_key = "leiden_scvi"
print(f"\nUsing cluster key: '{cluster_key}'")
print(f"Available clusters: {adata.obs[cluster_key].cat.categories.tolist()}")
print(f"Using layer: 'counts'")
# --- Step 1: Run Ligand-Receptor Interaction Analysis ---
print("\nRunning squidpy.gr.ligrec()...")
# n_perms: Number of permutations for significance testing.
# Increase to 1000 or more for publication-quality p-values.
# 100 is faster for exploration.
# layer: Specify the layer containing count data.
# use_raw=False: Important when using 'layer' argument.
# copy=False: Modifies the adata object in place, adding results to adata.uns['ligrec'].
# seed: For reproducible permutations.
sq.gr.ligrec(
adata,
n_perms=1000, # Use 100 for speed, 1000+ for robust results
cluster_key=cluster_key,
layer="counts", # Specify the counts layer
use_raw=False, # Set to False when using the 'layer' argument
copy=False, # Modify adata in place (results in adata.uns['ligrec'])
seed=42 # For reproducibility
)
print("Ligand-receptor analysis complete. Results stored in adata.uns['leiden_scvi_ligrec'].")
# --- Step 2: Access and Explore the Results ---
# The results are stored in adata.uns['leiden_scvi_ligrec'] as a dictionary containing
# DataFrames for 'means', 'pvalues', and potentially 'deconvoluted_expr'.
# Check the structure of the results
print("\nStructure of results in adata.uns['leiden_scvi_ligrec']: ", list(adata.uns['leiden_scvi_ligrec'].keys()))
# Access the p-values and mean interaction strengths
pvals_df = adata.uns['leiden_scvi_ligrec']['pvalues']
means_df = adata.uns['leiden_scvi_ligrec']['means'] # This contains the mean expression product (Ligand_expr * Receptor_expr)
print(f"\nShape of pvalues DataFrame: {pvals_df.shape}")
print(f"Shape of means DataFrame: {means_df.shape}")
# Display the first few rows of the means DataFrame
print("\nHead of means DataFrame (interaction strength):")
print(means_df.head())
# Display the first few rows of the p-values DataFrame
print("\nHead of pvalues DataFrame:")
print(pvals_df.head())
# --- Step 3: Filter Significant Interactions ---
# Often, we are interested in interactions with low p-values.
alpha = 0.05 # Significance threshold
print(f"\nFiltering interactions with p-value < {alpha}")
# Create a boolean mask for significant p-values
significant_mask = pvals_df < alpha
# Apply the mask to get significant mean interaction strengths
# Replace non-significant values (where mask is False) with NaN
significant_means = means_df.where(significant_mask)
# Count significant interactions per cluster pair
n_significant_interactions = significant_mask.sum()
print("\nNumber of significant interactions per cluster pair:")
print(n_significant_interactions)
# Get the total number of significant interactions overall
total_significant = n_significant_interactions.sum()
print(f"\nTotal significant interactions found (p < {alpha}): {total_significant}")
Squidpy version: 1.6.5
Scanpy version: 1.11.1
Omnipath version: 1.0.9
AnnData object info:
AnnData object with n_obs × n_vars = 850 × 18074
obs: 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'pass_qc', '_scvi_batch', '_scvi_labels', 'leiden_scvi', '_scvi_raw_norm_scaling'
var: 'gene_ids', 'feature_types', 'genome', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm'
uns: '_scvi_manager_uuid', '_scvi_uuid', 'leiden_scvi', 'leiden_scvi_colors', 'neighbors', 'spatial', 'umap'
obsm: 'X_scVI', 'X_umap', 'spatial'
layers: 'counts'
obsp: 'connectivities', 'distances'
Using cluster key: 'leiden_scvi'
Available clusters: ['0', '1', '2', '3', '4', '5', '6']
Using layer: 'counts'
Running squidpy.gr.ligrec()...
Ligand-receptor analysis complete. Results stored in adata.uns['leiden_scvi_ligrec'].
Structure of results in adata.uns['leiden_scvi_ligrec']: ['means', 'pvalues', 'metadata']
Shape of pvalues DataFrame: (8715, 49)
Shape of means DataFrame: (8715, 49)
Head of means DataFrame (interaction strength):
cluster_1 0 \
cluster_2 0 1 2 3 4 5
source target
FYN TRPC4 0 0.228062 0.105378 0.059821 0.07817 0.042148
TRPC6 0.045082 0.089826 0.057433 0.050647 0.087344 0.052675
TRPV4 0.040984 0.04865 0 0.050647 0.04606 0.052675
KL TRPV5 0 0.00704 0 0 0 0
S100A10 TRPV5 0 0.261138 0 0 0 0
cluster_1 1 ... 5 \
cluster_2 6 0 1 2 ... 4
source target ...
FYN TRPC4 0.045606 0 0.326471 0.203787 ... 0.104442
TRPC6 0.048513 0.143491 0.188235 0.155842 ... 0.113617
TRPV4 0 0.139392 0.147059 0 ... 0.072332
KL TRPV5 0 0 0.014706 0 ... 0
S100A10 TRPV5 0 0 0.873529 0 ... 0
cluster_1 6 \
cluster_2 5 6 0 1 2 3
source target
FYN TRPC4 0.068421 0.071879 0 0.278386 0.155702 0.110145
TRPC6 0.078947 0.074786 0.095406 0.14015 0.107757 0.100971
TRPV4 0.078947 0 0.091308 0.098974 0 0.100971
KL TRPV5 0 0 0 0.017476 0 0
S100A10 TRPV5 0 0 0 0.39829 0 0
cluster_1
cluster_2 4 5 6
source target
FYN TRPC4 0.128494 0.092472 0.09593
TRPC6 0.137668 0.102999 0.098837
TRPV4 0.096384 0.102999 0
KL TRPV5 0 0 0
S100A10 TRPV5 0 0 0
[5 rows x 49 columns]
Head of pvalues DataFrame:
cluster_1 0 1 ... 5 \
cluster_2 0 1 2 3 4 5 6 0 1 2 ... 4
source target ...
FYN TRPC4 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN
TRPC6 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN
TRPV4 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN
KL TRPV5 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN
S100A10 TRPV5 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN
cluster_1 6
cluster_2 5 6 0 1 2 3 4 5 6
source target
FYN TRPC4 0.999 NaN NaN NaN NaN NaN NaN NaN NaN
TRPC6 0.982 NaN NaN NaN NaN NaN NaN NaN NaN
TRPV4 0.923 NaN NaN NaN NaN NaN NaN NaN NaN
KL TRPV5 NaN NaN NaN NaN NaN NaN NaN NaN NaN
S100A10 TRPV5 NaN NaN NaN NaN NaN NaN NaN NaN NaN
[5 rows x 49 columns]
Filtering interactions with p-value < 0.05
Number of significant interactions per cluster pair:
cluster_1 cluster_2
0 0 False
1 False
2 False
3 False
4 False
5 False
6 False
1 0 False
1 False
2 False
3 False
4 False
5 False
6 False
2 0 False
1 False
2 True
3 False
4 False
5 True
6 False
3 0 False
1 False
2 False
3 False
4 False
5 False
6 False
4 0 False
1 False
2 False
3 False
4 False
5 False
6 False
5 0 False
1 False
2 True
3 False
4 False
5 True
6 False
6 0 False
1 False
2 False
3 False
4 False
5 False
6 False
dtype: Sparse[bool, False]
Total significant interactions found (p < 0.05): 4