scientific-epigenomics-chromatin

name	scientific-epigenomics-chromatin
description	エピゲノミクス・クロマチン生物学解析スキル。ChIP-seq ピーク呼び出し (MACS2/MACS3)、 ATAC-seq ヌクレオソームフリー領域検出、DNA メチル化パターン解析 (WGBS/RRBS)、 ヒストン修飾クロマチン状態モデリング (ChromHMM)、Hi-C 接触マップ・TAD 検出、 転写因子結合サイト予測 (モチーフ濃縮)、差次結合解析 (DiffBind) を統合した 計算エピゲノミクスパイプライン。ChIP-Atlas 43 万+実験との連携対応。 ToolUniverse 連携: chipatlas。
tu_tools	[{"key":"chipatlas","name":"ChIP-Atlas","description":"ChIP-Atlas エピゲノミクスエンリッチメント解析 (43万+実験)"}]

Scientific Epigenomics & Chromatin Biology

ChIP-seq・ATAC-seq・バイサルファイトシーケンシング・Hi-C データを対象に、 ピーク呼び出し→差次結合解析→クロマチン状態注釈→3D ゲノム構造解析の 統合エピゲノミクスパイプラインを提供する。

When to Use

• ChIP-seq データからヒストン修飾・転写因子結合部位を同定するとき
• ATAC-seq でクロマチンアクセシビリティを評価するとき
• DNA メチル化（WGBS/RRBS）パターンを解析するとき
• Hi-C データから TAD/ループ・3D ゲノム構造を推定するとき
• 複数エピゲノムマークを統合してクロマチン状態を分類するとき

---

Quick Start

1. ChIP-seq ピーク呼び出し (MACS2/MACS3)

import subprocess
import pandas as pd
import numpy as np


def chipseq_peak_calling(treatment_bam, control_bam, genome_size="hs",
                         outdir="results/chipseq", name="sample",
                         peak_type="narrow", qvalue=0.05):
    """
    MACS2/MACS3 による ChIP-seq ピーク呼び出し。

    Parameters:
        treatment_bam: 処理群 BAM ファイル
        control_bam: コントロール BAM ファイル (Input/IgG)
        genome_size: 有効ゲノムサイズ (hs/mm/ce/dm or int)
        peak_type: "narrow" (TF) or "broad" (ヒストン修飾 H3K27me3 等)
        qvalue: FDR 閾値
    """
    import os
    os.makedirs(outdir, exist_ok=True)

    cmd = [
        "macs3", "callpeak",
        "-t", treatment_bam,
        "-c", control_bam,
        "-g", str(genome_size),
        "--outdir", outdir,
        "-n", name,
        "-q", str(qvalue),
        "--keep-dup", "auto",
        "--call-summits",
    ]

    if peak_type == "broad":
        cmd.extend(["--broad", "--broad-cutoff", str(qvalue)])

    print(f"Running MACS3 peak calling ({peak_type} mode)...")
    subprocess.run(cmd, check=True)

    # ピークファイル読み込み
    suffix = "broadPeak" if peak_type == "broad" else "narrowPeak"
    peak_file = f"{outdir}/{name}_peaks.{suffix}"
    cols = ["chr", "start", "end", "name", "score", "strand",
            "signalValue", "pValue", "qValue"]
    if peak_type == "narrow":
        cols.append("summit")

    peaks = pd.read_csv(peak_file, sep="\t", header=None, names=cols)
    peaks["width"] = peaks["end"] - peaks["start"]

    print(f"  Called {len(peaks):,} {peak_type} peaks (q < {qvalue})")
    print(f"  Median peak width: {peaks['width'].median():.0f} bp")
    print(f"  Mean signal value: {peaks['signalValue'].mean():.2f}")

    return peaks


def chipseq_qc_metrics(peaks, frip_bam=None, total_reads=None):
    """
    ChIP-seq QC 指標の算出。

    Returns:
        dict: peak 数、中央値幅、FRiP (Fraction of Reads in Peaks)
    """
    metrics = {
        "n_peaks": len(peaks),
        "median_width_bp": float(peaks["width"].median()),
        "mean_signal": float(peaks["signalValue"].mean()),
        "mean_log10_qvalue": float(peaks["qValue"].mean()),
    }

    # ENCODE 品質基準との比較
    if metrics["n_peaks"] < 500:
        metrics["quality_flag"] = "LOW — < 500 peaks"
    elif metrics["n_peaks"] < 10000:
        metrics["quality_flag"] = "MODERATE"
    else:
        metrics["quality_flag"] = "HIGH"

    return metrics

