From c04d740d59c62c4e58cd6614cbc84423d2d29c1c Mon Sep 17 00:00:00 2001 From: Andreas Gschwind Date: Mon, 27 Jul 2026 09:06:14 -0700 Subject: [PATCH] Update direct vs. indirect effects analysis code Add additional analyses for supplementary figures added during revisions round 1. Add code to analyze indirect effect rates in the Gasperini et al., 2019 dataset to compare K562 DC-TAP-seq target elements and genes to genome-wide elements and genes. Add code to compute bootstrapped confidence intervals for the reported probability of direct effects. --- interpretation_analysis/README.md | 2 +- .../indirect_effects/README.md | 14 + .../bootstrapped_indirect_effects/README.md | 95 ++++ .../analyze_bootstrapped_directness.Rmd | 267 +++++++++ .../bootstrap_directness_CIs.R | 511 ++++++++++++++++++ .../dc_tapseq_indirect_effect_sizes.Rmd | 253 +++++++++ ...tapseq_indirect_effects_main_analyses.Rmd} | 91 +++- ...ndirect_effects_supplementary_analyses.Rmd | 286 ++++++++++ .../00_create_gene_universe.R | 54 ++ .../00_find_dctap_perts_in_gasperini.R | 51 ++ .../1_create_ondisc_sceptre_object.R | 65 +++ .../2_run_sceptre_nextflow.sh | 46 ++ .../3_calculate_indirect_effects.R | 103 ++++ .../4_analyze_indirect_effects.Rmd | 200 +++++++ .../nextflow.config | 5 + .../sceptre_formula.rds | Bin 0 -> 176 bytes 16 files changed, 2019 insertions(+), 24 deletions(-) create mode 100644 interpretation_analysis/indirect_effects/README.md create mode 100644 interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/README.md create mode 100644 interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/analyze_bootstrapped_directness.Rmd create mode 100644 interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/bootstrap_directness_CIs.R create mode 100644 interpretation_analysis/indirect_effects/dc_tapseq_indirect_effect_sizes.Rmd rename interpretation_analysis/{dc_tapseq_indirect_effects.Rmd => indirect_effects/dc_tapseq_indirect_effects_main_analyses.Rmd} (86%) create mode 100644 interpretation_analysis/indirect_effects/dc_tapseq_indirect_effects_supplementary_analyses.Rmd create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_create_gene_universe.R create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_find_dctap_perts_in_gasperini.R create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/1_create_ondisc_sceptre_object.R create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/2_run_sceptre_nextflow.sh create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/3_calculate_indirect_effects.R create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/4_analyze_indirect_effects.Rmd create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/nextflow.config create mode 100644 interpretation_analysis/indirect_effects/genome_wide_indirect_effects/sceptre_formula.rds diff --git a/interpretation_analysis/README.md b/interpretation_analysis/README.md index 11590620..5f2f6f3b 100644 --- a/interpretation_analysis/README.md +++ b/interpretation_analysis/README.md @@ -8,4 +8,4 @@ This directory contains code, intermediate files, and results used for the inter - `other_crispr_screens.categorized.tsv.gz`: same as above but for CRISPR perturbation data from 8 other studies - `all_gw_pairs.categorized.summary.tsv` - `thresholds.tsv`: quantiles for various chromatin features for distal elements in 5 cell types used for categorization - - `dc_tapseq_indirect_effects.Rmd`: explore and analyze probability of direct effects and create results shown in Figure 4 and Figure S3 + - `indirect_effects`: analyze probability of direct effects and create results shown in Figure 4 and Supplementary figures 3, 15 and 16 diff --git a/interpretation_analysis/indirect_effects/README.md b/interpretation_analysis/indirect_effects/README.md new file mode 100644 index 00000000..5b70fc50 --- /dev/null +++ b/interpretation_analysis/indirect_effects/README.md @@ -0,0 +1,14 @@ +### Direct versus indirect effects analyses + +This directory contains code used for direct vs. indirect effects analyses presented in Figure 4 and +Supplementary Figures 3, 15, and 16. Direct, indirect effect rates and probabilities of direct +effects are computed by the separate workflow https://github.com/EngreitzLab/CRISPR_indirect_effects. +Running this workflow is required to fully reproduce analysis in this directory. + +Code for the following analyses are in this directory: +- `dc_tapseq_indirect_effects_main_analyses.Rmd`: Main analyses shown in Figures 4 and S3 +- `dc_tapseq_indirect_effects_supplementary_analyses.Rmd`: Supplementary analyses and figures for modeling direct and indirect effect rates shown in Supplementary Figure 15 +- `dc_tapseq_indirect_effect_sizes.Rmd`: Effect size distributions of indirect effects shown in Supplementary Figure 16a +- `genome_wide_indirect_effects`: Comparison of indirect effects of DC-TAP-seq elements and target genes versus genome-wide elements and target genes shown in Supplementary Figures 16b,c. Scripts need to be run sequentially +- `bootstrapped_indirect_effects`: Bootstrap analysis to compute 95% confidence intervals for the modeled probability of direct effects + diff --git a/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/README.md b/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/README.md new file mode 100644 index 00000000..8b9a89c3 --- /dev/null +++ b/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/README.md @@ -0,0 +1,95 @@ +# Bootstrapped confidence intervals for the probability of direct effect + +This small workflow adds **per-pair confidence intervals** to the "probability +of direct effect" reported for the DC-TAP-seq results. It responds to the +reviewer comment that the framework "provides no quantification of per-pair +uncertainty." + +## What it does + +The published estimate (from the upstream workflow +`distal_regulation_paper/CRISPR_indirect_effects_FDR_check`) is: + +1. **indirect rate** per dataset = fraction of significant *trans* (other-chromosome) + perturbation–gene pairs; +2. **cis rate** per distance-to-TSS bin (50 kb) = fraction of significant cis pairs; +3. **direct rate** = `max(cis_rate − indirect_rate, 0)`, averaged across the 7 + datasets, fit as a power law `direct_rate = a · dist^b` (`lm(log ~ log)`); +4. **per-pair directness** = `direct_rate(dist) / (direct_rate(dist) + indirect_rate)`, + reported for positive and negative effects. + +This workflow re-implements that chain and wraps it in a **non-parametric cluster +(block) bootstrap that resamples perturbations** (it does **not** re-run +sceptre/MAST — it reuses the existing differential-expression results). Each of +`B = 1000` replicates resamples, within every dataset, its perturbations with +replacement (a perturbation's cis and trans pairs move together), recomputes the +rates, re-fits the power law, and recomputes each pair's directness. The 95% +confidence interval is the percentile interval across replicates. + +Resampling whole perturbations (rather than individual pairs) preserves the +within-perturbation correlation of trans effects and probes sensitivity to the +composition of the targeted panel. + +## Inputs (read-only, from upstream workflows) + +All 7 datasets that feed the fitted model are bootstrapped uniformly: + +- 6 "held-out" datasets (sceptre pipeline), per-dataset under + `.../CRISPR_indirect_effects_FDR_check/results//`: + `results_run_discovery_analysis.rds` (trans DE) and + `annotated_cis_results.tsv.gz` (cis DE; comma-separated despite the extension). +- 1 "training" dataset, Gasperini2019 (MAST pipeline) under + `.../ENCODE_CRISPR_data/results/Gasperini2019/`: + `diff_expr/output_MAST_perCRE.tsv.gz` (cis), + `trans_effects/output_trans_effects_MAST_perCRE.tsv.gz` (trans), and the + ENCODE-format `results/ENCODE/ENCODE_Gasperini2019_0.13gStd_MAST_perCRE_GRCh38.tsv.gz` + (for the ValidConnection enhancer filter). +- Shared: the gene-universe TSS BED (`CollapsedGeneBounds.hg38.TSS500bp.bed`) and + the DC-TAP results file to annotate + (`Final_DC_TAP_Seq_Results_..._fdr20_direct_effects.tsv`). +- For validation only: the published `direct_effects_model.rds` and + `trans_positive_hit_rates.tsv`. + +Paths are set in the config block at the top of `bootstrap_directness_CIs.R`. + +## How to run + +R (>=4.2.0) with the sceptre and tidyverse packages needs to be available. + +```bash +WD=$(pwd) # this directory + +# 1) bootstrap (≈ a few minutes on 8 cores; SLURM_CPUS_PER_TASK enables parallelism) +srun -p normal -n1 -c8 --mem=24G -t 00:45:00 --chdir="$WD" \ + env -u R_HOME BOOT_B=1000 SLURM_CPUS_PER_TASK=8 "Rscript" bootstrap_directness_CIs.R + +# 2) render the report +srun -p dev -n1 -c2 --mem=16G -t 00:20:00 --chdir="$WD" \ + env -u R_HOME "Rscript" -e 'rmarkdown::render("analyze_bootstrapped_directness.Rmd")' +``` + +Overridable environment variables: `BOOT_B` (replicates, default 1000), +`BOOT_SEED` (default 20250716), `SLURM_CPUS_PER_TASK` (cores). + +The script asserts that its re-implementation reproduces the published power-law +coefficients and per-dataset indirect rates before running the bootstrap; it +stops with an error if they do not match. + +## Outputs (`results/`) + +- `bootstrap_powerlaw_coefficients.tsv` — per replicate (`replicate` 0 = point + estimate): `type, intercept, slope, a, b`. +- `bootstrap_indirect_rates.tsv` — per replicate: `dataset, type, indirect_rate`. +- `directness_probability_CI_by_distance.tsv` — smooth curves for the ribbon plot: + `dist_to_tss, cell_type, type, point, median, lower, upper`. +- `Final_DC_TAP_Seq_Results_..._fdr20_direct_effects_with_CIs.tsv` — the DC-TAP + results file plus `direct_vs_indirect_{positive,negative}_{median,lower,upper}` + for every pair. +- `bootstrap_summary.rds` — compact object consumed by the report. +- `analyze_bootstrapped_directness.html` — report creating supplementary figure. + +## Files + +- `bootstrap_directness_CIs.R` — computation (loaders → normalized table → + cluster bootstrap → CIs). +- `analyze_bootstrapped_directness.Rmd` — report/plots. diff --git a/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/analyze_bootstrapped_directness.Rmd b/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/analyze_bootstrapped_directness.Rmd new file mode 100644 index 00000000..f4b1c334 --- /dev/null +++ b/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/analyze_bootstrapped_directness.Rmd @@ -0,0 +1,267 @@ +--- +title: "Confidence intervals for the per-pair probability of direct effect" +author: "DC-TAP paper -- reviewer response (per-pair uncertainty)" +date: "`r Sys.Date()`" +output: + html_document: + toc: true + toc_float: true + code_folding: hide + df_print: paged +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, + fig.width = 9, fig.height = 6, dpi = 100) +``` + +```{r attachPackages, message=FALSE, warning=FALSE} +library(tidyverse) +``` + +```{r loadData} +# load all required input data +res_dir <- file.path(getwd(), "results") +S <- readRDS(file.path(res_dir, "bootstrap_summary.rds")) +curves <- S$curves +coef_df <- S$coef +indirect <- S$indirect +ci_lab <- paste0(round(100 * S$ci_level), "%") + +# reformat curve data for plotting (done once, reused by all curve plots) +curves <- curves %>% + mutate(type = factor(type, levels = c("negative", "positive")), + cell_type = factor(cell_type)) + +# per-pair DC-TAP results + CIs, used in the per-pair analyses +dctap <- read_tsv(S$dctap_ci_file, show_col_types = FALSE) +``` + +## Motivation + +The "probability of direct effect" main analysis reports a per-pair point estimate but **no per-pair uncertainty**. Here we add confidence intervals by **bootstrapping the differential-expression +results** (not re-running sceptre/MAST). + +The published estimate is reconstructed exactly (see validation in the run log): an indirect-effect rate is estimated per dataset from other-chromosome (*trans*) perturbations, a distance-decaying +direct-effect rate is fit as a power law `direct_rate = a * dist^b`, and per pair: + +$$p_{\text{direct}}(d) = \frac{r_{\text{direct}}(d)}{r_{\text{direct}}(d) + r_{\text{indirect}}}$$ + +**Uncertainty model.** We use a non-parametric **cluster (block) bootstrap that resamples +perturbations with replacement** within each of the 7 datasets that feed the model +(`B = `r S$B``, seed `r S$seed`). Resampling whole perturbations (rather than individual pairs) +preserves the fact that a single perturbation could drives many *trans* pairs, and it probes +sensitivity to the composition of the candidate element selection. Each replicate recomputes the +direct/indirect rates, re-fits the power law, and recomputes the probability of direct effects. + +## Directness probability vs. distance, with `r ci_lab` CI + +```{r, fig.height=4, fig.width=8} +# plot probability of direct effects as function of distance to TSS including confidence intervals. +# the observational support behind these curves is shown separately below. +ggplot(curves, aes(x = dist_to_tss, group = interaction(cell_type, type))) + + facet_wrap(~cell_type, scales = "free") + + geom_ribbon(aes(ymin = lower, ymax = upper, fill = cell_type), alpha = 0.15) + + geom_line(aes(y = point, color = cell_type, lty = type), linewidth = 0.9) + + labs(x = "Distance to TSS", y = "P(direct effect)", + title = "Probability of direct effect with bootstrap CI", + subtitle = paste0("Solid = published point estimate; ribbon = ", ci_lab, " CI"), + fill = "Cell type", color = "Cell type", lty = "Direction of effect") + + scale_color_manual(values = c(K562 = "red", WTC11 = "blue")) + + scale_fill_manual(values = c(K562 = "red", WTC11 = "blue")) + + scale_x_continuous(labels = function(x) paste0(x/1e3, "kb")) + + scale_y_continuous(limits = c(0, 1)) + + theme_classic(base_size = 12) + + theme(strip.background = element_blank()) + +ggsave(filename = "../plots/prob_direct_effect_with_CI.pdf", height = 4, width = 8, create.dir = TRUE) +``` + +Include inferred median rate and show Distance to TSS in log10 space: +```{r ribbon} +ggplot(curves, aes(dist_to_tss)) + + geom_ribbon(aes(ymin = lower, ymax = upper, fill = cell_type), alpha = 0.25) + + geom_line(aes(y = point, color = cell_type), linewidth = 0.9) + + geom_line(aes(y = median, color = cell_type), linetype = "22", linewidth = 0.5) + + facet_grid(type ~ cell_type) + + scale_x_log10(labels = function(x) paste0(x/1e3, "kb")) + + labs(x = "Distance to TSS", y = "P(direct effect)", + title = "Probability of direct effect with bootstrap CI", + subtitle = paste0("Solid = published point estimate; dashed = bootstrap median; ribbon = ", + ci_lab, " CI"), + fill = "Cell type", color = "Cell type") + + theme_bw(base_size = 12) +``` + +The ribbon widens with distance to TSS -- where the direct-effect rate approaches +the indirect baseline, the probability that a hit is direct becomes most +uncertain. The published point estimate lies within the interval throughout. + +## Support behind the estimate: cis pairs and perturbations vs. distance + +Rather than annotate the curves with a rug (which could only show pairs from the +plotted cell types), we show the full observational support here: the cis pairs +and **distinct perturbations** that feed the direct-effect rate, across all +datasets used in the fit. Because the bootstrap resamples *perturbations* (not +pairs), the perturbation count is the unit of independent information -- and it is +what governs how stable each fitted rate is at a given distance. + +```{r support-counts, fig.height=6, fig.width=10} +cis_support <- read_tsv(file.path(res_dir, "cis_support_by_distance.tsv"), + show_col_types = FALSE) + +# bin distance the same way the direct-effect rate is computed (50 kb bins, <= 1 Mb) +bin_w <- 5e4; max_d <- 1e6 +support_binned <- cis_support %>% + filter(dist_to_tss <= max_d) %>% + mutate(bin_mid = (floor(dist_to_tss / bin_w) + 0.5) * bin_w) %>% + group_by(dataset, bin_mid) %>% + summarize(pairs = n(), perturbations = n_distinct(perturbation), .groups = "drop") %>% + pivot_longer(c(pairs, perturbations), names_to = "unit", values_to = "count") + +ggplot(support_binned, aes(bin_mid, count, color = unit)) + + geom_line(linewidth = 0.7) + + geom_point(size = 0.9) + + facet_wrap(~ dataset, scales = "free_y") + + scale_x_continuous(labels = function(x) paste0(x / 1e3, "kb")) + + labs(x = "Distance to TSS", y = "Count (cis)", color = NULL, + title = "Cis support behind the direct-effect rate", + subtitle = "Tested cis pairs and distinct perturbations per 50 kb bin (<= 1 Mb), per dataset") + + theme_bw(base_size = 12) +``` + +Breaking the **significant** hits down by effect direction shows directly why the +positive-effect CI blows up near the TSS: there are very few significant positive +hits at short distances, so that end of the positive direct-rate fit rests on a +handful of perturbations. + +```{r support-direction, fig.height=6, fig.width=10} +sig_binned <- cis_support %>% + filter(dist_to_tss <= max_d, significant) %>% + mutate(bin_mid = (floor(dist_to_tss / bin_w) + 0.5) * bin_w, + direction = if_else(regulated_positive, "positive", "negative")) %>% + group_by(dataset, bin_mid, direction) %>% + summarize(hits = n(), perturbations = n_distinct(perturbation), .groups = "drop") + +ggplot(sig_binned, aes(bin_mid, hits, color = direction)) + + geom_line(linewidth = 0.7) + + geom_point(size = 0.9) + + facet_wrap(~ dataset, scales = "free_y") + + scale_x_continuous(labels = function(x) paste0(x / 1e3, "kb")) + + labs(x = "Distance to TSS", y = "Significant cis hits", color = "Direction of effect", + title = "Significant cis hits by effect direction vs. distance", + subtitle = "Few positive hits near the TSS -> unstable positive direct-rate fit -> wide CI") + + theme_bw(base_size = 12) +``` + +## Power-law coefficients: bootstrap distribution + +```{r coefs} +coef_ci <- coef_df %>% + filter(replicate >= 1) %>% + group_by(type) %>% + summarize(a_median = median(a), a_lo = quantile(a, .025), a_hi = quantile(a, .975), + b_median = median(b), b_lo = quantile(b, .025), b_hi = quantile(b, .975), + .groups = "drop") %>% + left_join(coef_df %>% filter(replicate == 0) %>% select(type, a_point = a, b_point = b), + by = "type") +knitr::kable(coef_ci, digits = 4, + caption = paste0("Power-law coefficients (direct_rate = a * dist^b): point estimate and ", + ci_lab, " bootstrap CI")) + +ggplot(filter(coef_df, replicate >= 1), aes(b)) + + geom_histogram(bins = 40, fill = "grey70", color = "white") + + geom_vline(data = filter(coef_df, replicate == 0), aes(xintercept = b), + color = "firebrick", linewidth = 0.8) + + facet_wrap(~ type, scales = "free") + + labs(x = "power-law exponent b", y = "bootstrap replicates", + title = "Distribution of the distance-decay exponent", + subtitle = "red = point estimate") + + theme_bw(base_size = 12) +``` + +## Indirect (trans) rate uncertainty per dataset + +```{r indirect} +ind_ci <- indirect %>% + filter(type == "all") %>% + group_by(dataset) %>% + summarize(median = median(indirect_rate), + lo = quantile(indirect_rate, .025), hi = quantile(indirect_rate, .975), + .groups = "drop") %>% + arrange(median) +ggplot(ind_ci, aes(reorder(dataset, median), median)) + + geom_pointrange(aes(ymin = lo, ymax = hi)) + + coord_flip() + + labs(x = NULL, y = "indirect (trans) hit rate", + title = paste0("Indirect-rate ", ci_lab, " CI per dataset")) + + theme_bw(base_size = 12) +``` + +## Per-pair confidence intervals (DC-TAP results) + +```{r perpair} +# `dctap` (per-pair results + CIs) is loaded once in the 'support' chunk above + +# tidy the per-pair directness + CI into long form for the two effect types +long <- bind_rows(lapply(S$dctap_types, function(ty) { + tibble(cell_type = dctap$cell_type, + distance = abs(dctap$distance_to_gencode_gene_TSS), + significant = dctap$significant, + type = ty, + point = dctap[[paste0("direct_vs_indirect_", ty)]], + lower = dctap[[paste0("direct_vs_indirect_", ty, "_lower")]], + upper = dctap[[paste0("direct_vs_indirect_", ty, "_upper")]]) +})) %>% filter(!is.na(point)) + +long <- mutate(long, ci_width = upper - lower) + +knitr::kable( + long %>% group_by(type, cell_type) %>% + summarize(n = n(), median_width = median(ci_width, na.rm = TRUE), + median_point = median(point, na.rm = TRUE), .groups = "drop"), + digits = 3, caption = paste0("Per-pair ", ci_lab, " CI width summary")) +``` + +### CI width vs. distance + +```{r width-dist} +ggplot(long, aes(distance, ci_width, color = cell_type)) + + geom_point(alpha = 0.15, size = 0.6) + + facet_grid(type ~ cell_type) + + scale_x_log10(labels = function(x) paste0(x/1e3, "kb")) + + labs(x = "Distance to TSS", y = paste0(ci_lab, " CI width"), + title = "Per-pair CI width grows with distance to TSS") + + guides(color = "none") + theme_bw(base_size = 12) +``` + +### Example: significant pairs with directness + CI + +```{r examples} +ex <- long %>% + filter(significant == TRUE, type == "positive") %>% + arrange(distance) %>% + transmute(cell_type, distance_kb = round(distance/1e3, 1), + `P(direct)` = round(point, 3), + CI = sprintf("[%.3f, %.3f]", lower, upper)) +knitr::kable(utils::head(ex, 20), + caption = "First 20 significant pairs (positive effects) by distance, with directness CI") +``` + +## Summary for the response to reviewers + +- We now provide a **`r ci_lab` confidence interval for every pair's probability + of direct effect**, obtained by a perturbation-level cluster bootstrap of the + differential-expression results (no re-running of sceptre/MAST). +- The bootstrap **reproduces the published point estimate exactly** (power-law + coefficients and per-dataset indirect rates match to numerical tolerance; see + the run log) and adds uncertainty around it. +- Because whole perturbations are resampled, the interval reflects both sampling + noise in the hit rates **and** sensitivity to the composition of the targeted + panel -- the CI is widest at large distances to TSS, where the direct-effect + rate approaches the indirect baseline and directness is inherently least + certain. + +```{r session} +sessionInfo() +``` diff --git a/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/bootstrap_directness_CIs.R b/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/bootstrap_directness_CIs.R new file mode 100644 index 00000000..077db88a --- /dev/null +++ b/interpretation_analysis/indirect_effects/bootstrapped_indirect_effects/bootstrap_directness_CIs.R @@ -0,0 +1,511 @@ +## ================================================================================================= +## Confidence intervals for the per-pair "probability of direct effect" +## +## Addresses reviewer point (2): the published "probability of direct effect" +## framework reports a per-pair point estimate but no uncertainty. This script +## adds confidence intervals via a NON-PARAMETRIC CLUSTER (block) BOOTSTRAP that +## resamples perturbations (not by re-running sceptre/MAST). +## +## It faithfully re-implements the published estimation chain +## (see the upstream workflow CRISPR_indirect_effects_FDR_check): +## 1. indirect rate (per dataset) = fraction of significant trans pairs +## 2. cis rate (per distance-to-TSS bin) = fraction of significant cis pairs +## 3. direct rate = max(cis_rate - indirect_rate, 0), averaged across datasets, +## then a power law direct_rate = a * dist^b is fit by lm(log ~ log) +## 4. per-pair directness = direct_rate(dist) / (direct_rate(dist) + indirect_rate) +## +## The bootstrap resamples each dataset's perturbations with replacement (a +## perturbation's cis and trans pairs move together), recomputes 1-4, and records +## the power-law coefficients and indirect rates for every replicate. Percentiles +## of the resulting directness distribution give the confidence intervals. +## +## Run on a compute node (NOT the login node), e.g.: +## Rscript bootstrap_directness_CIs.R +## with the `sceptre` conda env (R 4.2.0 + tidyverse). See README.md. +## ================================================================================================= + +suppressPackageStartupMessages({ + library(readr) + library(dplyr) + library(tidyr) +}) + +## Define inputs and parameters -------------------------------------------------------------------- + +# upstream source workflows that hold the differential-expression (DE) results +SRC <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/CRISPR_indirect_effects_FDR_check" +ENC <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/ENCODE_CRISPR_data" + +# shared inputs +GENE_UNIV_FILE <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/CRISPR_benchmarks/resources/genome_annotations/CollapsedGeneBounds.hg38.TSS500bp.bed" +DCTAP_FILE <- file.path(SRC, "results/annotated_crispr_data/Final_DC_TAP_Seq_Results_w_Chromatin_Categories_on_resized_and_merged_elements_fdr20_direct_effects.tsv") + +# full direct effects models and indirect effect rates from DC-TAP-seq paper, used to validate that +# bootstrapping results match previous results +PUB_MODEL_FILE <- file.path(SRC, "results/direct_effect_models/direct_effects_model.rds") +PUB_TRANS_FILE <- file.path(SRC, "results/trans_positive_hit_rates.tsv") + +# the 6 "held-out" datasets processed by the sceptre pipeline +SCEPTRE_DATASETS <- c("K562_DC_TAPseq", "WTC11_DC_TAPseq", "Klann2021", + "Xie2019", "Morris2023_v1", "Morris2023_v2") + +# the 1 "training" dataset (Gasperini et al., 2019) processed by the (separate) MAST pipeline +GAS <- list( + dataset = "Gasperini2019", + cis = file.path(ENC, "results/Gasperini2019/diff_expr/output_MAST_perCRE.tsv.gz"), + trans = file.path(ENC, "results/Gasperini2019/trans_effects/output_trans_effects_MAST_perCRE.tsv.gz"), + encode = file.path(ENC, "results/ENCODE/ENCODE_Gasperini2019_0.13gStd_MAST_perCRE_GRCh38.tsv.gz") +) + +# DC-TAP cell type -> dataset whose indirect rate applies (mirrors dctap script) +CELLTYPE_TO_DATASET <- c(K562 = "K562_DC_TAPseq", WTC11 = "WTC11_DC_TAPseq") + +# parameters (bin_size / max_distance match the upstream workflow) +B <- as.integer(Sys.getenv("BOOT_B", "1000")) # bootstrap replicates +SEED <- as.integer(Sys.getenv("BOOT_SEED", "20250716")) +BIN_SIZE <- 5e4 +MAX_DIST <- 1e6 +CI_LEVEL <- 0.95 +CORES <- as.integer(Sys.getenv("SLURM_CPUS_PER_TASK", "1")) +DIST_COL <- "distance_to_gencode_gene_TSS" # DC-TAP distance column + +# effect types: "all" (=significant), "positive", "negative". DC-TAP directness +# is reported for positive & negative (matches predict_direct_vs_indirect_effect_dctap.R) +RATE_TYPES <- c("all", "positive", "negative") +DCTAP_TYPES <- c("positive", "negative") + +# create output directory if needed +OUTDIR <- file.path(getwd(), "results") +dir.create(OUTDIR, showWarnings = FALSE, recursive = TRUE) + +# lower, median and upper boundaries of confidence intervals +ci_probs <- c((1 - CI_LEVEL) / 2, 0.5, 1 - (1 - CI_LEVEL) / 2) # lower, median, upper + +set.seed(SEED) +message("Bootstrap CIs for per-pair probability of direct effects | B = ", B, ", seed = ", SEED, + ", cores = ", CORES) + +## Section 1: loaders -> one normalized per-pair table --------------------------------------------- +## Normalized columns: dataset, perturbation, arm ("cis"/"trans"), +## significant, regulated_positive, regulated_negative, dist_to_tss (cis only) +## Significance labels are computed once (frozen); the bootstrap resamples +## perturbations, not the significance threshold. + +# gene universe (Ensembl IDs) used to filter the sceptre datasets +gene_univ <- read_tsv(GENE_UNIV_FILE, show_col_types = FALSE) + +# function to load annotated cis-analysis results and sceptre trans-testing results for one dataset +load_sceptre_dataset <- function(input_dir, dataset) { + cis_file <- file.path(input_dir, "results", dataset, "annotated_cis_results.tsv.gz") + trans_file <- file.path(input_dir, "results", dataset, "results_run_discovery_analysis.rds") + + # NOTE: annotated cis file is comma-separated despite the .tsv.gz extension + cis <- read_csv(cis_file, show_col_types = FALSE) + trans <- as.data.frame(readRDS(trans_file)) + + # cis significance is precomputed (20% FDR); derive the p-value threshold and + # apply it to trans, exactly as calculate_positive_hit_rates.R does + pval_threshold_cis <- max(cis$p_value[cis$significant == TRUE], na.rm = TRUE) + + trans <- trans %>% + mutate(significant = p_value <= pval_threshold_cis, + regulated_negative = significant & log_2_fold_change < 0, + regulated_positive = significant & log_2_fold_change >= 0) + cis <- cis %>% + mutate(regulated_negative = significant & log_2_fold_change < 0, + regulated_positive = significant & log_2_fold_change >= 0) + + # valid distal elements (from cis) define which trans perturbations are valid + valid_elements <- unique(cis$grna_target[cis$valid_element == TRUE]) + trans <- mutate(trans, valid_connection = grna_target %in% valid_elements) + + # filter cis and trans results based on whether they are valid connection for the analysis + cis <- filter(cis, valid_connection == TRUE, pass_qc == TRUE, + response_id %in% gene_univ$Ensembl_ID) + trans <- filter(trans, valid_connection == TRUE, pass_qc == TRUE, + response_id %in% gene_univ$Ensembl_ID) + + # create tables with relevant data for bootstrapping analysis + cis_norm <- tibble(dataset = dataset, perturbation = cis$grna_target, arm = "cis", + significant = cis$significant, + regulated_positive = cis$regulated_positive, + regulated_negative = cis$regulated_negative, + dist_to_tss = abs(cis$dist_to_tss)) + trans_norm <- tibble(dataset = dataset, perturbation = trans$grna_target, arm = "trans", + significant = trans$significant, + regulated_positive = trans$regulated_positive, + regulated_negative = trans$regulated_negative, + dist_to_tss = NA_real_) + output <- bind_rows(cis_norm, trans_norm) + + return(output) + +} + +# function to load cis and trans results for Gasperini et al., 2019 dataset +load_gasperini_dataset <- function(g) { + cis <- read_tsv(g$cis, show_col_types = FALSE) + trans <- read_tsv(g$trans, show_col_types = FALSE) # already contains re-calculated significance + encode <- read_tsv(g$encode, show_col_types = FALSE) # contains ValidConnection column + + # get gene symbols for all genes + genes <- encode %>% + select(gene = measuredEnsemblID, gene_symbol = measuredGeneSymbol) %>% + distinct() + + # add gene symbols and ValidConnection lookups from the ENCODE-format table to cis results + cis <- cis %>% + left_join(genes, by = "gene") %>% + mutate(name = paste0(gene_symbol, "|", pert_chr, ":", pert_start, "-", pert_end, ":.")) %>% + left_join(distinct(select(encode, name, ValidConnection)), by = "name") + + # extract list of valid candidate enhancers not overlapping TSSs based on ValidConnection column + enh_filter <- c("overlaps potential promoter", "TSS targeting guide(s)") + valid_enh <- encode %>% + filter(!ValidConnection %in% enh_filter) %>% + pull(PerturbationTargetID) + + # label pairs with invalid candidate enhancers + trans <- trans %>% + mutate(pert_id = paste0(pert_chr, ":", pert_start, "-", pert_end, ":."), + ValidConnection = if_else(pert_id %in% valid_enh, "TRUE", "Invalid enhancer")) + + # cis significance = adjusted p < 0.05 (MAST pipeline); drop NA effect rows + cis <- cis %>% + filter(ValidConnection == "TRUE") %>% + mutate(significant = pval_adj < 0.05, + regulated_negative = significant & logFC < 0, + regulated_positive = significant & logFC >= 0) %>% + filter(!is.na(regulated_negative), !is.na(regulated_positive)) + trans <- filter(trans, ValidConnection == "TRUE") + + # create tables with relevant data for bootstrapping analysis + cis_norm <- tibble(dataset = g$dataset, perturbation = cis$perturbation, arm = "cis", + significant = cis$significant, + regulated_positive = cis$regulated_positive, + regulated_negative = cis$regulated_negative, + dist_to_tss = abs(cis$dist_to_tss)) + trans_norm <- tibble(dataset = g$dataset, perturbation = trans$perturbation, arm = "trans", + significant = trans$significant, + regulated_positive = trans$regulated_positive, + regulated_negative = trans$regulated_negative, + dist_to_tss = NA_real_) + output <- bind_rows(cis_norm, trans_norm) + + return(output) +} + +message("Loading and pre-processing DE results for all datasets ...") +norm_list <- c( + lapply(SCEPTRE_DATASETS, load_sceptre_dataset, input_dir = SRC), + list(load_gasperini_dataset(GAS)) +) +norm <- bind_rows(norm_list) +ALL_DATASETS <- c(SCEPTRE_DATASETS, GAS$dataset) + +## Export cis support (pairs + perturbations per distance) for the report -------------------------- +## `norm` already holds, per cis pair: dataset, perturbation, distance, significance and effect +## direction. This is the observational support behind the direct-effect rate. Since the bootstrap +## resamples the `perturbation` clusters, both the pair and perturbation counts are informative. +cis_support <- norm %>% + filter(arm == "cis") %>% + transmute(dataset, perturbation, dist_to_tss, + significant, regulated_positive, regulated_negative) +write_tsv(cis_support, file.path(OUTDIR, "cis_support_by_distance.tsv")) + +# regenerate just the support table (no bootstrap) with: SUPPORT_ONLY=1 Rscript bootstrap_directness_CIs.R +if (nzchar(Sys.getenv("SUPPORT_ONLY"))) { + message("SUPPORT_ONLY set -- wrote cis_support_by_distance.tsv and exiting before the bootstrap.") + quit(save = "no") +} + +## Section 2: precompute per-dataset bootstrap structures ------------------------------------------ +## For each dataset: perturbation index (shared by cis+trans), fixed distance +## bins (from full cis data), and the logical hit vectors. The bootstrap then +## only reweights rows by per-perturbation draw counts (cluster bootstrap). + +# coerce a hit indicator to numeric with NA -> 0. This reproduces the upstream +# rate scripts, which count hits with sum(..., na.rm = TRUE) while keeping the +# full pair count n() in the denominator (NA = not a hit, still a tested pair). +na0 <- function(x) { x <- as.numeric(x); x[is.na(x)] <- 0; x } + +build_dataset_struct <- function(d) { + sub <- filter(norm, dataset == d) + cis <- filter(sub, arm == "cis") + trans <- filter(sub, arm == "trans") + + perts <- unique(c(cis$perturbation, trans$perturbation)) # should be the same for cis and trans + + # fixed distance bins for this dataset (from full cis data), reused every rep + max_dist <- ceiling(max(cis$dist_to_tss, na.rm = TRUE) / 1e6) * 1e6 + breaks <- seq(0, max_dist, by = BIN_SIZE) + cis_bin <- cut(cis$dist_to_tss, breaks = breaks, include.lowest = TRUE) + + # construct object containing relevant information from the given dataset + output <- list( + dataset = d, + n_pert = length(perts), + cis = list(pert_idx = match(cis$perturbation, perts), + bin = cis_bin, + dist = cis$dist_to_tss, + sig = na0(cis$significant), + pos = na0(cis$regulated_positive), + neg = na0(cis$regulated_negative)), + trans = list(pert_idx = match(trans$perturbation, perts), + sig = na0(trans$significant), + pos = na0(trans$regulated_positive), + neg = na0(trans$regulated_negative)) + ) + + return(output) + +} + +ds_structs <- lapply(ALL_DATASETS, build_dataset_struct) +names(ds_structs) <- ALL_DATASETS + +## Section 2b: rate + power-law functions ---------------------------------------------------------- + +# per-dataset direct rates (per bin) + trans/indirect rates, given perturbation +# draw counts (weights). weights = rep(1, n_pert) reproduces the point estimate. +dataset_rates <- function(ds, weights) { + + # calculate how many times each E-G pair occurs in the bootstrap sample based on weights + wc <- weights[ds$cis$pert_idx] + wt <- weights[ds$trans$pert_idx] + + # cis: weighted sums per distance bin + M <- rowsum(cbind(w = wc, sig = wc * ds$cis$sig, pos = wc * ds$cis$pos, + neg = wc * ds$cis$neg, wd = wc * ds$cis$dist), + group = ds$cis$bin, reorder = FALSE) + tot <- M[, "w"] + keep <- tot > 0 + M <- M[keep, , drop = FALSE]; tot <- tot[keep] + + # cis hit rate across distance bins per rate type + cis_rate <- data.frame( + dataset = ds$dataset, + dist_bin = rownames(M), + mean_dist = M[, "wd"] / tot, + cis_all = M[, "sig"] / tot, + cis_pos = M[, "pos"] / tot, + cis_neg = M[, "neg"] / tot, + stringsAsFactors = FALSE, row.names = NULL + ) + + # distance independent trans hit rate per rate type + tt <- sum(wt) + ind <- c(all = sum(wt * ds$trans$sig) / tt, + positive = sum(wt * ds$trans$pos) / tt, + negative = sum(wt * ds$trans$neg) / tt) + + # direct rate = max(cis - indirect, 0), per type + cis_rate$direct_all <- pmax(cis_rate$cis_all - ind["all"], 0) + cis_rate$direct_pos <- pmax(cis_rate$cis_pos - ind["positive"], 0) + cis_rate$direct_neg <- pmax(cis_rate$cis_neg - ind["negative"], 0) + + # create output containing boootstrapped direct and indirect hit rates + output <- list(cis = cis_rate, indirect = ind) + + return(output) + +} + +# fit the three power-law models from all datasets' per-bin direct rates; +# returns intercept (= log a) and slope (= b) for each type +fit_models <- function(rates_list) { + + # get cis hit rates within given maximum distance + cis <- do.call(rbind, lapply(rates_list, `[[`, "cis")) + cis <- cis[cis$mean_dist <= MAX_DIST, , drop = FALSE] + + # average across datasets within each distance bin + avg <- cis %>% + group_by(dist_bin) %>% + summarize(dist_to_tss = mean(mean_dist), + direct_all = mean(direct_all), + direct_pos = mean(direct_pos), + direct_neg = mean(direct_neg), + .groups = "drop") + + # fit power law model to bootstrapped direct effects rates + fit_one <- function(y) { + df <- data.frame(direct_rate = avg[[y]], dist_to_tss = avg$dist_to_tss) + df <- df[df$direct_rate > 0, , drop = FALSE] + co <- coef(lm(log(direct_rate) ~ log(dist_to_tss), data = df)) + c(intercept = unname(co[1]), slope = unname(co[2])) + } + + # create output + output <- list(all = fit_one("direct_all"), + positive = fit_one("direct_pos"), + negative = fit_one("direct_neg")) + + return(output) + +} + +# one bootstrap replicate: resample each dataset's perturbations with replacement +one_replicate <- function(seed_i) { + set.seed(seed_i) + rates_list <- lapply(ds_structs, function(ds) { + counts <- tabulate(sample.int(ds$n_pert, ds$n_pert, replace = TRUE), nbins = ds$n_pert) + dataset_rates(ds, counts) + }) + models <- fit_models(rates_list) + indirect <- sapply(rates_list, `[[`, "indirect") # types x datasets + colnames(indirect) <- ALL_DATASETS + list(models = models, indirect = indirect) +} + +## Section 3: point estimate + validation ---------------------------------------------------------- + +point_rates <- lapply(ds_structs, function(ds) dataset_rates(ds, rep(1, ds$n_pert))) +point_models <- fit_models(point_rates) +point_indirect <- sapply(point_rates, `[[`, "indirect"); colnames(point_indirect) <- ALL_DATASETS + +# validate power-law coefficients against the published model +pub_model <- readRDS(PUB_MODEL_FILE) +pub_coef <- sapply(c(all = "direct_rate_all", positive = "direct_rate_positive", + negative = "direct_rate_negative"), + function(nm) coef(pub_model[[nm]])) +message("\n== Validation: recomputed vs published power-law coefficients ==") +for (ty in RATE_TYPES) { + rc <- point_models[[ty]] + pc <- pub_coef[, ty] + message(sprintf(" %-9s intercept %.5f (pub %.5f) slope %.5f (pub %.5f)", + ty, rc["intercept"], pc[1], rc["slope"], pc[2])) + stopifnot(abs(rc["intercept"] - pc[1]) < 1e-3, abs(rc["slope"] - pc[2]) < 1e-3) +} + +# validate indirect (trans) rates against the published table +pub_trans <- read_tsv(PUB_TRANS_FILE, show_col_types = FALSE) +message("== Validation: recomputed vs published indirect rates (significant) ==") +for (d in ALL_DATASETS) { + pubv <- pub_trans$positive_rate_significant[pub_trans$dataset == d] + recv <- point_indirect["all", d] + message(sprintf(" %-16s recomputed %.6f published %.6f", d, recv, pubv)) + stopifnot(abs(recv - pubv) < 1e-4) +} +message("Validation passed: re-implementation reproduces the published point estimate.\n") + +## Section 4: bootstrap loop ----------------------------------------------------------------------- + +message("Running ", B, " bootstrap replicates ...") +seeds <- SEED + seq_len(B) +boot <- if (CORES > 1) { + parallel::mclapply(seeds, one_replicate, mc.cores = CORES) +} else { + lapply(seeds, one_replicate) +} + +# tidy the bootstrap draws +coef_rows <- lapply(seq_along(boot), function(i) { + m <- boot[[i]]$models + do.call(rbind, lapply(RATE_TYPES, function(ty) { + data.frame(replicate = i, type = ty, + intercept = m[[ty]]["intercept"], slope = m[[ty]]["slope"], + a = exp(m[[ty]]["intercept"]), b = m[[ty]]["slope"], + row.names = NULL) + })) +}) +coef_df <- do.call(rbind, coef_rows) + +indirect_rows <- lapply(seq_along(boot), function(i) { + ind <- boot[[i]]$indirect + do.call(rbind, lapply(ALL_DATASETS, function(d) + data.frame(replicate = i, dataset = d, type = RATE_TYPES, + indirect_rate = ind[RATE_TYPES, d], row.names = NULL))) +}) +indirect_df <- do.call(rbind, indirect_rows) + +# add the point estimate as replicate 0 for reference +point_coef <- do.call(rbind, lapply(RATE_TYPES, function(ty) + data.frame(replicate = 0L, type = ty, + intercept = point_models[[ty]]["intercept"], slope = point_models[[ty]]["slope"], + a = exp(point_models[[ty]]["intercept"]), b = point_models[[ty]]["slope"], + row.names = NULL))) +coef_df <- rbind(point_coef, coef_df) + +write_tsv(coef_df, file.path(OUTDIR, "bootstrap_powerlaw_coefficients.tsv")) +write_tsv(indirect_df, file.path(OUTDIR, "bootstrap_indirect_rates.tsv")) + +## Section 5: propagate to directness CIs ---------------------------------------------------------- + +# helper: directness probability across all B replicates for a vector of +# distances, one cell type, one effect type -> returns a B x length(dist) matrix +directness_matrix <- function(dist, type, dataset) { + inter <- coef_df$intercept[coef_df$replicate >= 1 & coef_df$type == type] + slope <- coef_df$slope[coef_df$replicate >= 1 & coef_df$type == type] + ind <- indirect_df$indirect_rate[indirect_df$dataset == dataset & indirect_df$type == type] + logd <- log(dist) + # B x n: direct rate (capped at 1), then probability + direct <- pmin(exp(outer(inter, rep(1, length(dist))) + outer(slope, logd)), 1) + direct / (direct + ind) # ind (length B) recycles down each column +} + +# point directness for a distance vector (published coefficients + published indirect) +point_directness <- function(dist, type, dataset) { + co <- point_models[[type]] + ind <- point_indirect[type, dataset] + direct <- pmin(exp(co["intercept"] + co["slope"] * log(dist)), 1) + unname(direct / (direct + ind)) +} + +## 5a. smooth CI curves over a distance grid (for the ribbon plot) +grid <- 10 ^ seq(log10(1e3), log10(MAX_DIST), length.out = 200) +curve_rows <- list(); k <- 1 +for (ct in names(CELLTYPE_TO_DATASET)) { + d <- CELLTYPE_TO_DATASET[[ct]] + for (ty in DCTAP_TYPES) { + mat <- directness_matrix(grid, ty, d) + qs <- apply(mat, 2, quantile, probs = ci_probs, na.rm = TRUE) + curve_rows[[k]] <- data.frame( + dist_to_tss = grid, cell_type = ct, type = ty, + point = point_directness(grid, ty, d), + median = qs[2, ], lower = qs[1, ], upper = qs[3, ], row.names = NULL) + k <- k + 1 + } +} +curve_df <- do.call(rbind, curve_rows) +write_tsv(curve_df, file.path(OUTDIR, "directness_probability_CI_by_distance.tsv")) + +## 5b. per-pair CIs for the DC-TAP results file +message("Computing per-pair directness CIs for the DC-TAP results ...") +dctap <- read_tsv(DCTAP_FILE, show_col_types = FALSE) +dctap$.dist <- pmax(abs(dctap[[DIST_COL]]), 1) # guard log(0)/negatives +dctap$.dataset <- unname(CELLTYPE_TO_DATASET[dctap$cell_type]) + +for (ty in DCTAP_TYPES) { + lo <- md <- up <- rep(NA_real_, nrow(dctap)) + for (ct in names(CELLTYPE_TO_DATASET)) { + idx <- which(dctap$cell_type == ct) + if (!length(idx)) next + mat <- directness_matrix(dctap$.dist[idx], ty, CELLTYPE_TO_DATASET[[ct]]) + qs <- apply(mat, 2, quantile, probs = ci_probs, na.rm = TRUE) + lo[idx] <- qs[1, ]; md[idx] <- qs[2, ]; up[idx] <- qs[3, ] + } + dctap[[paste0("direct_vs_indirect_", ty, "_median")]] <- md + dctap[[paste0("direct_vs_indirect_", ty, "_lower")]] <- lo + dctap[[paste0("direct_vs_indirect_", ty, "_upper")]] <- up +} +dctap <- select(dctap, -c(.dist, .dataset)) + +out_dctap <- file.path(OUTDIR, sub("\\.tsv$", "_with_CIs.tsv", basename(DCTAP_FILE))) +write_tsv(dctap, out_dctap) + +## Save a compact object for the report + session info --------------------------------------------- +saveRDS(list(coef = coef_df, indirect = indirect_df, curves = curve_df, + point_models = point_models, point_indirect = point_indirect, + pub_coef = pub_coef, B = B, seed = SEED, ci_level = CI_LEVEL, + celltype_to_dataset = CELLTYPE_TO_DATASET, dctap_types = DCTAP_TYPES, + dctap_ci_file = out_dctap), + file.path(OUTDIR, "bootstrap_summary.rds")) + +message("\nDone. Outputs written to ", OUTDIR) +message(" - bootstrap_powerlaw_coefficients.tsv") +message(" - bootstrap_indirect_rates.tsv") +message(" - directness_probability_CI_by_distance.tsv") +message(" - ", basename(out_dctap)) +message(" - bootstrap_summary.rds") +print(sessionInfo()) diff --git a/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effect_sizes.Rmd b/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effect_sizes.Rmd new file mode 100644 index 00000000..d87c1f56 --- /dev/null +++ b/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effect_sizes.Rmd @@ -0,0 +1,253 @@ +--- +title: "DC-TAP-seq indirect effect sizes" +author: "Andreas R. Gschwind" +date: "July 23, 2026" +output: + html_document: + toc: true + toc_float: true + code_folding: show +--- + +```{r setupDocument, include=FALSE} +# set output html chunk options +knitr::opts_chunk$set(warning = FALSE, message = FALSE) + +# suppress dplyrs summarize info +options(dplyr.summarize.inform = FALSE) +``` + +```{r requiredPackages, warning=FALSE, message=FALSE} +library(tidyverse) +library(cowplot) +``` + +## Goal +Plot effect sizes of indirect effects in K562 and WTC11 DC-TAP-seq datasets. +```{r} +# path to indirect effects analysis output (adjust to reproduce) +indirect_effects_dir <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/CRISPR_indirect_effects_FDR_check/results" + +# path to Gasperini et al indirect effects analysis output (adjust to reproduce) +gasp_dir <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/ENCODE_CRISPR_data/results/Gasperini2019" +``` + +*** + +## Used data +Sceptre and MAST output used to calculated indirect effect rates across all datasets is used to +analyze effect sizes of detected significant indirect effects. +```{r} +# all output files for sceptre-based trans-acting effect tests +scpt_datasets <- c("K562_DC_TAPseq", "WTC11_DC_TAPseq", "Klann2021", "Xie2019", "Morris2023_v1", "Morris2023_v2") +scpt_files <- file.path(indirect_effects_dir, scpt_datasets, "results_run_discovery_analysis.rds") +names(scpt_files) <- scpt_datasets + +# load all sceptre result files +sceptre_trans_results <- lapply(scpt_files, FUN = readRDS) + +# all output files for sceptre-based cis analysis +cis_files <- file.path(indirect_effects_dir, scpt_datasets, "annotated_cis_results.tsv.gz") +names(cis_files) <- scpt_datasets + +# load all cis analysis results +sceptre_cis_results <- lapply(cis_files, FUN = read_csv, show_col_types = FALSE) + +# gasperini et al trans and cis result files +gasp_trans_file <- file.path(gasp_dir, "trans_effects/output_trans_effects_MAST_perCRE.tsv.gz") +#gasp_cis_file <- file.path(gasp_dir, "diff_expr/output_MAST_perCRE.tsv.gz") +gasp_encode_file <- file.path(gasp_dir, "../ENCODE/ENCODE_Gasperini2019_0.13gStd_MAST_perCRE_GRCh38.tsv.gz") + +# load gasperini files +gasp_trans_results <- read_tsv(gasp_trans_file, show_col_types = FALSE) +#gasp_cis_results <- read_tsv(gasp_cis_file, show_col_types = FALSE) +gasp_encode <- read_tsv(gasp_encode_file, show_col_types = FALSE) + +# load gene universe file +gene_univ <- read_tsv("/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/CRISPR_benchmarks/resources/genome_annotations/CollapsedGeneBounds.hg38.TSS500bp.bed", + show_col_types = FALSE) +``` + +*** + +## Extract effect sizes for sigificant indirect hits +Filter trans differential expression results using the same filters applied ot calculate indirect +effect rates and extract indirect effect sizes. + +### Sceptre based results +```{r} +# helper function to extract significant indirect effects for one dataset +get_indirect_effects <- function(trans_results, cis_results, gene_univ) { + + # get significance threshold for cis effects + pval_threshold_cis <- cis_results %>% + filter(significant == TRUE) %>% + pull(p_value) %>% + max() + + # add significance and regulated column to trans-effect results based on cis significance threshold + trans_results <- trans_results %>% + mutate(significant = p_value <= pval_threshold_cis, + regulated_negative = significant & log_2_fold_change < 0, + regulated_positive = significant & log_2_fold_change >= 0) + + # get all valid elements from cis results + valid_elements <- cis_results %>% + filter(valid_element == TRUE) %>% + pull(grna_target) %>% + unique() + + # add a valid connection column to trans-results based on whether perturbation is a valid element + trans_results <- mutate(trans_results, valid_connection = grna_target %in% valid_elements) + + # filter trans results for valid connections and only retain data on gene universe + trans_results_filt <- trans_results %>% + filter(valid_connection == TRUE, pass_qc == TRUE) %>% + filter(response_id %in% gene_univ$Ensembl_ID) + + # extract significantly regulated indirect effects + indirect_effects <- trans_results_filt %>% + filter(regulated_negative == TRUE | regulated_positive == TRUE) %>% + mutate(effect_direction = case_when( + regulated_negative & regulated_positive ~ "both (ERROR)", + regulated_negative ~ "negative", + regulated_positive ~ "positive", + TRUE ~ NA_character_ + )) + + return(indirect_effects) +} + +# helper function to convert log-fold changes to percentage changes +logFC_to_pctChange <- function(logFC, base = exp(1)) { + pct_change <- base ^ logFC - 1 + return(pct_change) +} +``` + +```{r} +# extract significant indirect effects for both dataset +indirect_effects <- mapply(FUN = get_indirect_effects, sceptre_trans_results, sceptre_cis_results, + MoreArgs = list(gene_univ = gene_univ), SIMPLIFY = FALSE) + +# combine into one table and convert effect sizes from log2FC to percent change in expression +indirect_effects <- indirect_effects %>% + bind_rows(.id = "dataset") %>% + mutate(pct_change = logFC_to_pctChange(log_2_fold_change, base = 2)) +``` + +### MAST based results (Gasperini dataset) +```{r} +# get valid enhancers for gasperini dataset based on ENCODE output +enh_filter <- c("overlaps potential promoter", "TSS targeting guide(s)") +gasp_valid_enh <- gasp_encode %>% + filter(!ValidConnection %in% enh_filter) %>% + pull(PerturbationTargetID) + +# add ValidConnection column to gasperini trans-results +gasp_trans_results <- gasp_trans_results %>% + mutate(pert_id = paste0(pert_chr, ":", pert_start, "-", pert_end, ":.")) %>% + mutate(ValidConnection = if_else(pert_id %in% gasp_valid_enh, true = "TRUE", + false = "Invalid enhancer")) + +# filter trans results based on ValidConnection +gasp_trans_results <- filter(gasp_trans_results, ValidConnection == "TRUE") + +# extract significantly regulated indirect effects +gasp_trans_results <- gasp_trans_results %>% + filter(regulated_negative == TRUE | regulated_positive == TRUE) %>% + mutate(effect_direction = case_when( + regulated_negative & regulated_positive ~ "both (ERROR)", + regulated_negative ~ "negative", + regulated_positive ~ "positive", + TRUE ~ NA_character_ + )) + +# convert effect sizes from log2FC to percent change in expression +gasp_trans_results <- gasp_trans_results %>% + mutate(pct_change = logFC_to_pctChange(logFC, base = exp(1))) +``` + +### Combine into one table +```{r} +# select relevant columns and combine from sceptre based trans results +indirect_effects <- indirect_effects %>% + select(dataset, perturbation = grna_target, gene = response_id, regulated_negative, + regulated_positive, effect_direction, pct_change) + +# select relevant columns and combine from gasperini trans results +gasp_trans_results <- gasp_trans_results %>% + mutate(dataset = "Gasperini2019") %>% + select(dataset, perturbation, gene, regulated_negative, regulated_positive, effect_direction, + pct_change) + +# combine into one table +indirect_effects <- bind_rows(indirect_effects, gasp_trans_results) +``` + +*** + +## Plot indirect effect sizes + +Group data for plots: 1 = K562 DC-TAP-seq, 2 = WTC11 DC-TAP-seq, 3 = Other datasets. +```{r} +# redefine datasets to group data for plots +indirect_effects <- indirect_effects %>% + mutate(dataset_plot = case_when( + dataset == "K562_DC_TAPseq" ~ "K562", + dataset == "WTC11_DC_TAPseq" ~ "WTC11", + TRUE ~ "Other datasets" + )) %>% + mutate(dataset_plot = factor(dataset_plot, levels = c("K562", "WTC11", "Other datasets"))) +``` + +```{r} +# label outliers in each plotted distribution for filtering to prevent long tails on violin plots +indirect_effects <- indirect_effects %>% + group_by(dataset_plot, effect_direction) %>% + mutate(outlier = !between(pct_change, + quantile(pct_change, 0.25) - 1.5 * IQR(pct_change), + quantile(pct_change, 0.75) + 1.5 * IQR(pct_change))) + +# calculate percentage of outlier points removed per distribution +indirect_effects <- indirect_effects %>% + mutate(total_pairs = n(), + outlier_pairs = sum(outlier), + pct_outliers = outlier_pairs / total_pairs) + +# filter out outlier points +indirect_effects_filt <- filter(indirect_effects, outlier == FALSE) + +# create labels with number of hits for plots +indirect_effects_filt <- indirect_effects_filt %>% + mutate(filt_pairs = n()) %>% + mutate(label = paste0(effect_direction, "\n(n = ", filt_pairs, " trans-pairs,\n", outlier_pairs, + " outliers removed)")) +``` + + +Plot effect sizes for negative and positive indirect effects as violin plots. +```{r, fig.height=4, fig.width=11} +# plot effect sizes of indirect effects as violin plots +ggplot(indirect_effects_filt, aes(x = label, y = pct_change, fill = effect_direction)) + + facet_wrap(~dataset_plot, drop = TRUE, scales = "free") + + geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5) + + geom_violin() + + geom_boxplot(color = "black", fill = NA, outlier.shape = NA, width = 0.2) + + labs(y = "CRISPR effect size\n(% change in expression)", x = "Indirect effect direction") + + scale_fill_manual(values = c("firebrick3", "steelblue")) + + scale_y_continuous(limits = range(indirect_effects_filt$pct_change), labels = scales::percent) + + theme_classic() + + theme(legend.position = "none", strip.background = element_blank()) + +# save plot to pdf +ggsave(file = "plots/indirect_effect_sizes_violin.pdf", height = 4, width = 11, create.dir = TRUE) +``` + +*** + +## Session information +This document was produced using following packages: +```{r sessionInfo} +sessionInfo() +``` diff --git a/interpretation_analysis/dc_tapseq_indirect_effects.Rmd b/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effects_main_analyses.Rmd similarity index 86% rename from interpretation_analysis/dc_tapseq_indirect_effects.Rmd rename to interpretation_analysis/indirect_effects/dc_tapseq_indirect_effects_main_analyses.Rmd index a19bfc68..351afac8 100644 --- a/interpretation_analysis/dc_tapseq_indirect_effects.Rmd +++ b/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effects_main_analyses.Rmd @@ -51,7 +51,7 @@ analysis workflow as input. Probabilities of direct effects for K562 and WTC11 D experiments calculated by this workflow can be found in Table S3. ```{r} # adjust path to indirect effects analysis output as needed to reproduce results -indirect_effects_dir <- "../../CRISPR_indirect_effects_FDR_check/results" +indirect_effects_dir <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/CRISPR_indirect_effects_FDR_check/results" # adjust path to table S3 as needed to reproduce results table_s3_file <- "./Table_S3_Final_DC_TAP_Seq_Results_w_Chromatin_Categories_on_resized_and_merged_elements_250908_append.tsv" @@ -181,42 +181,87 @@ as functions of distance to TSS Here we show the overall (direct + indirect) cis hit rate of observed significant DE-gene pairs for 50kb bins of distance to TSS. The distance to TSS independent indirect effects rates are shown as lines. -```{r cisVsIndirectEffectRate, fig.height=4, fig.width=8} + +```{r plotObservedRates} # reformat cis hit rates for plot cis_rates_plot <- cis_rates %>% select(dataset, mean_dist, positive_rate_negative, positive_rate_positive) %>% pivot_longer(cols = -c(dataset, mean_dist), names_to = "effect_direction", values_to = "positive_rate") %>% mutate(cell_type = sub("_DC_TAPseq", "", dataset)) %>% - mutate(effect_direction = str_to_title(sub("positive_rate_", "", effect_direction))) + mutate(effect_direction = str_to_title(sub("positive_rate_", "", effect_direction))) %>% + mutate(effect_direction = factor(effect_direction)) + +# helper function to plot observed hit rates as function of distance to TSS +plot_observed_rates <- function(observed_rates, ylim) { + ggplot(observed_rates, aes(x = mean_dist / 1e6, y = positive_rate, color = cell_type, + shape = effect_direction)) + + geom_hline(yintercept = 0, linetype = "dashed", color = "gray") + + geom_point() + + labs(title = paste(unique(observed_rates$effect_direction), "effects on expression"), + x = "Distance to TSS (Mb)", + y = "Rate of observed effects", color = "Cell type", shape = "Direction of effect") + + scale_color_manual(values = c(K562 = "#b00404", WTC11 = "#0202cc")) + + scale_shape_manual(values = c(Negative = 19, Positive = 1), drop = FALSE) + + scale_y_continuous(limits = ylim, labels = scales::percent) + + theme_classic() + + theme(legend.position = "top") +} + +# plot observed hit rates for negative and positive effects on expression as function of distance +ylim <- c(NA, max(cis_rates_plot$positive_rate)) +p_obs_neg <- plot_observed_rates(filter(cis_rates_plot, effect_direction == "Negative"), ylim = ylim) +p_obs_pos <- plot_observed_rates(filter(cis_rates_plot, effect_direction == "Positive"), ylim = ylim) +``` +```{r plotIndirectRates} # reformat trans hit rates for plot trans_rates_plot <- trans_rates %>% select(dataset, positive_rate_negative, positive_rate_positive) %>% pivot_longer(cols = -dataset, names_to = "effect_direction", values_to = "positive_rate") %>% mutate(cell_type = sub("_DC_TAPseq", "", dataset)) %>% - mutate(effect_direction = str_to_title(sub("positive_rate_", "", effect_direction))) + mutate(effect_direction = str_to_title(sub("positive_rate_", "", effect_direction))) -# plot cis and trans hit rates as function of distance to TSS -ggplot(cis_rates_plot, - aes(x = mean_dist / 1e6, y = positive_rate, color = cell_type, shape = effect_direction)) + - facet_wrap(~cell_type, scales = "free") + - geom_point() + - geom_hline(data = trans_rates_plot, aes(yintercept = positive_rate, linetype = effect_direction, - color = cell_type)) + - labs(title = "Positive hit rates cis and trans analyses", x = "Distance to TSS (Mb)", - y = "Positive hite rate", color = "Cell type", shape = "Effect direction", - linetype = "Effect direction") + - scale_color_manual(values = c(K562 = "#b00404", WTC11 = "#0202cc")) + - scale_shape_manual(values = c(Negative = 19, Positive = 1)) + - scale_y_continuous(limits = c(NA, max(cis_rates_plot$positive_rate)), labels = scales::percent) + - theme_classic() + - theme(strip.background = element_blank()) +# helper function to plot indirect effects rates +plot_indirect_rates <- function(indirect_rates, ylim) { + ggplot(indirect_rates, aes(x = "Other\nchromosomes\n(Rate of\nindirect effects)", + y = positive_rate, color = cell_type, shape = effect_direction)) + + geom_hline(yintercept = 0, linetype = "dashed", color = "gray") + + geom_point() + + scale_color_manual(values = c(K562 = "#b00404", WTC11 = "#0202cc")) + + scale_shape_manual(values = c(Negative = 19, Positive = 1)) + + scale_y_continuous(limits = ylim, labels = scales::percent) + + theme_classic() + + theme(legend.position = "none", axis.title.x = element_blank(), + axis.line.y = element_blank(), axis.title.y = element_blank(), + axis.text.y = element_blank(), axis.ticks.y = element_blank()) +} + +# plot indirect hit rates for negative and positive effects on expression +p_ind_neg <- plot_indirect_rates(filter(trans_rates_plot, effect_direction == "Negative"), ylim = ylim) +p_ind_pos <- plot_indirect_rates(filter(trans_rates_plot, effect_direction == "Positive"), ylim = ylim) +``` + +```{r assembleFigure, fig.height=4, fig.width=7} +# extract legend from one plot +legend <- get_legend(p_obs_neg) + +# arrange all plots +plots <- plot_grid( + p_obs_neg + theme(legend.position = "none"), + p_ind_neg, + p_obs_pos + theme(legend.position = "none"), + p_ind_pos, + nrow = 1, rel_widths = c(1, 0.1, 1, 0.1), align = "h" +) + +# combine legend and plots into one figure +plot_grid(legend, plots, nrow = 2, rel_heights = c(0.1, 1)) ``` ```{r savePlot1, echo=FALSE} # save plot to pdf -ggsave(file = "./plots/cis_and_trans_effect_rates.pdf", height = 4, width = 8) +ggsave(file = "plots/cis_and_trans_effect_rates.pdf", height = 4, width = 7, create.dir = TRUE) ``` #### Probability of direct effect @@ -262,7 +307,7 @@ ggplot(random_pairs_plot, aes(x = dist_to_tss / 1e3, y = direct_vs_indirect, col ```{r savePlot2, echo=FALSE} # save plot to pdf -ggsave(file = "./plots//probability_direct_vs_indirect.pdf", height = 4, width = 6) +ggsave(file = "plots/probability_direct_vs_indirect.pdf", height = 4, width = 6) ``` *** @@ -442,7 +487,7 @@ ggplot(pairs_p90_combined_plot, ```{r, savePlot3, echo=FALSE} # save plot to pdf -ggsave(file = "./plots//negative_vs_positive_effects_combined.pdf", height = 4, width = 3.5) +ggsave(file = "plots/negative_vs_positive_effects_combined.pdf", height = 4, width = 3.5) ``` Alternatively, we can also look at this by plotting the observed CRISPR effect size on target gene @@ -466,7 +511,7 @@ ggplot(random_pairs_hits, aes(x = cell_type, y = log_2_FC_effect_size, color = f ```{r savePlot4, echo=FALSE} # save plot to pdf -ggsave(file = "./plots//probability_direct_vs_effect_size.pdf", height = 4, width = 6) +ggsave(file = "plots/probability_direct_vs_effect_size.pdf", height = 4, width = 6) ``` *** diff --git a/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effects_supplementary_analyses.Rmd b/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effects_supplementary_analyses.Rmd new file mode 100644 index 00000000..3d443afe --- /dev/null +++ b/interpretation_analysis/indirect_effects/dc_tapseq_indirect_effects_supplementary_analyses.Rmd @@ -0,0 +1,286 @@ +--- +title: "DC-TAP-seq direct vs indirect effects supplementary analyses" +author: "Andreas R. Gschwind" +date: "July 23, 2026" +output: + html_document: + toc: true + toc_float: true + code_folding: show +--- + +```{r setupDocument, include=FALSE} +# set output html chunk options +knitr::opts_chunk$set(warning = FALSE, message = FALSE) + +# suppress dplyrs summarize info +options(dplyr.summarize.inform = FALSE) +``` + +```{r requiredPackages, warning=FALSE, message=FALSE} +library(tidyverse) +library(cowplot) +``` + +## Goal +Analyze direct vs. indirect effects across datasets and plot modeled direct effect rates as function +of distance to TSS. +```{r} +# adjust path to indirect effects analysis output as needed to reproduce results +indirect_effects_dir <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/CRISPR_indirect_effects_FDR_check/results" +``` + +*** + +## Used data +Estimated cis and trans effect rates, as well as predicted direct rates and probability of direct vs +indirect effects are used to plot direct and indirect rates. +```{r} +# load cis and trans effect rates +cis_rates <- read_tsv(file.path(indirect_effects_dir, "cis_positive_hit_rates.tsv"), + show_col_types = FALSE) +trans_rates <- read_tsv(file.path(indirect_effects_dir, "trans_positive_hit_rates.tsv"), + show_col_types = FALSE) + +# load direct effect rates per dataset and averaged across datasets +direct_rates_datasets <- read_tsv(file.path(indirect_effects_dir, "direct_effect_models/direct_rates_per_datasets.tsv"), + show_col_types = FALSE) +direct_rates_average <- read_tsv(file.path(indirect_effects_dir, "direct_effect_models/direct_rates_average_across_datasets.tsv"), + show_col_types = FALSE) + +# load direct effect models +models <- readRDS(file.path(indirect_effects_dir, "direct_effect_models/direct_effects_model.rds")) + +# load dataset ids +dataset_ids <- read_tsv(file.path(indirect_effects_dir, "../config/crispr_dataset_ids.tsv"), + show_col_types = FALSE) + +# add dataset ids to use in plots to cis, trans and direct rates +dataset_ids_rates <- select(dataset_ids, dataset, dataset_plots) +cis_rates <- left_join(cis_rates, dataset_ids_rates, by = "dataset") +trans_rates <- left_join(trans_rates, dataset_ids_rates, by = "dataset") +direct_rates_datasets <- left_join(direct_rates_datasets, dataset_ids_rates, by = "dataset") +``` + +*** + +## Model direct effect rate +The cis positive hit rate is a combination of the direct and indirect hit rate. The expected +indirect hit rate in these CRISPR experiments was estimated by the trans positive hit rate. We can +now calculate the estimated direct effect hit rates by subtracting the indirect effects rate +(trans hit rate) from the cis positive hit rate for each distance bin. + +Some distance bins do have a lower (or 0) cis hit rate in some datasets due to the overall sparsity +of positive hits in these experiments. Thus the resulting negative direct hit rate will can be +considered an artifact. + +To model the direct hit rate as a function of distance to TSS, we average the calculated direct +effects hit rate across all dataset per distance bin and fit a powerlaw model. The following plots +show the powerlaw fit visualized on the direct effect hit rates averaged across datasets. +```{r} +# plot direct effects model fit +plot_powerlaw_model_fit <- function(rate_type = c("all", "negative", "positive"), + models, direct_rates_average, log = FALSE, ylim = NA) { + + rate_type <- match.arg(rate_type) + + # get direct rates for given rate type + rate_type_data <- paste0("direct_rate_", rate_type) # rate type in data and model list + direct_rates_avg_type <- filter(direct_rates_average, type == rate_type_data) + + # get model coefficients for given rate type + model_type <- models[[rate_type_data]] + a <- exp(coef(model_type)[1]) + b <- coef(model_type)[2] + + # plot model fit vs average direct rate per dataset + p <- ggplot(direct_rates_avg_type, aes(x = dist_to_tss, y = direct_rate)) + + geom_point(color = "steelblue") + + stat_function(fun = function(x) a * x^b, color = "black") + + labs(title = str_to_sentence(paste("Direct hit rate", rate_type, "effects")), + x = "Distance to TSS (bp)", y = "Direct hit rate") + + scale_x_continuous(labels = scales::label_scientific(digits = 1)) + + scale_y_continuous(limits = c(NA, ylim)) + + theme_classic() + + theme(text = element_text(size = 13)) + + # transform plots to log space if specified + if (log == TRUE) { + p <- p + + scale_x_log10() + + scale_y_log10(limits = c(NA, ylim)) + + annotation_logticks() + } + + return(p) + +} +``` + +Plot effect rates in log space. Note that this introduces infinite values in the direct hit rates +because some distance bins have no positives due to overall data sparsity. +```{r, fig.height=4, fig.width=8} +# plot all model fits +rate_types <- c(negative = "negative", positive = "positive") +model_fit_plots <- lapply(rate_types, FUN = plot_powerlaw_model_fit, models = models, + direct_rates_average = direct_rates_average, + ylim = 0.2, log = TRUE) + +# arrange all plots into on figure +p1 <- plot_grid(plotlist = model_fit_plots, nrow = 1) +p1 +``` + +*** + +## Predicted direct effects rate +Using the powerlaw model, we can predict the direct effects rate as function of distance to TSS. The +following plots show the predicted direct effect rates as function of distance to TSS for both +negative and positive effects on expression. Empirically calculated indirect effect rates are +overlayed as horizontal lines. +```{r} +# function to predict the direct hit rate for a given distance +predict_direct_hit_rate <- function(models, new_data, type) { + + # get model for specified type + type <- paste0("direct_rate_", type) + fit <- models[[type]] + + # predict direct hit rate for provided distances and exponentiate to get back to original scale + predicted_log_direct_rate <- predict(fit, newdata = new_data) + predicted_direct_rate <- exp(predicted_log_direct_rate) + + # cap predicted values at 1, because the rate can't be higher than that and add to input data + new_data$predicted_direct_rate <- pmin(predicted_direct_rate, 1) + + return(new_data) + +} +``` + +```{r, fig.height=4, fig.width=10} +# convert trans rates to long format and calculate average indirect rate across datasets +trans_rates <- trans_rates %>% + select(dataset_plots, all = positive_rate_significant, negative = positive_rate_negative, + positive = positive_rate_positive) %>% + pivot_longer(cols = -dataset_plots, names_to = "effect_direction", values_to = "indirect_rate") %>% + filter(effect_direction %in% rate_types) + +# predict direct effects rate from 1bp to 1Mb distance to TSS and combine into one table +direct_effects_1mb <- lapply(rate_types, FUN = predict_direct_hit_rate, models = models, + new_data = tibble(dist_to_tss = 1:1e6)) %>% + bind_rows(.id = "effect_direction") +``` + +```{r, fig.height=3, fig.width=8} +# prettier effect directions for plots +direct_effects_1mb <- direct_effects_1mb %>% + mutate(effect_direction = paste(str_to_title(effect_direction), "effects on expression")) +trans_rates <- trans_rates %>% + mutate(effect_direction = paste(str_to_title(effect_direction), "effects on expression")) + +# compute average indirect effect rates per effect size direction +avg_trans_rates <- trans_rates %>% + group_by(effect_direction) %>% + summarize(hit_rate = mean(indirect_rate)) + +# add rate labels for plots +direct_effects_1mb <- mutate(direct_effects_1mb, rate_label = "Predicted direct effects") +avg_trans_rates <- mutate(avg_trans_rates, rate_label = "Avg. indirect effects rate") + +# plot predicted effects rate as function of distance to TSS plus average indirect rate +p2 <- ggplot(direct_effects_1mb, + aes(x = dist_to_tss, y = predicted_direct_rate, color = rate_label)) + + facet_wrap(~effect_direction, scales = "free_y") + + geom_line() + + geom_hline(data = avg_trans_rates, aes(yintercept = hit_rate, color = rate_label)) + + labs(x = "Distance to TSS (bp)", y = "Effects rate", color = "Effects rate:") + + scale_color_manual(values = c(`Avg. indirect effects rate` = "red", + `Predicted direct effects` = "black")) + + scale_x_continuous(labels = scales::label_scientific()) + + scale_y_continuous(limits = c(0, 1)) + + theme_classic() + + theme(strip.background = element_blank(), text = element_text(size = 13), + strip.text = element_text(hjust = 0, size = 16, margin = margin(b = 8))) + +p2 +``` + +```{r, fig.height=3, fig.width=8} +# compute average indirect effect rates only for datasets other than DC-TAP-seq +trans_rates_plot <- trans_rates %>% + mutate(rate_label = case_when( + dataset_plots == "DC-TAP-seq (K562)" ~ "K562 indirect effects rate", + dataset_plots == "DC-TAP-seq (WTC11)" ~ "WTC11 indirect effects rate", + TRUE ~ "Avg. indirect effects rate other datasets" + )) %>% + group_by(rate_label, effect_direction) %>% + summarize(hit_rate = mean(indirect_rate)) + +# colors for the different rates +rate_colors <- c( + `Predicted direct effects` = "black", + `K562 indirect effects rate` = "red", + `WTC11 indirect effects rate` = "blue", + `Avg. indirect effects rate other datasets` = "#637a91" + +) + +# plot rates and clip y-axis at 0.1 for better visualization +p3 <- ggplot(direct_effects_1mb, + aes(x = dist_to_tss, y = predicted_direct_rate, color = rate_label)) + + facet_wrap(~effect_direction, scales = "free_y") + + geom_line() + + geom_hline(data = trans_rates_plot, aes(yintercept = hit_rate, color = rate_label)) + + labs(x = "Distance to TSS (bp)", y = "Effects rate\n(Displayed range: [0, 0.1])", + color = "Effects rate:") + + scale_color_manual(values = rate_colors) + + scale_x_continuous(labels = scales::label_scientific()) + + coord_cartesian(ylim = c(0, 0.1)) + + theme_classic() + + theme(strip.background = element_blank(), text = element_text(size = 13), + strip.text = element_text(hjust = 0, size = 16, margin = margin(b = 8))) + +p3 +``` + +*** + +## Arrange supplementary figure +```{r, fig.height=8, fig.width=8} +# arrange panels into figure +plot_grid( + p1, + p2 + theme(legend.position = "bottom"), + labels = c("a", "b"), + ncol = 1, + scale = 0.9, + rel_heights = c(1, 1.1) + ) + +# save figure to pdf +ggsave(file = "plots/FigS15_direct_indirect_rates_avg.pdf", height = 8, width = 8, create.dir = TRUE) +``` + + +```{r, fig.height=8, fig.width=8} +# arrange panels into figure +plot_grid( + p1, + p3 + theme(legend.position = "bottom") + guides(color = guide_legend(nrow = 2)), + labels = c("a", "b"), + ncol = 1, + scale = 0.9, + rel_heights = c(1, 1.2) + ) + +# save figure to pdf +ggsave(file = "plots/FigS15_direct_indirect_rates_k562_wtc11.pdf", height = 8, width = 8, create.dir = TRUE) +``` +*** + +## Session information +This document was produced using following packages: +```{r sessionInfo} +sessionInfo() +``` diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_create_gene_universe.R b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_create_gene_universe.R new file mode 100644 index 00000000..61c2b5d3 --- /dev/null +++ b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_create_gene_universe.R @@ -0,0 +1,54 @@ +## create gene universe for all genes in Gasperini et al., 2019 + +suppressPackageStartupMessages({ + library(tidyverse) + library(rtracklayer) + library(sceptre) +}) + +# load table with main screen results +screen_results_file <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/analyses/dc_tapseq_paper/Table_S3_Final_DC_TAP_Seq_Results_w_Chromatin_Categories_on_resized_and_merged_elements_250908_append.tsv" +screen_results <- read_tsv(screen_results_file, show_col_types = FALSE) + +# filter for K562 results only +k562_screen_results <- filter(screen_results, cell_type == "K562") + +# load genome annotations +annot <- import("../../../results/genome_annotation_files/gencode.v26lift37.annotation.gtf.gz", + format = "gtf") + +# extract gene symbols, ids (without version) and chromosomes for all genes +genes <- annot[annot$type == "gene"] %>% + as.data.frame() %>% + select(gene_id, gene_name, gene_type, gene_chr = seqnames) %>% + mutate(gene_id = sub("\\..+", "", gene_id)) %>% + distinct() + +# only retain annotations on regular chromosomes and from protein-coding or lincRNA genes +regular_chrs <- paste0("chr", c(1:22, "X", "M")) +genes_filt <- genes %>% + filter(gene_chr %in% regular_chrs) %>% + filter(gene_type %in% c("protein_coding", "lincRNA")) + +# load Gasperini et al., 2019 sceptre object containing UMI counts used to compute TPM values +gasp_sceptre_file <- "/oak/stanford/groups/engreitz/Users/jgalante/DC_TAP_Paper/results/main_figure_1_and_2/duplicate_pairs_analysis/differential_expression/sceptre_diffex_input.rds" +gasp_sceptre <- readRDS(gasp_sceptre_file) + +# calculate TPM for each gene +umi_per_gene <- rowSums(gasp_sceptre@response_matrix[[1]]) +tpm_per_gene <- umi_per_gene * 1e6 / sum(umi_per_gene) +tpm_per_gene <- enframe(tpm_per_gene, name = "gene_id", value = "tpm") + +# add TPM to genes table +genes_filt <- left_join(genes_filt, tpm_per_gene, by = "gene_id") + +# label genes that were part of the K562 DC-TAP-seq enhancer screen +k562_dc_tap_genes <- unique(k562_screen_results$gene_id) +dc_tap_genes_in_gasp_genes <- sum(k562_dc_tap_genes %in% genes_filt$gene_id) +message(dc_tap_genes_in_gasp_genes, " out of ", length(k562_dc_tap_genes), + " K562 DC-TAP-seq genes found in Gasperini et al., 2019 data") +genes_filt <- mutate(genes_filt, k562_dc_tap = gene_id %in% k562_dc_tap_genes) + +# save filtered and annotated gene list to output file +write_tsv(genes_filt, file = "gasperini_genes_annotated.tsv.gz") + \ No newline at end of file diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_find_dctap_perts_in_gasperini.R b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_find_dctap_perts_in_gasperini.R new file mode 100644 index 00000000..e34f249f --- /dev/null +++ b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/00_find_dctap_perts_in_gasperini.R @@ -0,0 +1,51 @@ +## Intersect coordinates of K562 DC-TAP-seq candidate elements with Gasperini et al., 2019 +## candidate elements to find elements perturbed in both screens + +suppressPackageStartupMessages({ + library(data.table) + library(dplyr) + library(tidyr) + library(sceptre) + library(GenomicRanges) +}) + +# load DC-TAP-seq cis-analysis results +screen_results_file <- "/oak/stanford/groups/engreitz/Users/agschwin/distal_regulation_paper/analyses/dc_tapseq_paper/Table_S3_Final_DC_TAP_Seq_Results_w_Chromatin_Categories_on_resized_and_merged_elements_250908_append.tsv" +screen_results <- fread(screen_results_file) + +# filter for K562 results only +k562_screen_results <- filter(screen_results, cell_type == "K562") + +# load Gasperini et al., 2019 sceptre object and extract grna targets table with candidate elements +gasp <- read_ondisc_backed_sceptre_object( + sceptre_object_fp = "sceptre_object_ondisc/sceptre_object.rds", + response_odm_file_fp = "sceptre_object_ondisc/response.odm", + grna_odm_file_fp = "sceptre_object_ondisc/grna.odm" +) +gasp_grna_targets <- gasp@grna_target_data_frame + +# create GenomicRanges for DC-TAP-seq candidate elements +dctap_elements <- k562_screen_results %>% + select(intended_target_name_hg38, chr = targeting_chr_hg19, start = targeting_start_hg19, + end = targeting_end_hg19) %>% + distinct() %>% + makeGRangesFromDataFrame(keep.extra.columns = TRUE, starts.in.df.are.0based = TRUE) + +# create GenomicRanges for Gasperini candidate elements +gasp_elements <- gasp_grna_targets %>% + select(grna_target) %>% + distinct() %>% + separate(grna_target, into = c("chr", "start", "end"), sep = ":|-", remove = FALSE) %>% + makeGRangesFromDataFrame(keep.extra.columns = TRUE, starts.in.df.are.0based = TRUE) + +# intersect K562 DC-TAP-seq and Gasperini elements +ovl <- as.data.frame(findOverlapPairs(dctap_elements, gasp_elements)) + +# reformat table with DC TAP-seq elements overlapping Gasperini elements +output <- ovl %>% + select(dc_tapseq_element = `first.intended_target_name_hg38`, + gasperini_element = `second.grna_target`, gasp_chr = second.X.seqnames, + gasp_start = second.X.start, gasp_end = second.X.end) + +# write to output file +fwrite(output, file = "overlapping_elements.tsv", sep = "\t") diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/1_create_ondisc_sceptre_object.R b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/1_create_ondisc_sceptre_object.R new file mode 100644 index 00000000..a08868e0 --- /dev/null +++ b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/1_create_ondisc_sceptre_object.R @@ -0,0 +1,65 @@ +## Create sceptre on-disc object for Gasperini et al., 2019 data for use with Nextflow pipeline + +suppressPackageStartupMessages({ + library(sceptre) + library(ondisc) +}) + +# load in-memory sceptre object from main analyses of Gasperini et al dataset +sceptre_object_file <- "/oak/stanford/groups/engreitz/Users/jgalante/DC_TAP_Paper/results/main_figure_1_and_2/duplicate_pairs_analysis/differential_expression/sceptre_diffex_input.rds" +sceptre_object <- readRDS(sceptre_object_file) + +# extract used formula and save to formula object file for sceptre Nextflow pipeline +formula <- sceptre_object@formula_object +saveRDS(formula, file = "sceptre_formula.rds") + +# get used MOI setting +moi <- ifelse(isTRUE(sceptre_object@low_moi), yes = "low", no = "high") + +# extract gene (response) and grna count matrices +response_matrix <- sceptre_object@response_matrix[[1]] +grna_matrix <- sceptre_object@grna_matrix[[1]] + +# extract required columns from gRNA targets data frame +grna_target_df <- sceptre_object@grna_target_data_frame[, c("grna_id", "grna_target")] + +# extract covariates data frame +covariate_df <- sceptre_object@covariate_data_frame + +# remove any covariates computed internally when creating the new sceptre object +computed_covars <- c("response_n_nonzero", "response_n_umis", "grna_n_nonzero", "grna_n_umis", + "response_p_mito") +extra_cols <- setdiff(colnames(covariate_df), computed_covars) +if (length(extra_cols) > 0L) { + extra_covars <- covariate_df[, extra_cols, drop = FALSE] +} else { + extra_covars <- data.frame() +} + +# create new on-disc sceptre object +output_dir <- "sceptre_object_ondisc" +dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) +ondisc_sceptre_object <- import_data( + response_matrix = response_matrix, + grna_matrix = grna_matrix, + grna_target_data_frame = grna_target_df, + moi = moi, + extra_covariates = extra_covars, + use_ondisc = TRUE, + directory_to_write = output_dir +) + +# get side parameter for differential expression testsfrom input sceptre object +side <- switch(as.character(sceptre_object@side_code), "-1" = "left", "0" = "both", "1" = "right", + stop("Unexpected side_code: ", sceptre_object@side_code)) + +# set analysis parameters +ondisc_sceptre_object <- set_analysis_parameters( + sceptre_object = ondisc_sceptre_object, + side = side, + formula_object = formula, + resampling_mechanism = "permutations" +) + +# save ondisc sceptre object to file +write_ondisc_backed_sceptre_object(ondisc_sceptre_object, directory_to_write = output_dir) diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/2_run_sceptre_nextflow.sh b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/2_run_sceptre_nextflow.sh new file mode 100644 index 00000000..4bd9a8dc --- /dev/null +++ b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/2_run_sceptre_nextflow.sh @@ -0,0 +1,46 @@ +#!/bin/bash + +## Run SCEPTRE Nextflow pipeline to perform all-against-all trans analysis for the Gasperini et al., +## 2019 dataset. Requires that Nextflow and R are available in the PATH, and the sceptre and ondisc +## R-packages are installed. + +# load required modules +ml R/4.2.0 +ml java/21.0.4 +ml biology +ml nextflow/25.04.7 + +########################## +# REQUIRED INPUT ARGUMENTS +########################## +data_directory="./sceptre_object_ondisc" +# sceptre object +sceptre_object_fp=$data_directory"/sceptre_object.rds" +# response ODM +response_odm_fp=$data_directory"/response.odm" +# grna ODM +grna_odm_fp=$data_directory"/grna.odm" +# object containing model formula +formula_object="./sceptre_formula.rds" + +################### +# OUTPUT DIRECTORY: +################## +output_directory="./sceptre_outputs" + +################# +# Invoke pipeline +################# +nextflow run timothy-barry/sceptre-pipeline -r main \ +--sceptre_object_fp $sceptre_object_fp \ +--response_odm_fp $response_odm_fp \ +--grna_odm_fp $grna_odm_fp \ +--output_directory $output_directory \ +--formula_object $formula_object \ +--grna_assignment_method "thresholding" \ +--threshold 1 \ +--pair_pod_size 250000 \ +--run_association_analysis_time_per_pair 0.25s \ +--run_association_analysis_memory 32G \ +--discovery_pairs trans \ +-resume diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/3_calculate_indirect_effects.R b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/3_calculate_indirect_effects.R new file mode 100644 index 00000000..ce5c98bf --- /dev/null +++ b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/3_calculate_indirect_effects.R @@ -0,0 +1,103 @@ +## Process massive trans analysis output from sceptre, and compute and analyze indirect effect rates + +suppressPackageStartupMessages({ + library(tidyverse) + library(arrow) +}) + +# load Gasperini et al., 2019 gene universe, including chromosome, TPM and DC_TAP-seq screen annotations +genes <- read_tsv("gasperini_genes_annotated.tsv.gz", show_col_types = FALSE) + +# load elements overlapping between K562 DC-TAP-seq and Gasperini et al. +overlapping_elements <- read_tsv("overlapping_elements.tsv", show_col_types = FALSE) + +# load cis-analysis results +cis_results <- read_tsv("/oak/stanford/groups/engreitz/Users/jgalante/DC_TAP_Paper/results/main_figure_1_and_2/duplicate_pairs_analysis/results_with_element_gene_pair_categories.tsv", + show_col_types = FALSE) + +# filter for distal element-gene pairs only (remove promoters) +cis_results <- filter(cis_results, DistalElement_Gene == TRUE) + +# extract p-value cutoff for distal element-gene interactions +pval_threshold <- cis_results %>% + filter(DistalElement_Gene == TRUE, abs(distance_to_abc_canonical_TSS) <= 1e6) %>% + filter(significant == TRUE) %>% + slice_max(sceptre_p_value) %>% + pull(sceptre_p_value) + +# Schema defining column formats for sceptre trans analysis results. Each parquet file carries its +# own dictionary, and Arrow cannot unify differing dictionaries across files during a +# compute/join. If column formats aren't set before any operation that materializes columns +# (e.g. the distinct() below), collect() or compute() fails with "Unifying differing dictionaries" +sceptre_schema <- schema( + response_id = utf8(), + grna_target = utf8(), + n_nonzero_trt = int32(), + n_nonzero_cntrl = int32(), + pass_qc = bool(), + p_value = float64(), + fold_change = float64(), + se_fold_change = float64(), + log_2_fold_change = float64() +) + +# open sceptre results +trans_results <- open_dataset("sceptre_outputs/trans_results", schema = sceptre_schema) + +# cast the dictionary-encoded id columns to plain strings. Each parquet file carries its +# own dictionary, and Arrow cannot unify differing dictionaries across files during a +# compute/join. This must happen before any operation that materializes these columns +# (e.g. the distinct() below), otherwise collect() fails with "Unifying differing dictionaries" +trans_results <- trans_results %>% + mutate(grna_target = cast(grna_target, string()), + response_id = cast(response_id, string())) + +# only retain unique tested pairs that did pass sceptre QC +trans_results <- trans_results %>% + filter(pass_qc == TRUE) %>% + distinct() + +# extract chromosome for each perturbation from grna_target id (format: "{chr}:{start}-{end}") +trans_results <- mutate(trans_results, pert_chr = sub("(chr.+):.+", "\\1", grna_target)) + +# add gene annotations to sceptre output +trans_results <- left_join(trans_results, genes, by = c("response_id" = "gene_id")) + +# filter for trans-acting interactions and only retain pairs involving genes with TPM values +# (or DC-TAP-seq genes) or distal elements +distal_elements <- unique(pull(cis_results, intended_target_name_hg19)) +trans_results <- trans_results %>% + filter(pert_chr != gene_chr) %>% + filter(!is.na(tpm) | k562_dc_tap == TRUE) %>% + filter(grna_target %in% distal_elements) + +# label elements that overlap K562 DC-TAP-seq candidate elements +trans_results <- trans_results %>% + mutate(k562_dc_tap_gene = k562_dc_tap, + k562_dc_tap_element = grna_target %in% overlapping_elements$gasperini_element) + +# re-calculate significance based on nominal p-value cutoff from cis-analysis +trans_results <- mutate(trans_results, significant = p_value <= pval_threshold) + +# perform all computations and materialize the filtered results into memory +trans_results <- collect(trans_results) + +# compute indirect effects rate per gene +indirect_effects_genes <- trans_results %>% + group_by(response_id, gene_name, tpm, k562_dc_tap_gene) %>% + summarize(tested_trans_pairs = n(), + significant_trans_pairs = sum(significant), + trans_hit_rate = significant_trans_pairs / tested_trans_pairs, + .groups = "drop") + +# compute indirect effects rate per element +indirect_effects_elements <- trans_results %>% + group_by(grna_target, k562_dc_tap_element) %>% + summarize(tested_trans_pairs = n(), + significant_trans_pairs = sum(significant), + trans_hit_rate = significant_trans_pairs / tested_trans_pairs, + .groups = "drop") + +# save indirect effects tables to output file +write_tsv(indirect_effects_genes, file = "gasperini2019_indirect_effects_per_gene.tsv.gz") +write_tsv(indirect_effects_elements, file = "gasperini2019_indirect_effects_per_element.tsv.gz") diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/4_analyze_indirect_effects.Rmd b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/4_analyze_indirect_effects.Rmd new file mode 100644 index 00000000..a0634b44 --- /dev/null +++ b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/4_analyze_indirect_effects.Rmd @@ -0,0 +1,200 @@ +--- +title: "Genome-wide indirect effects" +author: "Andreas R. Gschwind" +date: "July 23, 2026" +output: + html_document: + toc: true + toc_float: true + code_folding: show +--- + +```{r setupDocument, include=FALSE} +# set output html chunk options +knitr::opts_chunk$set(warning = FALSE) + +# suppress dplyrs summarize info +options(dplyr.summarize.inform = FALSE) +``` + +```{r attachPackages, message=FALSE, warning=FALSE} +library(tidyverse) +library(cowplot) +``` + +## Goal +Analyze genome-wide indirect effects calculated for all perturbations and genes in the Gasperini et +al., 2019 dataset. + +*** + +## Used data +Load indirect effect rates for all genes and elements in the Gasperini et al., 2019 datasets +calculated from sceptre differential expression tests for trans-acting effects of all enhancer perturbations. +```{r} +# load indirect effect rates per gene +indirect_effects_genes <- read_tsv("gasperini2019_indirect_effects_per_gene.tsv.gz", show_col_types = FALSE) +indirect_effects_elements <- read_tsv("gasperini2019_indirect_effects_per_element.tsv.gz", show_col_types = FALSE) +``` + +*** + +## Indirect effects vs. TPM +Plot relationship between indirect effect rate vs. target gene TPM. +```{r, fig.height=4, fig.width=6} +# create plotting labels for genes in K562 DC-TAP-seq target gene panels +indirect_effects_genes <- indirect_effects_genes %>% + mutate(dc_tap_group = if_else(k562_dc_tap_gene == TRUE, true = "K562 DC-TAP-seq targets", + false = "Other genes")) %>% + mutate(dc_tap_group = fct_relevel(dc_tap_group, "Other genes")) %>% + arrange(dc_tap_group) + +# filter for genes within custom TPM range for plot +indirect_effects_tpm_range <- indirect_effects_genes %>% + filter(tpm >= 10, tpm <= 1e4) + +# plot indirect effects rate vs target gene TPM +ggplot(indirect_effects_tpm_range, aes(x = tpm, y = trans_hit_rate, color = dc_tap_group)) + + geom_point() + + geom_smooth(aes(group = 1), method = "lm", se = FALSE, color = "black") + + labs(title = "Indirect effects rate vs. TPM", x = "Target gene TPM", + y = "Indirect effects rate", color = "Gene type") + + scale_color_manual(values = c("K562 DC-TAP-seq targets" = "steelblue", "Other genes" = "gray")) + + scale_x_log10() + + theme_classic() +``` + +The positive relationship between the observed indirect effects rate and target gene expression +level is likely driven by higher statistical power to detect differential expression effects for +genes with higher UMI counts. + +Therefore to compare indirect effect rates between DC-TAP-seq target genes and other genes, we need +to match TPM distributions of the two groups. + +## Compare indirect effects +Compare indirect effects rates between K562 DC-TAP-seq target genes versus a sample of other genes +expressed with similar TPM distribution. +```{r} +# split genes into the two groups to compare +target_genes <- filter(indirect_effects_genes, k562_dc_tap_gene == TRUE) +other_genes <- filter(indirect_effects_genes, k562_dc_tap_gene == FALSE) + +# restrict other genes to the TPM range spanned by the DC TAP-seq target genes, so the matched +# sample can only contain genes expressed within that range +tpm_range <- range(target_genes$tpm) +other_genes <- filter(other_genes, tpm >= tpm_range[1], tpm <= tpm_range[2]) + +# define TPM bins (log10 scale) from deciles of the target genes' TPM distribution; the outer +# breaks equal the min/max target TPM, which now also bound the other genes +tpm_breaks <- quantile(log10(target_genes$tpm), probs = seq(0, 1, length.out = 11)) + +target_genes <- mutate(target_genes, tpm_bin = cut(log10(tpm), tpm_breaks, include.lowest = TRUE)) +other_genes <- mutate(other_genes, tpm_bin = cut(log10(tpm), tpm_breaks, include.lowest = TRUE)) + +# number and proportion of target genes per TPM bin (the TPM distribution to match), plus the number +# of other genes available in each bin +n_per_bin <- target_genes %>% + count(tpm_bin, name = "n_target") %>% + mutate(prop_target = n_target / sum(n_target)) %>% + left_join(count(other_genes, tpm_bin, name = "n_other_avail"), by = "tpm_bin") %>% + mutate(n_other_avail = coalesce(n_other_avail, 0L)) + +# sample as many "Other genes" as possible while keeping their per-bin proportions identical to the +# target genes' TPM distribution. The largest achievable total is set by the bin that runs out of +# available genes first: scale each bin's target proportion up by the maximum factor for which no bin +# exceeds its available genes, then take that many per bin. +scale_factor <- min(n_per_bin$n_other_avail / n_per_bin$prop_target) +n_per_bin <- mutate(n_per_bin, n_other = floor(scale_factor * prop_target)) + +# sample the computed number of "Other genes" per TPM bin +set.seed(20260710) +matched_other_genes <- other_genes %>% + group_by(tpm_bin) %>% + group_modify(~ slice_sample(.x, n = n_per_bin$n_other[n_per_bin$tpm_bin == .y$tpm_bin])) %>% + ungroup() + +# combine target genes with the TPM-matched sample of other genes +indirect_effects_tpm_matched <- bind_rows(target_genes, matched_other_genes) +``` + +Sanity check: the TPM distributions of the two groups should overlap after matching. +```{r, fig.height=4, fig.width=6} +# plot TPM density of both groups to confirm the matching worked +ggplot(indirect_effects_tpm_matched, aes(x = tpm, fill = dc_tap_group, color = dc_tap_group)) + + geom_density(alpha = 0.3) + + labs(title = "TPM distributions after matching", x = "Target gene TPM", y = "Density", + fill = "Gene type", color = "Gene type") + + scale_color_manual(values = c("K562 DC-TAP-seq targets" = "steelblue", "Other genes" = "gray")) + + scale_fill_manual(values = c("K562 DC-TAP-seq targets" = "steelblue", "Other genes" = "gray")) + + scale_x_log10() + + theme_classic() +``` + +```{r, fig.height=4, fig.width=4} +# test for a difference in indirect effect rates between the two groups. The groups are independent +# and the rates are non-normal, so use a two-sided Wilcoxon rank-sum (Mann-Whitney U) test +rate_test <- wilcox.test(trans_hit_rate ~ dc_tap_group, data = indirect_effects_tpm_matched) +rate_test + +# add sample sizes to labels +indirect_effects_tpm_matched <- indirect_effects_tpm_matched %>% + add_count(k562_dc_tap_gene, name = "n") %>% + mutate(dc_tap_group = paste0(dc_tap_group, "\n(n = ", n, " )" )) + +# compare indirect effect rates between K562 target genes any other 30+ TPM genes +p1 <- ggplot(indirect_effects_tpm_matched, + aes(x = dc_tap_group, y = trans_hit_rate, color = dc_tap_group)) + + geom_jitter(width = 0.2) + + geom_boxplot(outlier.shape = NA, fill = NA, color = "black", width = 0.25, notch = TRUE) + + labs(y = "Indirect effects rate", title = "Indirect effects per gene", + subtitle = paste0("Wilcoxon rank-sum p = ", signif(rate_test$p.value, 2))) + + scale_color_manual(values = c("steelblue", "gray")) + + theme_classic() + + theme(legend.position = "none", axis.title.x = element_blank()) + +p1 +``` + +*** + +## Per-element indirect effect rates +```{r, fig.height=4, fig.width=4} +# test for a difference in indirect effect rates between the two groups. The groups are independent +# and the rates are non-normal, so use a two-sided Wilcoxon rank-sum (Mann-Whitney U) test +rate_test <- wilcox.test(trans_hit_rate ~ k562_dc_tap_element, data = indirect_effects_elements) +rate_test + +# add new labels for elements +indirect_effects_elements <- indirect_effects_elements %>% + mutate(element_group = if_else(k562_dc_tap_element == TRUE, true = "K526 DC-TAP-seq element", + false = "Other elements")) %>% + add_count(k562_dc_tap_element, name = "n") %>% + mutate(element_group = paste0(element_group, "\n(n = ", n, " )" )) + +# make plot +p2 <- ggplot(indirect_effects_elements, + aes(x = element_group, y = trans_hit_rate, color = element_group)) + + geom_jitter(width = 0.2) + + geom_boxplot(outlier.shape = NA, fill = NA, color = "black", width = 0.25, notch = TRUE) + + labs(y = "Indirect effects rate", title = "Indirect effects per element", + subtitle = paste0("Wilcoxon rank-sum p = ", signif(rate_test$p.value, 2))) + + scale_y_sqrt() + + scale_color_manual(values = c("steelblue", "gray")) + + theme_classic() + + theme(legend.position = "none", axis.title.x = element_blank()) + +p2 +``` + +```{r, fig.height=3, fig.width=7} +# arrange plots into figure +plot_grid(p1, p2, nrow = 1) +ggsave(filename = "../plots/indirect_rates_gasperini.pdf", height = 3, width = 7, create.dir = TRUE) +``` + + +## Session information +This document was produced using following packages: +```{r sessionInfo} +sessionInfo() +``` diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/nextflow.config b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/nextflow.config new file mode 100644 index 00000000..c73f442a --- /dev/null +++ b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/nextflow.config @@ -0,0 +1,5 @@ +process.executor = 'slurm' +executor.queueSize = 300 +executor.submitRateLimit = '10 sec' +workDir = 'work' +process.queue = 'normal' diff --git a/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/sceptre_formula.rds b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/sceptre_formula.rds new file mode 100644 index 0000000000000000000000000000000000000000..6d51bb89007ef914f6a7448e2261a63b6f44f90b GIT binary patch literal 176 zcmV;h08jrPiwFP!000001B>8dU|?WoU||B1tUx9MYiNj@t_6_C#sC&&U}Run-~>q| z=Oh*v0~rEHg6wJeMY*LpiD+th^j!1GGK=!_a#Qn4fc(Ee{0~W