Update direct vs. indirect effects analysis code

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

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
+

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/<dataset>/`:
+  `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.

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()
+```