染色质重塑因子分析¶

染色质重塑 = "搬运工移动路障"

DNA缠绕在核小体上 = 路被路障挡住
染色质重塑因子 = 搬运工，用ATP的能量移动路障

四种搬运方式：
  ┌─ 滑动(Sliding)：把核小体沿DNA滑动
  ├─ 弹出(Eviction)：把核小体直接移走
  ├─ 替换(Editing)：换掉核小体中的组蛋白变体
  └─ 间距调节(Spacing)：调整核小体之间的距离

四大重塑因子家族：
  SWI/SNF (BAF/PBAF)  → 打开染色质，激活基因（肿瘤抑制子）
  ISWI (SNF2H/SNF2L)  → 排列核小体，维持规则间距
  CHD (CHD1-9/NuRD)   → 双功能：激活或抑制基因
  INO80 (INO80/SWR1)  → DNA修复、组蛋白变体替换

与疾病关系：
  SWI/SNF基因突变 → ~20%人类肿瘤
  CHD突变 → 自闭症、先天性心脏病
  NuRD异常 → 白血病

# ========== MNase-seq分析流程 ==========
# MNase酶消化暴露的DNA，只保留被核小体保护的片段

# 比对（与普通ChIP-seq相同）
bowtie2 -x hg38_index -1 R1.fq.gz -2 R2.fq.gz --no-mixed -X 300 -p 8 | \
    samtools sort -o mnase.bam

# 去重
picard MarkDuplicates I=mnase.bam O=mnase_dedup.bam M=dup.txt REMOVE_DUPLICATES=true

# 只保留单核小体片段（147bp ± 20bp）
samtools view -b mnase_dedup.bam | \
    python3 -c "
import pysam
bam_in = pysam.AlignmentFile('-', 'rb')
bam_out = pysam.AlignmentFile('mnase_mono.bam', 'wb', template=bam_in)
for read in bam_in:
    if 127 <= abs(read.template_length) <= 167:  # 147±20bp
        bam_out.write(read)
bam_out.close()
"

# 生成核小体定位信号
bamCoverage -b mnase_mono.bam -o nucleosome.bw --binSize 1 --normalizeUsing RPKM -p 8

#!/usr/bin/env python3
"""核小体定位分析"""

import pyBigWig  # 读取BigWig文件
import numpy as np  # 数值计算
import matplotlib.pyplot as plt  # 画图

# ========== TSS处核小体定位图 ==========
def nucleosome_profile_at_tss(bigwig_file, tss_file, upstream=2000, downstream=2000):
    """计算TSS附近的核小体定位信号"""

    bw = pyBigWig.open(bigwig_file)  # 打开BigWig文件

    tss_list = []
    with open(tss_file) as f:
        for line in f:
            parts = line.strip().split('\t')
            tss_list.append({
                "chr": parts[0],
                "pos": int(parts[1]),
                "strand": parts[5] if len(parts) > 5 else "+"
            })

    # 计算平均信号
    window_size = upstream + downstream
    signal_sum = np.zeros(window_size)
    count = 0

    for tss in tss_list[:5000]:  # 取前5000个基因
        chrom = tss["chr"]
        pos = tss["pos"]
        start = pos - upstream
        end = pos + downstream

        if start < 0:
            continue

        try:
            values = bw.values(chrom, start, end)  # 获取信号值
            values = np.array(values)
            values = np.nan_to_num(values)  # NaN替换为0

            if tss.get("strand") == "-":  # 负链反转
                values = values[::-1]

            signal_sum += values
            count += 1
        except:
            continue

    signal_avg = signal_sum / count  # 平均信号
    bw.close()

    # 画图
    positions = np.arange(-upstream, downstream)
    plt.figure(figsize=(12, 6))
    plt.plot(positions, signal_avg, color="steelblue", linewidth=1.5)
    plt.axvline(0, color="red", linestyle="--", alpha=0.5, label="TSS")
    plt.xlabel("Distance from TSS (bp)", fontsize=14)
    plt.ylabel("Average MNase-seq Signal", fontsize=14)
    plt.title("Nucleosome Positioning at TSS", fontsize=16)

    # 标注核小体位置
    plt.text(-150, max(signal_avg)*0.9, "-1", fontsize=12, ha="center", color="gray")
    plt.text(150, max(signal_avg)*0.9, "+1", fontsize=12, ha="center", color="gray")
    plt.text(0, max(signal_avg)*0.1, "NFR", fontsize=12, ha="center", color="red")

    plt.legend()
    plt.tight_layout()
    plt.savefig("nucleosome_tss.png", dpi=300)
    plt.close()

    print(f"分析了 {count} 个基因的TSS区域")
    return signal_avg

