Scripts summary

To reuse some scripts, here is a summary for myself organized by topics

#3D genomics

1. from .paired.gz to hic

#!/usr/bin/env bash
set -euo pipefail

# ─── User settings ──────────────────────────────────────────────────────────
data_dir="../data"
sample="HAP-WT"

chromsizes="../annotations/hg38.chrom.sizes"
build="hg38"
juicer_tools="../tools/juicer_tools.2.20.00.jar"

merged_pairs="${data_dir}/${sample}_330M.pairs.gz"
sorted_merged="${merged_pairs%.pairs.gz}.sorted.pairs.gz"
dedup_pairs="${data_dir}/${sample}_330M.nodups.pairs.gz"
sorted_dedup="${dedup_pairs%.pairs.gz}.sorted.pairs.gz"
pre_pairs="${data_dir}/${sample}_330M.forPre.pairs.gz"
hic_file="${data_dir}/${sample}_330M.hic"

# comma‐sep list of resolutions (no spaces!)
res_list="2500000,1000000,500000,250000,100000,50000,25000,10000,5000,1000"

threads=12
# ────────────────────────────────────────────────────────────────────────────

echo "2 sort merged pairs"
pairtools sort --nproc $threads -o "$sorted_merged" "$merged_pairs"

echo "3 global deduplication"
pairtools dedup --nproc-in  $threads --nproc-out $threads --output "$dedup_pairs" "$sorted_merged"

echo "4 sort deduplicated pairs"
pairtools sort --nproc $threads -o "$sorted_dedup" "$dedup_pairs"

echo "5 prepare for juicebox pre:"
gunzip -c "$sorted_dedup" | grep -E '^(## pairs format|^#columns|^[^#])' | gzip -c > "$pre_pairs"

echo "6 building multi-res .hic"
java -Xmx16g -jar "$juicer_tools" pre -t $threads -r "$res_list" "$pre_pairs" "$hic_file" "$build"

#Binformatics

1. change bedgraph to tdf for igv visualization

# conda activate tools_igvtools

genome="hg19"
src_dir="../3_prepare_macs2_peaks_forIGV"

for bg in "$src_dir"/*.bedgraph; do
  base="${bg%.bedgraph}"
  echo ${base}
  tdf=${base}.tdf

  igvtools toTDF ${bg} ${tdf} ${genome}
done

2. draw plot using Python

import pathlib as p
import numpy as np
import pandas as pd

import matplotlib
# matplotlib.use('Qt5Agg')
# matplotlib.rcParams['font.family'] = 'Arial'
import matplotlib.pyplot as plt
plt.style.use('ggplot')

from scipy.stats import pearsonr, spearmanr, gaussian_kde

pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
pd.set_option('display.width', 500)

if __name__ == '__main__':
    root = p.Path(__file__).absolute().parent

    resolution = 10000
    src_dir = root.parent / f'data_2_cool_{resolution}_targetChroms_normalizeSmallest'
    dest_dir = root.parent / 'result_scatterPlot_of_twoReplicates_allChroms_noLabels_withColorCode'
    dest_dir.mkdir(parents=True, exist_ok=True)

    groups = [['ctrl1', 'ctrl2', 'ctrl', '#4D7731']]

    for group in groups:
        rep1 = group[0]
        rep2 = group[1]
        group_name = group[2]
        color_code = group[3]
        print(group_name)
        print(rep1)
        print(rep2)
        c1 = cooler.Cooler(str(src_dir / f'{rep1}_{resolution}_targetChroms_normalizeSmallest.cool'))
        c2 = cooler.Cooler(str(src_dir / f'{rep2}_{resolution}_targetChroms_normalizeSmallest.cool'))

        matrix1_flat_raw = c1.matrix(balance=False)[:].flatten()
        matrix2_flat_raw = c2.matrix(balance=False)[:].flatten()

        matrix1_flat, matrix2_flat = filter_zero_indices(matrix1_flat_raw, matrix2_flat_raw)

        corr_pearson, p_pearson = pearsonr(matrix1_flat, matrix2_flat)

        fig, ax = plt.subplots(figsize=(5, 5), dpi=300)
        ax.set_facecolor('white')
        for spine in ax.spines.values():
            spine.set_color('black')
            spine.set_linewidth(1)

        plt.scatter(np.log10(matrix1_flat), np.log10(matrix2_flat), alpha=0.3, s=1.5, c=color_code)
        # plt.scatter(np.log2(matrix1_flat), np.log2(matrix2_flat), alpha=0.3, s=1.5, c='black')
        # sc = ax.scatter(x, y, c=z, s=1.5, alpha=0.3, cmap='inferno')

        ax.set_xticks([1, 2, 3])
        ax.set_xticklabels([10, 100, 1000])
        ax.set_yticks([1, 2, 3])
        ax.set_yticklabels([10, 100, 1000])

        # ax.set_xlabel(f"{group_name}_rep1", fontsize=10)
        # ax.set_ylabel(f"{group_name}_rep2", fontsize=10)
        ax.set_xlabel(f"")
        ax.set_ylabel(f"")

        ax.set_title(f"")

        plt.savefig(dest_dir/f'{group_name}_resolution{resolution}_correlation_twoReps_allChroms_pearsonCorr{corr_pearson:.2f}_noLabels_withColorCode.png', dpi=300)