This page generates analysis code in your browser. It does not upload data or execute the workflow on the Dormanomics server. Review the generated parameters, software references and input requirements, then run and validate the code in your own computing environment.
#!/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")