#!/bin/bash
# Step 1.1: MACS2 peak calling (ChIP-seq/ATAC-seq)
mkdir -p 01_peaks
macs2 callpeak -t sample_treatment.bam -c sample_input.bam \
-f BAM -g {genome} -n 01_peaks/treatment \
--broad --broad-cutoff 0.1 --qvalue 0.05
# For ATAC-seq:
# macs2 callpeak -t sample_atac.bam -f BAMPE -g mm -n 01_peaks/atac --nomodel --shift -100 --extsize 200
# Step 1.2: Differential binding with DiffBind (R)
Rscript -e '
library(DiffBind)
samples <- dba(sampleSheet="{samplesheet}")
samples <- dba.count(samples, summits=250)
samples <- dba.normalize(samples)
samples <- dba.contrast(samples, categories=DBA_CONDITION)
samples <- dba.analyze(samples)
pdf("02_diffbind_heatmap.pdf")
dba.plotHeatmap(samples)
dev.off()
deg_peaks <- dba.report(samples)
write.csv(as.data.frame(deg_peaks), "02_differential_peaks.csv")
'
# Step 1.3: MethylKit for WGBS
Rscript -e '
library(methylKit)
file_list <- list("C1.txt","C2.txt","T1.txt","T2.txt")
obj <- methRead(file_list, sample.id=list("C1","C2","T1","T2"),
assembly="mm10", treatment=c(0,0,1,1), context="CpG")
obj_filt <- filterByCoverage(obj, lo.count=10, lo.perc=NULL, hi.count=NULL, hi.perc=99.9)
obj_norm <- normalizeCoverage(obj_filt)
obj_meth <- unite(obj_norm, destrand=FALSE)
myDiff <- calculateDiffMeth(obj_meth)
dmrs <- getMethylDiff(myDiff, difference=25, qvalue=0.01)
write.csv(as.data.frame(dmrs), "02_DMRs.csv")
library(chromVAR)
library(JASPAR2020)
library(motifmatchr)
library(ggplot2)
# Step 2.1: Read peaks and count fragments
peaks <- getPeaks("{peaks_bed}")
counts <- getCounts("{bigwig}", peaks, paired=FALSE)
# Step 2.2: Compute deviations
motifs <- getJasparMotifs(species=9606, collection="CORE")
motif_ix <- matchMotifs(motifs, peaks, genome=BSgenome.Mmusculus.UCSC.mm10)
dev <- computeDeviations(object=counts, annotations=motif_ix)
# Step 2.3: TF-specific footprinting
target_tfs <- unlist(strsplit("{target_tfs}", ","))
plot_list <- list()
for(tf in target_tfs) {
if(tf %in% names(motifs)) {
p <- plotFootprint(dev, motifs[[tf]], window=500) +
ggtitle(paste(tf, "Footprint"))
plot_list[[tf]] <- p
}
}
ggsave("03_tf_footprints.pdf", marrangeGrob(plot_list, nrow=2, ncol=3),
width=15, height=10)
# Step 2.4: Motif enrichment
if(require("monaLisa")) {
se <- calcBinnedMotifEnrR(seqs=getSeq(BSgenome.Mmusculus.UCSC.mm10, peaks),
bins=factor(rep(c("open","closed"), each=length(peaks)/2)),
pwm=motifs)
plotMotifHeatmap(se, show_motif_GC=TRUE)
ggsave("03_motif_enrichment.pdf", width=10, height=8)
}library(ChIPseeker)
library(GenomicRanges)
library(ggplot2)
library(clusterProfiler)
# Step 3.1: Annotate peaks to genes
peaks <- readPeakFile("{peak_regions}")
txdb <- makeTxDbFromGFF("{gene_annot}", format="gff3")
peak_annot <- annotatePeak(peaks, TxDb=txdb, annoDb="org.Mm.eg.db",
addFlankGeneInfo=TRUE, flankDistance={link_dist})
write.csv(as.data.frame(peak_annot), "04_peak_annotation.csv")
# Step 3.2: Plot peak distribution
plotAnnoPie(peak_annot)
ggsave("04_peak_distribution.pdf", width=8, height=6)
# Step 3.3: Correlate with DEGs
degs <- read.csv("{rna_deg}")
peak_df <- as.data.frame(peak_annot)
linked <- merge(peak_df, degs, by.x="geneId", by.y="gene_id")
# 4-quadrant plot: ATAC signal vs RNA fold-change
ggplot(linked, aes(x=log2FoldChange_ATAC, y=log2FoldChange_RNA)) +
geom_point(aes(color=annotation), alpha=0.5, size=1) +
geom_vline(xintercept=0, linetype="dashed") +
geom_hline(yintercept=0, linetype="dashed") +
theme_bw(base_size=14) +
labs(title="ATAC vs RNA Changes",
x="ATAC Log2FC", y="RNA Log2FC")
ggsave("04_epi_rna_correlation.pdf", width=8, height=6)
# Step 3.4: Chromatin switch analysis
# Categorize genes by epigenetic-transcription state
linked$chrom_state <- case_when(
linked$ATAC_logFC > 0 & linked$log2FoldChange > 0 ~ "Active_Active",
linked$ATAC_logFC > 0 & linked$log2FoldChange < 0 ~ "Active_Repressed",
linked$ATAC_logFC < 0 & linked$log2FoldChange > 0 ~ "Repressed_Active",
TRUE ~ "Stable"
)
write.csv(linked, "04_chromatin_switch.csv")library(methylKit)
library(ggplot2)
library(tidyverse)
# Step 4.1: Load multi-stage methylation
stages <- unlist(strsplit("{stages}", ","))
target_genes <- unlist(strsplit("{target_genes}", ","))
obj <- readRDS("{methyl_obj}")
obj_united <- unite(obj, destrand=FALSE)
# Step 4.2: Extract promoter methylation for target genes
promoter_meth <- getMethylationStats(obj_united, plot=FALSE)
promoter_df <- as.data.frame(promoter_meth)
# Step 4.3: Methylation trajectory
traj_df <- promoter_df %>%
filter(gene %in% target_genes) %>%
gather(key="stage", value="methylation", -gene, -chr, -start, -end) %>%
mutate(stage=factor(stage, levels=stages))
ggplot(traj_df, aes(x=stage, y=methylation, color=gene, group=gene)) +
geom_line(linewidth=1) + geom_point(size=3) +
theme_bw(base_size=14) +
labs(title="Promoter Methylation Trajectory",
x="Dormancy Stage", y="Methylation Level (%)")
ggsave("05_methylation_trajectory.pdf", width=10, height=5)
# Step 4.4: Reversibility score
# Compare Torpor vs Post to check if marks return to baseline
reversibility <- traj_df %>%
spread(key="stage", value="methylation") %>%
mutate(reversibility_score = abs(Torpor - Pre) - abs(Post - Pre))
write.csv(reversibility, "05_reversibility_scores.csv")
# Heatmap of all reversible marks
rev_matrix <- reversibility %>%
select(gene, reversibility_score) %>%
column_to_rownames("gene")
pheatmap::pheatmap(as.matrix(rev_matrix),
main="Epigenetic Reversibility Score",
filename="05_reversibility_heatmap.pdf")