2. ATAC-seq アクセシビリティ解析

import numpy as np
import pandas as pd


def atacseq_nucleosome_free_regions(fragments_file, output_dir="results/atacseq"):
    """
    ATAC-seq フラグメントサイズ分布に基づくヌクレオソーム占有解析。

    フラグメントサイズによる分類:
    - < 150 bp: Nucleosome-Free Region (NFR)
    - 150-300 bp: Mono-nucleosome
    - 300-500 bp: Di-nucleosome
    - > 500 bp: Tri-nucleosome+
    """
    import os
    os.makedirs(output_dir, exist_ok=True)

    # フラグメントサイズ分布
    fragments = pd.read_csv(fragments_file, sep="\t",
                            names=["chr", "start", "end", "barcode", "count"])
    fragments["length"] = fragments["end"] - fragments["start"]

    # サイズカテゴリ分類
    bins = [0, 150, 300, 500, 10000]
    labels = ["NFR (<150)", "Mono-nuc (150-300)",
              "Di-nuc (300-500)", "Tri-nuc+ (>500)"]
    fragments["category"] = pd.cut(fragments["length"], bins=bins, labels=labels)

    size_dist = fragments["category"].value_counts(normalize=True)
    nfr_ratio = size_dist.get("NFR (<150)", 0)

    print(f"  Fragment size distribution:")
    for cat, pct in size_dist.items():
        print(f"    {cat}: ")
    ()

     fragments, size_dist


 ():
    
     pybedtools  BedTool

    peaks_bt = BedTool.from_dataframe(
        peaks[[, , , , ]]
    )
    
    

    ()
    ()

     peaks_bt

3. DNA メチル化パターン解析

import numpy as np
import pandas as pd


def bisulfite_methylation_analysis(methylation_file, min_coverage=10,
                                    output_prefix="results/methylation"):
    """
    WGBS/RRBS バイサルファイトシーケンシングデータのメチル化解析。

    入力: Bismark methylation extractor 出力 (CpG context)
    処理:
    1. カバレッジフィルタリング
    2. メチル化レベル算出 (β 値)
    3. CpG アイランド/ショア/シェルフ注釈
    4. 差次メチル化領域 (DMR) 検出
    """
    import os
    os.makedirs(os.path.dirname(output_prefix), exist_ok=True)

    # Bismark 出力の読み込み
    df = pd.read_csv(methylation_file, sep="\t",
                     names=["chr", "pos", "strand", "count_m", "count_u"])
    df["coverage"] = df["count_m"] + df["count_u"]
    df["beta"] = df["count_m"] / df["coverage"]

    # カバレッジフィルタ
    n_before = len(df)
    df = df[df["coverage"] >= min_coverage].copy()
    print(f"  Coverage filter (≥{min_coverage}x): {n_before:,} → {len(df):,} CpGs")

    # グローバルメチル化統計
    mean_beta = df["beta"].mean()
    median_beta = df["beta"].median()
    print(f"  Global methylation: mean β = {mean_beta:.3f}, median β = {median_beta:.3f}")

    # メチル化状態分類
    df[] = pd.cut(df[],
                          bins=[, , , ],
                          labels=[, , ])
    status_counts = df[].value_counts(normalize=)
    ()
    ()
    ()

     df


 ():
    
     scipy.stats  mannwhitneyu

    results = []
    mean_g1 = group1_betas.mean(axis=)
    mean_g2 = group2_betas.mean(axis=)
    delta_beta = mean_g2 - mean_g1

     i  ((positions)):
        stat, pval = mannwhitneyu(
            group1_betas[i, :], group2_betas[i, :], alternative=
        )
        results.append({
            : positions[i][],
            : positions[i][],
            : (delta_beta[i]),
            : pval,
            : (mean_g1[i]),
            : (mean_g2[i]),
        })

    df = pd.DataFrame(results)

    
     statsmodels.stats.multitest  multipletests
    df[] = multipletests(df[], method=)[]

    
    sig = df[(df[] < pvalue_cutoff) &
             (df[].() >= delta_beta_cutoff)]
    ()

     df, sig

4. クロマチン状態モデリング (ChromHMM)

import subprocess
import pandas as pd
import numpy as np


