ATAC-seq Cross-Species
Comparative Analysis

# Decompress mouse peaks
cd "${BASE}/MouseAtac/AdrenalGland/peak/idr_reproducibility"
zcat idr.conservative_peak.narrowPeak.gz > idr.conservative_peak.narrowPeak

# Decompress human peaks
cd "${BASE}/HumanAtac/peak/idr_reproducibility"
zcat idr.conservative_peak.narrowPeak.gz > idr.conservative_peak.narrowPeak

bash "$SCRIPT_DIR/integrate_halper.sh" \
  --base       "$SCRIPT_DIR" \
  --hal-file   "$HAL_FILE" \
  --halper-map "$HALPER_MAP"

export PATH="${HAL_BIN}:${PATH}"
export PYTHONPATH="${HALPER_PP}:${PYTHONPATH:-}"
mkdir -p "${HALPER_OUT}"

bash "${HALPER_MAP_SH}" \
  -b "${MOUSE_CONSERVATIVE_NARROWPEAK}" \  # input mouse peaks
  -o "${HALPER_OUT}/"                    \  # output directory
  -s Mouse                              \  # source species
  -t Human                              \  # target species
  -c "${HAL_FILE}"                       # whole-genome alignment

# Sort both files by chromosome then start position
sort -k1,1 -k2,2n mouse_to_human_conservative.narrowPeak \
     > mouse_to_human_conservative.sorted.narrowPeak

sort -k1,1 -k2,2n "${HUMAN_CONSERVATIVE_NARROWPEAK}" \
     > human_conservative.sorted.narrowPeak

# Shared: mouse peaks that overlap human peaks
bedtools intersect \
  -a mouse_to_human_conservative.sorted.narrowPeak \
  -b human_conservative.sorted.narrowPeak \
  -u > shared_peaks_conservative.narrowPeak

# Mouse-specific: lifted mouse peaks NOT in human
bedtools intersect \
  -a mouse_to_human_conservative.sorted.narrowPeak \
  -b human_conservative.sorted.narrowPeak \
  -v > mouse_specific_peaks_conservative.narrowPeak

# Human-specific: human peaks NOT in lifted mouse
bedtools intersect \
  -a human_conservative.sorted.narrowPeak \
  -b mouse_to_human_conservative.sorted.narrowPeak \
  -v > human_specific_peaks_conservative.narrowPeak

awk 'BEGIN{OFS="\t"} {
    chr=$1; start=$2; end=$3
    name=$4; score=$5; strand=$6; summit=$10
    if (summit != -1 && summit != ".") {
        center = start + summit
        print chr, center-1, center, name, score, \
              (strand=="." ? "+" : strand)
    } else {
        print chr, start, end, name, score, \
              (strand=="." ? "+" : strand)
    }
}' "$NARROWPEAK_FILE" > "$BED_FILE"

# [1/2] Find enriched motifs
findMotifsGenome.pl \
    "$BED_FILE" "$GENOME" "$MOTIF_DIR" \
    -size 200 -mask

# [2/2] Annotate peaks with motif instances
annotatePeaks.pl \
    "$BED_FILE" "$GENOME" \
    -m  "${MOTIF_DIR}/nonRedundant.motifs" \
    -mbed "${SAMPLE_DIR}/motif_instances.bed" \
    -size 200 \
    > "$ANNO_FILE"

awk -F'\t' 'BEGIN {OFS="\t"}
NR==1 {
    # Detect column positions dynamically from header
    for (i=2; i<=NF; i++) {
        if ($i=="Annotation")      annot = i
        if ($i=="Distance to TSS") dist  = i
        if ($i ~ /Distance From Peak/) motif_cols[n_motifs++] = i
    }
    # Build output header: PeakID + Annotation + motif cols
    header = "PeakID\tAnnotation"
    for (m=0; m<n_motifs; m++) header = header "\t" $motif_cols[m]
    print header; next
}
{
    new_annot = ""
    # Enhancer: intergenic/intronic AND > 2000 bp from TSS, or labeled enhancer
    if (((tolower($annot) ~ /intergenic/ || tolower($annot) ~ /intron/)
         && $dist >= 2000) || (tolower($annot) ~ /enhancer/))
        new_annot = "enhancer"
    # Promoter: labeled "promoter"
    else if (tolower($annot) ~ /promoter/)
        new_annot = "promoter"
    # Only print classified peaks
    if (new_annot != "") { ... print ... }
}'

python enhancer_promoter_analysis.py \
  --input-dir  ./filtered_annotations \
  --output-dir ./motif_annotation_split

bash run_GO_pipeline.sh \
  ./filtered_annotations \   # directory scanned for the 4 expected peak files
  ./rGREAT_results           # output directory (default: rGREAT_results)

# Read narrowPeak as GRanges (3-col BED, standard chroms only)
gr <- read_bed3(file)

# GO:BP enrichment
job_bp <- great(gr = gr, tss_source = genome, gene_sets = "GO:BP")
tbl_bp <- getEnrichmentTable(job_bp)

# Also run GO:MF and GO:CC
job_mf <- great(gr = gr, tss_source = genome, gene_sets = "GO:MF")
job_cc <- great(gr = gr, tss_source = genome, gene_sets = "GO:CC")

# Publication dot-plot (top 20 terms)
p <- ggplot(plot_tbl, aes(x = fold_enrichment, y = description,
          size = observed_region_hits, fill = log10_p_adjust)) +
  geom_point(alpha = 0.9, shape = 21) +
  scale_fill_distiller(palette = "YlOrRd", direction = 1) +
  theme_classic(base_size = 14)

cd RUN_FULL_PIPELINE

sbatch complete_analysis_pipeline.sh \
  /path/to/10plusway-master.hal \            # ARG 1: .hal alignment file
  /path/to/halper_map_peak_orthologs.sh      # ARG 2: HALPER mapping script