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 00000000..6d51bb89 Binary files /dev/null and b/interpretation_analysis/indirect_effects/genome_wide_indirect_effects/sceptre_formula.rds differ