def chromhmm_learn_model(binarized_dir, output_dir, n_states=15,
                         assembly="hg38"):
    """
    ChromHMM によるクロマチン状態モデリング。

    複数のヒストン修飾マーク (H3K4me1/me3, H3K27ac, H3K27me3,
    H3K36me3, H3K9me3 等) を入力として、ゲノムをクロマチン状態に分類。

    Roadmap Epigenomics 15-state モデル:
    1-TssA, 2-TssAFlnk, 3-TxFlnk, 4-Tx, 5-TxWk,
    6-EnhG, 7-Enh, 8-ZNF/Rpts, 9-Het, 10-TssBiv,
    11-BivFlnk, 12-EnhBiv, 13-ReprPC, 14-ReprPCWk, 15-Quies
    """
    import os
    os.makedirs(output_dir, exist_ok=True)

    cmd = [
        "java", "-mx8G", "-jar", "ChromHMM.jar", "LearnModel",
        "-b", "200",
        binarized_dir, output_dir, str(n_states), assembly
    ]
    print(f"Running ChromHMM LearnModel with {n_states} states...")
    subprocess.run(cmd, check=True)

    # 遷移確率マトリクス読み込み
    trans_file = f"{output_dir}/transitions_{n_states}.txt"
    if os.path.exists(trans_file):
        trans = pd.read_csv(trans_file, sep="\t", index_col=0)
        print(f"  Transition matrix: {trans.shape}")

    # エミッション確率読み込み
    emit_file = f"{output_dir}/emissions_{n_states}.txt"
    if os.path.exists(emit_file):
        emit = pd.read_csv(emit_file, sep="\t", index_col=)
        ()

     {: n_states, : output_dir}


 ():
    
    default_labels = {
        : , : ,
        : , : ,
        : , : ,
        : , : ,
        : , : ,
        : , : ,
        : , : ,
        : ,
    }
    labels = state_labels  default_labels

    segments = pd.read_csv(segments_bed, sep=,
                           names=[, , , ])
    segments[] = segments[] - segments[]
    segments[] = segments[].(labels)

    
    total_bp = segments[].()
    state_coverage = segments.groupby()[].() / total_bp
    ()
     label, pct  state_coverage.sort_values(ascending=).items():
        ()

     segments, state_coverage

5. Hi-C 3D ゲノム構造解析

import numpy as np
import pandas as pd


def hic_contact_matrix_analysis(cool_file, resolution=10000,
                                 chromosome="chr1"):
    """
    Hi-C 接触マップ解析 (.cool/.mcool 形式)。

    1. ICE 正規化
    2. A/B コンパートメント同定 (PCA)
    3. TAD 呼び出し (Insulation Score)
    """
    import cooler

    # クールファイル読み込み
    clr = cooler.Cooler(f"{cool_file}::resolutions/{resolution}")
    matrix = clr.matrix(balance=True).fetch(chromosome)

    print(f"  Contact matrix shape: {matrix.shape}")
    print(f"  Resolution: {resolution:,} bp")
    print(f"  Non-zero entries: {np.count_nonzero(~np.isnan(matrix)):,}")

    return matrix


def call_tads_insulation_score(matrix, resolution=10000, window_size=500000):
    """
    Insulation Score 法による TAD (Topologically Associating Domain) 呼び出し。

    Parameters:
        window_size: Insulation window サイズ (bp)
    """
    window_bins = window_size // resolution

    n = matrix.shape[0]
    insulation = np.zeros(n)

    for i in range(window_bins, n - window_bins):
        submat = matrix[i - window_bins:i, i:i + window_bins]
        insulation[i] = np.nanmean(submat)

    # log2 正規化
    mean_val = np.nanmean(insulation[insulation > 0])
    log_insulation = np.log2(insulation / mean_val + 1e-10)

    # TAD 境界 = Insulation Score の極小値
     scipy.signal  argrelextrema
    minima = argrelextrema(log_insulation, np.less, order=)[]

    tad_boundaries = minima * resolution
    n_tads = (tad_boundaries) - 

    ()
    ()
    ()

     log_insulation, tad_boundaries


 ():
    
     sklearn.decomposition  PCA

    
    matrix_clean = np.nan_to_num(matrix, nan=)
    expected = np.zeros_like(matrix_clean)
     d  (matrix_clean.shape[]):
        diag_vals = np.diag(matrix_clean, d)
        mean_val = np.mean(diag_vals)  (diag_vals) >   
        np.fill_diagonal(expected[d:, :], mean_val)
        np.fill_diagonal(expected[:, d:], mean_val)

    oe_matrix = matrix_clean / (expected + )

    
    corr_matrix = np.corrcoef(oe_matrix)
    corr_matrix = np.nan_to_num(corr_matrix)

    pca = PCA(n_components=)
    components = pca.fit_transform(corr_matrix)
    pc1 = components[:, ]

    
    compartment = np.where(pc1 > , , )
    a_frac = np.mean(compartment == )

    ()
    ()
    ()

     pc1, compartment

