|
#!/bin/bash |
|
#=============================================================================== |
|
# |
|
# VEDS ATAC-seq Complete Analysis Pipeline |
|
# |
|
# Associated Publication: |
|
# Juzwiak EE, Bowen CJ, Edwards R, Restrepo L, et al. (2025) |
|
# Steroid hormone antagonism affords vascular protection in a mouse |
|
# model of vascular Ehlers-Danlos syndrome. |
|
# JCI Insight. https://doi.org/10.1172/jci.insight.198202 |
|
# |
|
# Description: |
|
# Complete ATAC-seq analysis pipeline including: |
|
# - Section 1: Peak calling with MACS2 |
|
# - Section 2: Quality control metrics (FRIP, TSS enrichment) |
|
# - Section 3: Peak annotation with ChIPseeker |
|
# - Section 4: Sequence modeling with gkm-SVM |
|
# - Section 5: Motif analysis with gkmPWM |
|
# - Section 6: Transcription factor binding site prediction |
|
# |
|
# Data Availability: |
|
# Raw FASTQ files: GEO GSE315744 |
|
# |
|
# Code References: |
|
# - gkmPWM: https://github.com/shigakiD/gkmPWM |
|
# - gkmSVM: https://github.com/Dongwon-Lee/lsgkm |
|
# |
|
# Authors: |
|
# Analysis performed by: Leda Restrepo (LR) |
|
# Code generated by: Dustin T. Shigaki (DTS), Michael A. Beer (MAB) |
|
# |
|
# Software Requirements: |
|
# - MACS2 (v2.2.7.1) |
|
# - Samtools |
|
# - Bedtools |
|
# - deepTools |
|
# - R (v4.0+) with packages: |
|
# - ChIPseeker (v1.40.0) |
|
# - TxDb.Mmusculus.UCSC.mm10.knownGene |
|
# - org.Mm.eg.db |
|
# - GenomicRanges |
|
# - ggplot2 |
|
# - lsgkm (gkmSVM) |
|
# - gkmPWM (Python package) |
|
# - Python 3 with numpy, scipy |
|
# |
|
#=============================================================================== |
|
|
|
#=============================================================================== |
|
# USER-DEFINED PARAMETERS |
|
#=============================================================================== |
|
|
|
# Input directories |
|
BAM_DIR="path/to/aligned/bams" |
|
BLACKLIST="path/to/ENCODE_mm10_blacklist.bed" |
|
TSS_BED="path/to/mm10_tss.bed" |
|
GENOME_FA="path/to/mm10.fa" |
|
|
|
# gkm-SVM parameters |
|
GKM_L=11 # full gapped k-mer length |
|
GKM_K=7 # number of non-gapped positions |
|
GKM_T=4 # number of threads |
|
|
|
# Number of threads |
|
THREADS=4 |
|
|
|
#=============================================================================== |
|
# SECTION 1: PEAK CALLING WITH MACS2 |
|
#=============================================================================== |
|
|
|
mkdir -p ${PEAK_DIR} |
|
|
|
# Call peaks with MACS2 |
|
|
|
for SAMPLE in "${SAMPLES[@]}"; do |
|
|
|
macs2 callpeak \ |
|
-t ${BAM_DIR}/${SAMPLE}.bam \ |
|
-f BAMPE \ |
|
-g mm \ |
|
-n ${SAMPLE} \ |
|
--outdir ${PEAK_DIR} \ |
|
-q 0.05 \ |
|
--nomodel \ |
|
--shift -100 \ |
|
--extsize 200 \ |
|
--keep-dup all \ |
|
-B --SPMR \ |
|
2> ${PEAK_DIR}/${SAMPLE}_macs2.log |
|
done |
|
|
|
|
|
for SAMPLE in "${SAMPLES[@]}"; do |
|
bedtools intersect \ |
|
-a ${PEAK_DIR}/${SAMPLE}_peaks.narrowPeak \ |
|
-b ${BLACKLIST} \ |
|
-v \ |
|
> ${PEAK_DIR}/${SAMPLE}_peaks.filtered.narrowPeak |
|
done |
|
|
|
|
|
cat ${PEAK_DIR}/Control_*_peaks.filtered.narrowPeak | \ |
|
sort -k1,1 -k2,2n | \ |
|
bedtools merge -i - \ |
|
> ${PEAK_DIR}/control_consensus_peaks.bed |
|
|
|
|
|
cat ${PEAK_DIR}/VEDS_*_peaks.filtered.narrowPeak | \ |
|
sort -k1,1 -k2,2n | \ |
|
bedtools merge -i - \ |
|
> ${PEAK_DIR}/veds_consensus_peaks.bed |
|
|
|
|
|
cat ${PEAK_DIR}/*_peaks.filtered.narrowPeak | \ |
|
sort -k1,1 -k2,2n | \ |
|
bedtools merge -i - \ |
|
> ${PEAK_DIR}/all_consensus_peaks.bed |
|
|
|
#=============================================================================== |
|
# SECTION 2: QUALITY CONTROL METRICS |
|
#=============================================================================== |
|
|
|
mkdir -p ${QC_DIR} |
|
mkdir -p ${BIGWIG_DIR} |
|
|
|
#------------------------------------------------------------------------------- |
|
# 2A: Calculate FRIP (Fraction of Reads in Peaks) |
|
#------------------------------------------------------------------------------- |
|
|
|
for SAMPLE in "${SAMPLES[@]}"; do |
|
# Total reads in BAM (excluding unmapped) |
|
TOTAL_READS=$(samtools view -c -F 4 ${BAM_DIR}/${SAMPLE}.bam) |
|
|
|
# Reads overlapping consensus peaks |
|
READS_IN_PEAKS=$(bedtools intersect \ |
|
-a ${BAM_DIR}/${SAMPLE}.bam \ |
|
-b ${PEAK_DIR}/all_consensus_peaks.bed \ |
|
-u -ubam | samtools view -c) |
|
|
|
# Calculate FRIP |
|
FRIP=$(echo "scale=4; ${READS_IN_PEAKS} / ${TOTAL_READS}" | bc) |
|
|
|
echo -e "${SAMPLE}\t${TOTAL_READS}\t${READS_IN_PEAKS}\t${FRIP}" >> ${QC_DIR}/frip_scores.txt |
|
echo " ${SAMPLE}: FRIP = ${FRIP}" |
|
done |
|
|
|
cat ${QC_DIR}/frip_scores.txt |
|
|
|
#------------------------------------------------------------------------------- |
|
# 2B: Generate CPM-normalized bigWigs |
|
#------------------------------------------------------------------------------- |
|
|
|
for SAMPLE in "${SAMPLES[@]}"; do |
|
|
|
bamCoverage \ |
|
-b ${BAM_DIR}/${SAMPLE}.bam \ |
|
-o ${BIGWIG_DIR}/${SAMPLE}.cpm.bw \ |
|
--normalizeUsing CPM \ |
|
--binSize 10 \ |
|
--extendReads \ |
|
-p ${THREADS} |
|
done |
|
|
|
#------------------------------------------------------------------------------- |
|
# 2C: Calculate TSS Enrichment |
|
#------------------------------------------------------------------------------- |
|
|
|
computeMatrix reference-point \ |
|
-S ${BIGWIG_DIR}/*.cpm.bw \ |
|
-R ${TSS_BED} \ |
|
--referencePoint TSS \ |
|
-a 2000 -b 2000 \ |
|
--binSize 10 \ |
|
-o ${QC_DIR}/tss_matrix.gz \ |
|
-p ${THREADS} |
|
|
|
plotProfile \ |
|
-m ${QC_DIR}/tss_matrix.gz \ |
|
-o ${QC_DIR}/tss_enrichment_profile.pdf \ |
|
--plotTitle "TSS Enrichment - VEDS ATAC-seq" \ |
|
--perGroup |
|
|
|
plotHeatmap \ |
|
-m ${QC_DIR}/tss_matrix.gz \ |
|
-o ${QC_DIR}/tss_enrichment_heatmap.pdf \ |
|
--colorMap RdYlBu_r \ |
|
--whatToShow 'heatmap and colorbar' |
|
|
|
#=============================================================================== |
|
# SECTION 3: PEAK ANNOTATION WITH CHIPSEEKER (R SCRIPT) |
|
#=============================================================================== |
|
|
|
mkdir -p ${ANNOT_DIR} |
|
|
|
export PEAK_DIR=${PEAK_DIR} |
|
export ANNOT_DIR=${ANNOT_DIR} |
|
|
|
Rscript - <<'RSCRIPT' |
|
|
|
library(ChIPseeker) |
|
library(TxDb.Mmusculus.UCSC.mm10.knownGene) |
|
library(org.Mm.eg.db) |
|
library(GenomicRanges) |
|
library(ggplot2) |
|
|
|
peak_dir <- Sys.getenv("PEAK_DIR", "path/to/peaks") |
|
output_dir <- Sys.getenv("ANNOT_DIR", "path/to/annotation") |
|
|
|
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE) |
|
|
|
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene |
|
|
|
# Promoter region ±3kb from TSS |
|
promoter_region <- c(-3000, 3000) |
|
|
|
|
|
cat("Loading consensus peaks...\n") |
|
|
|
consensus_peaks <- read.table( |
|
file.path(peak_dir, "all_consensus_peaks.bed"), |
|
header = FALSE, |
|
col.names = c("chr", "start", "end") |
|
) |
|
|
|
peaks_gr <- GRanges( |
|
seqnames = consensus_peaks$chr, |
|
ranges = IRanges(start = consensus_peaks$start, end = consensus_peaks$end) |
|
) |
|
|
|
peaks_gr$peak_id <- paste0( |
|
seqnames(peaks_gr), ":", |
|
start(peaks_gr), "-", |
|
end(peaks_gr) |
|
) |
|
|
|
cat("Total peaks to annotate:", length(peaks_gr), "\n") |
|
|
|
cat("\nAnnotating peaks with ChIPseeker...\n") |
|
|
|
peak_anno <- annotatePeak( |
|
peaks_gr, |
|
tssRegion = promoter_region, |
|
TxDb = txdb, |
|
annoDb = "org.Mm.eg.db", |
|
verbose = TRUE |
|
) |
|
|
|
anno_df <- as.data.frame(peak_anno) |
|
|
|
# Classify as promoter-proximal vs distal |
|
anno_df$peak_type <- ifelse( |
|
abs(anno_df$distanceToTSS) <= 3000, |
|
"Promoter-proximal", |
|
"Distal" |
|
) |
|
|
|
cat("\n=== Annotation Summary ===\n") |
|
print(peak_anno) |
|
|
|
cat("\n=== Peak Classification ===\n") |
|
print(table(anno_df$peak_type)) |
|
|
|
#------------------------------------------------------------------------------- |
|
# Save results |
|
#------------------------------------------------------------------------------- |
|
|
|
write.csv(anno_df, file.path(output_dir, "peaks_annotated.csv"), row.names = FALSE) |
|
write.csv(anno_df[anno_df$peak_type == "Promoter-proximal", ], |
|
file.path(output_dir, "peaks_promoter_proximal.csv"), row.names = FALSE) |
|
write.csv(anno_df[anno_df$peak_type == "Distal", ], |
|
file.path(output_dir, "peaks_distal.csv"), row.names = FALSE) |
|
|
|
#------------------------------------------------------------------------------- |
|
# Generate plots |
|
#------------------------------------------------------------------------------- |
|
|
|
pdf(file.path(output_dir, "annotation_piechart.pdf"), width = 10, height = 8) |
|
plotAnnoPie(peak_anno) |
|
dev.off() |
|
|
|
pdf(file.path(output_dir, "distance_to_tss.pdf"), width = 10, height = 6) |
|
plotDistToTSS(peak_anno, title = "Distribution of ATAC-seq Peaks Relative to TSS") |
|
dev.off() |
|
|
|
pdf(file.path(output_dir, "annotation_barplot.pdf"), width = 8, height = 6) |
|
plotAnnoBar(peak_anno) |
|
dev.off() |
|
|
|
cat("\nAnnotation complete. Results saved to:", output_dir, "\n") |
|
|
|
print(sessionInfo()) |
|
|
|
RSCRIPT |
|
|
|
#=============================================================================== |
|
# SECTION 4: SEQUENCE MODELING WITH GKM-SVM |
|
# Reference: https://github.com/Dongwon-Lee/lsgkm |
|
#=============================================================================== |
|
|
|
mkdir -p ${GKMSVM_DIR}/control |
|
mkdir -p ${GKMSVM_DIR}/veds |
|
|
|
#------------------------------------------------------------------------------- |
|
# 4A: Prepare sequences for gkm-SVM training |
|
#------------------------------------------------------------------------------- |
|
|
|
extract_summit_sequences() { |
|
local PEAKS=$1 |
|
local OUTPUT_PREFIX=$2 |
|
local WINDOW=150 |
|
|
|
awk -v w=${WINDOW} 'BEGIN{OFS="\t"} { |
|
summit = $2 + $10; # narrowPeak format: start + summit offset |
|
start = summit - w; |
|
end = summit + w; |
|
if (start < 0) start = 0; |
|
print $1, start, end, $4, $5, "+" |
|
}' ${PEAKS} > ${OUTPUT_PREFIX}_summits.bed |
|
|
|
bedtools getfasta \ |
|
-fi ${GENOME_FA} \ |
|
-bed ${OUTPUT_PREFIX}_summits.bed \ |
|
-fo ${OUTPUT_PREFIX}_pos.fa \ |
|
-name |
|
|
|
} |
|
|
|
extract_summit_sequences \ |
|
"${PEAK_DIR}/Control_1_peaks.filtered.narrowPeak" \ |
|
"${GKMSVM_DIR}/control/control" |
|
|
|
|
|
extract_summit_sequences \ |
|
"${PEAK_DIR}/VEDS_1_peaks.filtered.narrowPeak" \ |
|
"${GKMSVM_DIR}/veds/veds" |
|
|
|
#------------------------------------------------------------------------------- |
|
# 4B: Generate GC-matched negative sequences |
|
#------------------------------------------------------------------------------- |
|
|
|
# Function to generate nullseq (GC-matched negatives) |
|
# Requires genNullSeqs from lsgkm or equivalent |
|
generate_negatives() { |
|
local POS_FA=$1 |
|
local OUTPUT_FA=$2 |
|
local GENOME=$3 |
|
|
|
# Using gkmSVM's genNullSeqs |
|
genNullSeqs \ |
|
-i ${POS_FA} \ |
|
-g ${GENOME} \ |
|
-o ${OUTPUT_FA} \ |
|
-r 1 |
|
} |
|
|
|
generate_negatives \ |
|
"${GKMSVM_DIR}/control/control_pos.fa" \ |
|
"${GKMSVM_DIR}/control/control_neg.fa" \ |
|
"${GENOME_FA}" |
|
|
|
generate_negatives \ |
|
"${GKMSVM_DIR}/veds/veds_pos.fa" \ |
|
"${GKMSVM_DIR}/veds/veds_neg.fa" \ |
|
"${GENOME_FA}" |
|
|
|
#------------------------------------------------------------------------------- |
|
# 4C: Train gkm-SVM models |
|
# Parameters: l=11, k=7 |
|
#------------------------------------------------------------------------------- |
|
|
|
# Train Control model |
|
gkmtrain \ |
|
-l ${GKM_L} \ |
|
-k ${GKM_K} \ |
|
-T ${GKM_T} \ |
|
${GKMSVM_DIR}/control/control_pos.fa \ |
|
${GKMSVM_DIR}/control/control_neg.fa \ |
|
${GKMSVM_DIR}/control/control_model |
|
|
|
# Train VEDS model |
|
gkmtrain \ |
|
-l ${GKM_L} \ |
|
-k ${GKM_K} \ |
|
-T ${GKM_T} \ |
|
${GKMSVM_DIR}/veds/veds_pos.fa \ |
|
${GKMSVM_DIR}/veds/veds_neg.fa \ |
|
${GKMSVM_DIR}/veds/veds_model |
|
|
|
|
|
#=============================================================================== |
|
# SECTION 5: MOTIF ANALYSIS WITH GKMPWM |
|
# Reference: https://github.com/shigakiD/gkmPWM |
|
#=============================================================================== |
|
|
|
mkdir -p ${MOTIF_DIR}/control |
|
mkdir -p ${MOTIF_DIR}/veds |
|
|
|
#------------------------------------------------------------------------------- |
|
# 5A: Derive PWMs using gkmPWM lasso mode |
|
#------------------------------------------------------------------------------- |
|
|
|
# Run gkmPWM on Control model |
|
gkmPWMlasso \ |
|
-m ${GKMSVM_DIR}/control/control_model.model.txt \ |
|
-l ${GKM_L} \ |
|
-k ${GKM_K} \ |
|
-o ${MOTIF_DIR}/control/control_pwm \ |
|
--meme # Output in MEME format for downstream analysis |
|
|
|
# Run gkmPWM on VEDS model |
|
gkmPWMlasso \ |
|
-m ${GKMSVM_DIR}/veds/veds_model.model.txt \ |
|
-l ${GKM_L} \ |
|
-k ${GKM_K} \ |
|
-o ${MOTIF_DIR}/veds/veds_pwm \ |
|
--meme |
|
|
|
#------------------------------------------------------------------------------- |
|
# 5B: Match derived PWMs to known motif databases |
|
# - Redundancy (R): max Pearson correlation to other PWMs |
|
# - Weight (W): model contribution score |
|
# - Z-score (Z): standardized mean of highest-weighted gapped k-mers |
|
# - Importance (I): relative increase in model error when removing PWM |
|
#------------------------------------------------------------------------------- |
|
|
|
|
|
# Compare to MEME motif database |
|
MEME_DB="path/to/JASPAR2022_CORE_vertebrates_non-redundant_pfms_meme.txt" |
|
|
|
# Control motifs |
|
tomtom \ |
|
-o ${MOTIF_DIR}/control/tomtom_results \ |
|
-thresh 0.05 \ |
|
${MOTIF_DIR}/control/control_pwm.meme \ |
|
${MEME_DB} |
|
|
|
# VEDS motifs |
|
tomtom \ |
|
-o ${MOTIF_DIR}/veds/tomtom_results \ |
|
-thresh 0.05 \ |
|
${MOTIF_DIR}/veds/veds_pwm.meme \ |
|
${MEME_DB} |
|
|
|
|
|
#=============================================================================== |
|
# SECTION 6: TRANSCRIPTION FACTOR BINDING SITE PREDICTION |
|
#=============================================================================== |
|
|
|
mkdir -p ${MOTIF_DIR}/control/tfbs |
|
mkdir -p ${MOTIF_DIR}/veds/tfbs |
|
|
|
#------------------------------------------------------------------------------- |
|
# 6A: Map TF binding sites using mapTF |
|
# Parameters: same l, k as gkm-SVM |
|
# - Average k-mer posterior probability >= 0.90 |
|
# - Model motif correlation >= 0.80 |
|
#------------------------------------------------------------------------------- |
|
|
|
# Get top 2 TFs by importance score for each genotype |
|
|
|
# Map TFs on Control peaks |
|
mapTF \ |
|
-m ${GKMSVM_DIR}/control/control_model.model.txt \ |
|
-p ${MOTIF_DIR}/control/control_pwm.meme \ |
|
-s ${GKMSVM_DIR}/control/control_pos.fa \ |
|
-l ${GKM_L} \ |
|
-k ${GKM_K} \ |
|
-o ${MOTIF_DIR}/control/tfbs/control_tfbs.bed \ |
|
--prob_thresh 0.90 \ |
|
--corr_thresh 0.80 |
|
|
|
# Map TFs on VEDS peaks |
|
mapTF \ |
|
-m ${GKMSVM_DIR}/veds/veds_model.model.txt \ |
|
-p ${MOTIF_DIR}/veds/veds_pwm.meme \ |
|
-s ${GKMSVM_DIR}/veds/veds_pos.fa \ |
|
-l ${GKM_L} \ |
|
-k ${GKM_K} \ |
|
-o ${MOTIF_DIR}/veds/tfbs/veds_tfbs.bed \ |
|
--prob_thresh 0.90 \ |
|
--corr_thresh 0.80 |
|
|
|
#------------------------------------------------------------------------------- |
|
# 6B: Assign TFBS to nearest genes |
|
#------------------------------------------------------------------------------- |
|
|
|
# Export for R |
|
export MOTIF_DIR=${MOTIF_DIR} |
|
export ANNOT_DIR=${ANNOT_DIR} |
|
|
|
Rscript - <<'RSCRIPT' |
|
|
|
#------------------------------------------------------------------------------- |
|
# Assign TFBS to nearest genes |
|
# Peaks assigned by nearest-TSS linkage on mm10 |
|
# Labeled promoter-proximal if distance <= 3kb, distal otherwise |
|
#------------------------------------------------------------------------------- |
|
|
|
library(GenomicRanges) |
|
library(TxDb.Mmusculus.UCSC.mm10.knownGene) |
|
library(org.Mm.eg.db) |
|
library(ChIPseeker) |
|
|
|
motif_dir <- Sys.getenv("MOTIF_DIR", "path/to/motifs") |
|
annot_dir <- Sys.getenv("ANNOT_DIR", "path/to/annotation") |
|
|
|
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene |
|
|
|
# Function to annotate TFBS |
|
annotate_tfbs <- function(tfbs_file, output_prefix) { |
|
|
|
# Read TFBS BED file |
|
tfbs <- read.table(tfbs_file, header = FALSE, sep = "\t", |
|
col.names = c("chr", "start", "end", "tf_name", "score", "strand")) |
|
|
|
if (nrow(tfbs) == 0) { |
|
cat(" No TFBS found in", tfbs_file, "\n") |
|
return(NULL) |
|
} |
|
|
|
# Create GRanges |
|
tfbs_gr <- GRanges( |
|
seqnames = tfbs$chr, |
|
ranges = IRanges(start = tfbs$start, end = tfbs$end), |
|
tf_name = tfbs$tf_name, |
|
score = tfbs$score |
|
) |
|
|
|
# Annotate with nearest gene |
|
tfbs_anno <- annotatePeak( |
|
tfbs_gr, |
|
tssRegion = c(-3000, 3000), |
|
TxDb = txdb, |
|
annoDb = "org.Mm.eg.db", |
|
verbose = FALSE |
|
) |
|
|
|
tfbs_df <- as.data.frame(tfbs_anno) |
|
|
|
# Add promoter-proximal vs distal classification |
|
tfbs_df$peak_type <- ifelse( |
|
abs(tfbs_df$distanceToTSS) <= 3000, |
|
"Promoter-proximal", |
|
"Distal" |
|
) |
|
|
|
# Save |
|
write.csv(tfbs_df, paste0(output_prefix, "_annotated.csv"), row.names = FALSE) |
|
|
|
cat(" Annotated", nrow(tfbs_df), "TFBS\n") |
|
cat(" Promoter-proximal:", sum(tfbs_df$peak_type == "Promoter-proximal"), "\n") |
|
cat(" Distal:", sum(tfbs_df$peak_type == "Distal"), "\n") |
|
|
|
return(tfbs_df) |
|
} |
|
|
|
# Annotate Control TFBS |
|
cat("\nAnnotating Control TFBS...\n") |
|
control_tfbs <- annotate_tfbs( |
|
file.path(motif_dir, "control/tfbs/control_tfbs.bed"), |
|
file.path(motif_dir, "control/tfbs/control_tfbs") |
|
) |
|
|
|
# Annotate VEDS TFBS |
|
cat("\nAnnotating VEDS TFBS...\n") |
|
veds_tfbs <- annotate_tfbs( |
|
file.path(motif_dir, "veds/tfbs/veds_tfbs.bed"), |
|
file.path(motif_dir, "veds/tfbs/veds_tfbs") |
|
) |
|
|
|
cat("\nTFBS annotation complete.\n") |
|
|
|
RSCRIPT |
|
|
|
#=============================================================================== |
|
# SECTION 7: COMPARATIVE ANALYSIS (CONTROL VS VEDS) |
|
#=============================================================================== |
|
|
|
export MOTIF_DIR=${MOTIF_DIR} |
|
export ANNOT_DIR=${ANNOT_DIR} |
|
export PEAK_DIR=${PEAK_DIR} |
|
|
|
Rscript - <<'RSCRIPT' |
|
#------------------------------------------------------------------------------- |
|
# Compare enriched motifs and target genes between Control and VEDS |
|
#------------------------------------------------------------------------------- |
|
|
|
library(dplyr) |
|
library(ggplot2) |
|
|
|
motif_dir <- Sys.getenv("MOTIF_DIR", "path/to/motifs") |
|
annot_dir <- Sys.getenv("ANNOT_DIR", "path/to/annotation") |
|
|
|
output_dir <- file.path(annot_dir, "comparative") |
|
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE) |
|
|
|
#------------------------------------------------------------------------------- |
|
# Load TFBS annotations |
|
#------------------------------------------------------------------------------- |
|
|
|
control_tfbs <- read.csv(file.path(motif_dir, "control/tfbs/control_tfbs_annotated.csv")) |
|
veds_tfbs <- read.csv(file.path(motif_dir, "veds/tfbs/veds_tfbs_annotated.csv")) |
|
|
|
#------------------------------------------------------------------------------- |
|
# Compare TF enrichment |
|
#------------------------------------------------------------------------------- |
|
|
|
cat("=== TF Binding Site Summary ===\n") |
|
|
|
cat("\nControl - Top TFs:\n") |
|
control_tf_summary <- control_tfbs %>% |
|
group_by(tf_name) %>% |
|
summarise( |
|
n_sites = n(), |
|
n_promoter = sum(peak_type == "Promoter-proximal"), |
|
n_distal = sum(peak_type == "Distal"), |
|
n_genes = n_distinct(SYMBOL) |
|
) %>% |
|
arrange(desc(n_sites)) |
|
|
|
print(head(control_tf_summary, 10)) |
|
|
|
cat("\nVEDS - Top TFs:\n") |
|
veds_tf_summary <- veds_tfbs %>% |
|
group_by(tf_name) %>% |
|
summarise( |
|
n_sites = n(), |
|
n_promoter = sum(peak_type == "Promoter-proximal"), |
|
n_distal = sum(peak_type == "Distal"), |
|
n_genes = n_distinct(SYMBOL) |
|
) %>% |
|
arrange(desc(n_sites)) |
|
|
|
print(head(veds_tf_summary, 10)) |
|
|
|
#------------------------------------------------------------------------------- |
|
# Save summary tables |
|
#------------------------------------------------------------------------------- |
|
|
|
write.csv(control_tf_summary, file.path(output_dir, "control_tf_summary.csv"), row.names = FALSE) |
|
write.csv(veds_tf_summary, file.path(output_dir, "veds_tf_summary.csv"), row.names = FALSE) |
|
|
|
cat("\nComparative analysis complete. Results saved to:", output_dir, "\n") |
|
|
|
RSCRIPT |
