---
title: "COMMOTR: Cell-Cell Communication Analysis in Spatial Transcriptomics"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{COMMOTR Tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  eval = FALSE
)
```

## Introduction

COMMOTR (Cell-cell cOMmunication inference via cOllective opTimal tRansport) is an R package for inferring cell-cell communication (CCC) in spatial transcriptomics data. It implements the COMMOT algorithm from [Cang et al., Nature Methods 2023](https://www.nature.com/articles/s41592-022-01728-4).

### Key Features

- **Optimal Transport-based CCC Inference**: Models ligand-receptor interactions with spatial distance constraints
- **Communication Direction Analysis**: Infers spatial vector fields for signaling directions
- **Cluster-level Analysis**: Aggregates to cluster-cluster communication with statistical testing
- **Comprehensive Visualization**: Publication-ready plots including vector fields, networks, and chord diagrams
- **Seurat Integration**: Works seamlessly with Seurat V4/V5 workflow

## Installation

```{r install, eval=FALSE}
# Install from GitHub
devtools::install_github("Zaoqu-Liu/COMMOTR")
```

## Quick Start

```{r quickstart}
library(COMMOTR)
library(Seurat)
library(ggplot2)

# Load your spatial transcriptomics data
# seurat_obj <- readRDS("your_spatial_data.rds")
```

## Step 1: Load Ligand-Receptor Database

COMMOTR includes built-in databases from CellChat and CellPhoneDB:

```{r load_database}
# Load CellChat database for mouse
df_ligrec <- ligand_receptor_database(
    database = "CellChat",
    species = "mouse",
    signaling_type = "Secreted Signaling"
)

head(df_ligrec)
```

### Filter LR Pairs by Expression

Filter to keep only LR pairs with sufficient expression in your data:

```{r filter_lr}
df_ligrec_filtered <- filter_lr_database(
    df_ligrec = df_ligrec,
    seurat_obj = seurat_obj,
    min_cell_pct = 0.05,
    heteromeric = TRUE
)

cat("LR pairs after filtering:", nrow(df_ligrec_filtered), "\n")
```

## Step 2: Infer Spatial Communication

The main function `spatial_communication()` computes cell-cell communication matrices:

```{r spatial_comm}
seurat_obj <- spatial_communication(
    seurat_obj,
    df_ligrec = df_ligrec_filtered,
    database_name = "CellChat",
    dis_thr = 500,              # Distance threshold in spatial units
    cot_eps_p = 0.1,            # Entropy regularization
    cot_rho = 10,               # Penalty for unmatched mass
    pathway_sum = TRUE,         # Compute pathway-level totals
    n_workers = 4,              # Parallel workers
    verbose = TRUE
)
```

### Access Results

```{r access_results}
# Get all results
results <- get_communication_results(seurat_obj, "CellChat")

# Get specific communication matrix
total_comm <- get_communication_matrix(seurat_obj, "CellChat", pathway_name = "total")

# Get sender/receiver summaries
sender_df <- get_sender_receiver_df(seurat_obj, "CellChat", "sender")
```

## Step 3: Communication Direction Analysis

Compute spatial vector fields showing communication directions:
```{r comm_direction}
# Compute direction for a specific pathway
seurat_obj <- communication_direction(
    seurat_obj,
    database_name = "CellChat",
    pathway_name = "TGFb",
    k = 5  # Neighbors for smoothing
)

# Visualize with vector field
plot_cell_communication(
    seurat_obj,
    database_name = "CellChat",
    pathway_name = "TGFb",
    plot_method = "grid",
    color_by = "signal",
    arrow_scale = 1.5
)
```

### Spatial Autocorrelation of Direction

Test if communication directions are spatially coherent:

```{r spatial_autocorr}
moranI_result <- communication_spatial_autocorrelation(
    seurat_obj,
    database_name = "CellChat",
    key = "TGFb",
    field = "sender",
    n_permutations = 999
)

cat("Moran's I:", moranI_result$I, "\n")
cat("P-value:", moranI_result$p_value, "\n")
```

## Step 4: Cluster-Level Analysis

Aggregate communication to cluster-cluster level:

```{r cluster_comm}
seurat_obj <- cluster_communication(
    seurat_obj,
    database_name = "CellChat",
    clustering = "seurat_clusters",
    pathway_name = "TGFb",
    n_permutations = 100
)
```

### Visualize Cluster Communication

```{r plot_cluster}
# Network plot
plot_cluster_communication(
    seurat_obj,
    database_name = "CellChat",
    clustering = "seurat_clusters",
    key = "TGFb",
    layout = "circle",
    p_threshold = 0.05
)

# Heatmap
plot_communication_heatmap(
    seurat_obj,
    database_name = "CellChat",
    clustering = "seurat_clusters",
    key = "TGFb"
)

# Chord diagram (requires circlize)
plot_chord_diagram(
    seurat_obj,
    database_name = "CellChat",
    clustering = "seurat_clusters",
    key = "TGFb"
)
```

### Dotplot for Multiple Pathways

```{r dotplot}
# Compare multiple pathways
plot_communication_dotplot(
    seurat_obj,
    database_name = "CellChat",
    clustering = "seurat_clusters",
    keys = c("TGFb", "WNT", "FGF"),
    top_n_pairs = 15
)
```

## Step 5: Downstream Analysis

### Communication-Dependent Gene Detection

Identify genes whose expression varies with communication intensity:

```{r deg_detection}
deg_results <- communication_deg_detection(
    seurat_obj,
    database_name = "CellChat",
    key = "TGFb",
    signal_type = "sender",
    n_workers = 4
)

# Top significant genes
head(deg_results[deg_results$padj < 0.05, ], 20)
```

### Communication Impact Analysis

Analyze which signaling pathways impact specific genes:

```{r comm_impact}
impact_results <- communication_impact(
    seurat_obj,
    database_name = "CellChat",
    genes = c("Acta2", "Col1a1", "Fn1"),  # Example fibrosis markers
    method = "partial_correlation"
)

# View top impacts
head(impact_results[order(-abs(impact_results$impact)), ])
```

### Group Cells by Communication Patterns

```{r group_cells}
seurat_obj <- group_cell_communication(
    seurat_obj,
    database_name = "CellChat",
    key = "total",
    n_groups = 5
)

# Visualize groups
DimPlot(seurat_obj, group.by = "commgroup_CellChat_total")
```

## Parameters Guide

### Key Parameters for `spatial_communication()`

| Parameter | Description | Recommended |
|-----------|-------------|-------------|
| `dis_thr` | Maximum communication distance | Based on tissue type (100-1000) |
| `cot_eps_p` | Entropy regularization | 0.1 (default) |
| `cot_rho` | Unmatched mass penalty | 10 (default) |
| `cot_strategy` | COT strategy | "cot" for most cases |

### COT Strategies

- **cot**: Standard collective OT, aggregates all marginals
- **row**: Row-wise, independent optimization per sender
- **col**: Column-wise, independent optimization per receiver
- **blk**: Block-wise, per cell pair

## Citation

If you use COMMOTR in your research, please cite:

> Cang, Z., Zhao, Y., Almet, A.A. et al. Screening cell–cell communication in spatial transcriptomics via collective optimal transport. *Nat Methods* 20, 218–228 (2023). https://doi.org/10.1038/s41592-022-01728-4

## Session Info

```{r session}
sessionInfo()
```