6. 転写因子モチーフ濃縮解析

import pandas as pd
import numpy as np
from scipy.stats import fisher_exact


def motif_enrichment_analysis(peak_sequences, background_sequences,
                               jaspar_db="JASPAR2024_CORE_vertebrates",
                               pvalue_cutoff=0.01):
    """
    ピーク領域における転写因子結合モチーフの濃縮解析。

    Parameters:
        peak_sequences: FASTA ファイル (ピーク中心 ±250 bp)
        background_sequences: ランダムゲノム領域 FASTA
        jaspar_db: JASPAR データベースバージョン
    """
    from Bio import motifs

    # JASPAR PWM スキャン (概念的コード)
    results = []

    # Homer / MEME-ChIP 呼び出し (実運用)
    print(f"  Scanning {jaspar_db} motifs against peak sequences...")
    print("  (Using FIMO from MEME Suite for motif scanning)")

    # Fisher exact test による濃縮
    # peak_hits / peak_total vs bg_hits / bg_total
    # for each motif in JASPAR database

    return results


def differential_binding_analysis(sample_sheet, peaks_dir,
                                   contrast=("Treatment", "Control"),
                                   fdr_cutoff=0.05, fold_change_cutoff=2):
    """
    DiffBind による差次結合解析。

    Parameters:
        sample_sheet: DiffBind サンプルシート CSV
        contrast: (treatment, control) 比較群
        fdr_cutoff: FDR 閾値
        fold_change_cutoff: log2FC 閾値
    """
    # R/rpy2 経由で DiffBind を呼び出し
    import subprocess
    r_script = f"""
    library(DiffBind)
    samples <- read.csv("{sample_sheet}")
    dba <- dba(sampleSheet=samples)
    dba <- dba.count(dba)
    dba <- dba.contrast(dba, categories=DBA_CONDITION)
    dba <- dba.analyze(dba)
    db_sites <- dba.report(dba, th=, fold=)
    write.csv(as.data.frame(db_sites), "results/diffbind_results.csv")
    """

    ()
    ()

     r_script

References

Output Files

ファイル	形式
results/chipseq/{name}_peaks.narrowPeak	BED/narrowPeak
results/chipseq/{name}_peaks.broadPeak	BED/broadPeak
results/atacseq/fragment_size_dist.csv	CSV
results/methylation/dmr_results.csv	CSV
results/chromhmm/emissions_{n}.txt	TSV
results/hic/tad_boundaries.bed	BED
results/hic/compartments.csv	CSV
results/diffbind_results.csv	CSV
figures/chromatin_state_heatmap.png	PNG
figures/hic_contact_map.png	PNG

利用可能ツール

ToolUniverse SMCP 経由で利用可能な外部ツール。

カテゴリ	主要ツール	用途
ChIP-Atlas	ChIPAtlas_enrichment_analysis	TF/ヒストン修飾エンリッチメント解析
ChIP-Atlas	ChIPAtlas_get_experiments	実験メタデータ取得 (43 万+実験)
ChIP-Atlas	ChIPAtlas_get_peak_data	ピークコールデータ取得
ChIP-Atlas	ChIPAtlas_search_datasets	データセット検索 (抗原/細胞種)
4DN	FourDN_search_data	Hi-C/ChIA-PET 3D ゲノムデータ検索
JASPAR	jaspar_search_matrices	転写因子結合モチーフ (PWM) 検索
JASPAR	jaspar_get_matrix	PWM (Position Weight Matrix) 取得
JASPAR	jaspar_list_collections	JASPAR コレクション一覧
SCREEN	SCREEN_get_regulatory_elements	cCRE (候補シス調節エレメント) 取得
ENCODE	ENCODE_search_experiments	ENCODE ChIP-seq/ATAC-seq 実験検索
ENCODE	ENCODE_get_experiment	ENCODE 実験詳細取得
ENCODE	ENCODE_list_files	ENCODE ファイル一覧

参照スキル

スキル	関連
scientific-single-cell-genomics	scATAC-seq 連携
scientific-sequence-analysis	ゲノム配列操作
scientific-bioinformatics	BAM/VCF 処理
scientific-population-genetics	eQTL・調節バリアント
scientific-gene-expression-transcriptomics	発現-エピゲノム統合

依存パッケージ

macs3, cooler, pybedtools, deeptools, scikit-learn, scipy, pandas, numpy, biopython