Skip to content

Add link validation; rfc full workflow #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
JavenTyr wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add link validation; rfc full workflow #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
162 changes: 60 additions & 102 deletions linkage_analysis/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,109 +1,67 @@
# Prerequisites
### Clone projects and restore environment:
```r
git clone https://github.com/stuart-lab/metapeak-analysis.git
cd metapeak-analysis/linkage_analysis
renv::restore()
# REMO analysis

This repository contains code to reproduce results shown in *Regulatory element modules as universal features for single-cell chromatin analysis*:

Chrysania Lim, Javen Tan Yih Ruay, and Tim Stuart, 2025. Regulatory element modules as universal features for single-cell chromatin analysis. *bioRxiv*. https://doi.org/10.64898/2025.12.10.692786

For code used to generate the REMO modules, please see this repository: https://github.com/stuart-lab/REMO

For the R package containing REMO modules and associated metadata, please see this repository: https://github.com/stuart-lab/REMO.v1.GRCh38

## Expected inputs

Recommended project layout:

```text
resources/
gencode.v32.basic.annotation.gtf.gz
GTEx_49tissues_release1.tsv.gz
multiome/
<dataset>_multiome_linkpeaks.rds
scripts/
config/
objects_and_links.csv
link_filtering_conds.csv
data/
multiome_links/
<dataset>/
chr1.tsv
chr2.tsv
...
all_chromosomes.tsv
combined_metadata/
chr1.tsv
chr2.tsv
...
filtered_links/
<condition>/
pgl_metadata/pgls_metadata_all_chr.tsv
genes_per_peak/genes_per_peak_all_chr.tsv
gtex_validation/
```
For restoring renv refer to https://rstudio.github.io/renv/articles/renv.html#collaboration

### Prepare resources folder
See below for expected structure:
```bash
./linkage_analysis/resources
├── example
│   └── pgls_metadata_all_chr.tsv.gz (from supplementary data)
├── gencode.v32.basic.annotation.gtf.gz (from https://www.gencodegenes.org/human/release_32.html)
└── multiome (your Seurat multiome datasets[.rds])
```

# Script setup
##### Modify slurm parameters as needed, e.g.(`cpus-per-task`, `mem-per-cpu`)
- `/scripts/01_link_peaks.sh`
- `/scripts/02_run_aggregate_links.sh`
- `/scripts/03_run_probability_analysis.sh`

##### copy pwd to `PROJECT_ROOT` of .sh files
```bash
pwd
$ /metapeak-analysis/linkage_analysis/
```

# Reproduce results
In `/scripts/03_run_probability_analysis.sh`, change `SIG_PGL_FPATH` variable to `"${PROJECT_ROOT}/resources/example/pgl_metadata_all_chr.tsv"`

Run the following script
```bash
./scripts/03_run_probability_analysis.sh
```

# Run full analysis workflow on user-supplied data
### 1. Replace examples files in `/resources/multiome/` with your intended Seurat multiome objects (`.rds`)

### 2. Modify the following variables in `01_link_peaks.sh` according to each seurat multiome object
- `OUTPUT_DIR` (links for each chromosome are stored here, in `/data/multiome_links/dataset_name`)
- `SEURAT_RDS` (filepath for Seurat multiome object, in `/resources/multiome`)

E.g.
```bash
OUTPUT_DIR="${PROJECT_ROOT}/data/multiome_links/jejunum"
SEURAT_RDS="${PROJECT_ROOT}/resources/multiome/jejunum_multiome_linkpeaks.rds"
```

### 3. Run the follow script:
```bash
./01_link_peaks.sh
```
Repeat for each Seurat multiome object.
`objects_and_links.csv` is used by step 2 for both aggregation and GTEx/SGL validation. It must contain a `links_dir` column pointing to each dataset's `data/multiome_links/<dataset>` directory. A `dataset` column is optional; if omitted, the dataset label is inferred from the link directory name.

### 4. Create the following config files in `/scripts/configs`
##### `objects_and_links.csv` (modify example table of seurat obj filepaths and their mutiome_link dirs)
```csv
seurat_obj_fpath,links_dir
/resources/multiome/jejunum_multiome_linkpeaks.rds,/data/multiome_links/jejunum/
/resources/multiome/bile_duct_multiome_linkpeaks.rds,/data/multiome_links/bile_duct/
```

##### `link_filtering_conds.csv` (specify p-value and score thresholds to filter peak-gene links for downstream analysis)
```csv
condition,pvalue,score
sig_pvalue,0.05,0
hsig_pvalue_high_coeff,0.01,0.1
```
The GTEx/SGL fine-mapping file used by the validation part of step 2 should be placed at `resources/GTEx_49tissues_release1.tsv.gz` by default. It is expected to contain at least the columns `gene`, `chromosome`, `end`, `tissue`, and `pip`.

### 5. Run the following script:
```bash
./scripts/02_run_aggregate_links.sh
```
## Workflow

### 6. Modify `/scripts/03_run_probability_analysis.sh`
Specify the `SIG_PGL_FPATH` variable as the path to the combined pgl_metadata for the condition of choice (e.g. `sig_pvalue`), example below:
```bash
${PROJECT_ROOT}/data/filtered_links/sig_pvalue/pgl_metadata/pgls_metadata_all_chr.tsv
```

### 7. Run the following script:
```bash
./scripts/03_run_probability_analysis.sh
```
1. `01_link_peaks.R` / `01_link_peaks.sh`
- Run Signac V2-compatible peak-gene linking per dataset.
- The supplied GENCODE GTF is attached to the peak assay.
- Genes with non-unique symbols are excluded before candidate link discovery.
- Per-chromosome link TSVs, `all_chromosomes.tsv`, and a linked Seurat RDS are written.

# Results
2. `02_aggregate_links.R`, `02a_validate_links_against_gtex.R` / `02_aggregate_and_validate_links.sh`
- Run aggregation first, then GTEx/SGL validation using the aggregated outputs.
- `02_aggregate_links.R` combines per-chromosome link TSVs across datasets, adds REMO metadata, creates per-condition filtered PGL metadata, and builds proximal/distal/all gene lists per peak.
- `02a_validate_links_against_gtex.R` validates link metrics against GTEx/SGL fine-mapping positives.
- Uses `scripts/config/objects_and_links.csv` for per-dataset validation.
- Uses `data/combined_metadata/chr*.tsv`, written by the aggregation step, for all-dataset validation.
- Writes per-dataset outputs under `data/gtex_validation/per_dataset/<dataset>/`.
- Writes all-dataset outputs under `data/gtex_validation/all_datasets/`.
- For each dataset and for all aggregated links, writes ROC plots, AUC tables, GTEx-labelled link tables, and diagnostics.
- Metrics evaluated: `score`, `zscore`, `exp(-distance / 200 kb)`, and `score * exp(-distance / 200 kb)`. The redundant `1 / distance` metric is not used.

Output structure:
```bash
./linkage_analysis/data
├── combined_metadata
├── filtered_links
├── multiome_links
└── probability_analysis
```

A table of conditional probabilities can be found at:
```
data/probability_analysis/conditional_prob_unlinked_singleCRE.tsv
```

The per-chromosome odds ratio for the likelihood that CRE pairs linked to the same gene are within the same REMO module can be found at:
```
data/probability_analysis/logistic_regression_results/glm_results_all.tsv
```
3. `03_get_conditional_probs.R`, `03a_get_glm_coeffs.R`, `03b_build_glm_summary.R` / `03_run_probability_analysis.sh`
- Run conditional probability and GLM analyses on a selected filtered PGL metadata file.
Loading