# 使用
# nucleosome_profile_at_tss("nucleosome.bw", "hg38_tss.bed")

#!/usr/bin/env Rscript
# 分析重塑因子结合位点与核小体变化的关系

library(GenomicRanges)
library(rtracklayer)

# ========== 读取数据 ==========
# 重塑因子（如BRG1/SMARCA4）ChIP-seq peaks
brg1_peaks <- import("BRG1_peaks.narrowPeak", format = "BED")

# ATAC-seq peaks（开放区域）
atac_peaks <- import("ATAC_peaks.narrowPeak", format = "BED")

# ========== 分析重塑因子与开放区域的重叠 ==========
overlap <- subsetByOverlaps(brg1_peaks, atac_peaks)
cat("BRG1与ATAC-seq重叠区域:", length(overlap), "\n")
cat("重叠比例:", length(overlap) / length(brg1_peaks) * 100, "%\n")

# ========== 差异分析：KO vs WT ==========
# 比较重塑因子敲除前后的ATAC-seq变化
wt_peaks <- import("WT_ATAC_peaks.narrowPeak", format = "BED")
ko_peaks <- import("KO_ATAC_peaks.narrowPeak", format = "BED")

# WT特有（KO中丢失）= 重塑因子依赖的开放区域
wt_only <- subsetByOverlaps(wt_peaks, ko_peaks, invert = TRUE)
cat("重塑因子依赖的开放区域:", length(wt_only), "\n")

# KO特有（WT中没有）= 重塑因子抑制的区域
ko_only <- subsetByOverlaps(ko_peaks, wt_peaks, invert = TRUE)
cat("重塑因子抑制的区域:", length(ko_only), "\n")

========================================
染色质重塑因子分析 速查表
========================================

【四大重塑因子家族】
SWI/SNF              → BAF/PBAF，打开染色质（~20%肿瘤突变）
ISWI                 → 核小体间距调节
CHD                  → 染色质域，双功能
INO80                → DNA修复，组蛋白变体替换

【核小体定位分析】
MNase-seq            → 经典方法，检测核小体保护的DNA
ATAC-seq             → 检测开放区域（间接看核小体缺失）
NOMe-seq             → 同时检测甲基化和核小体

【关键概念】
NFR                  → 无核小体区域（TSS处）
+1核小体            → TSS下游第一个核小体
-1核小体            → TSS上游第一个核小体
核小体间距           → 约200bp（含147bp DNA + linker）

【SWI/SNF肿瘤相关基因】
SMARCA4(BRG1)        → 催化亚基，多种肿瘤突变
SMARCB1(SNF5)        → 核心亚基，横纹肌样瘤
ARID1A               → BAF特异亚基，子宫内膜/卵巢癌
PBRM1                → PBAF特异亚基，肾癌

【面试考点】
Q: SWI/SNF和ISWI的区别？
A: SWI/SNF主要弹出/滑动核小体打开染色质，ISWI主要调节核小体间距

Q: 为什么SWI/SNF突变与肿瘤有关？
A: SWI/SNF打开肿瘤抑制基因的染色质，突变后这些基因被沉默
========================================