转录因子结合位点预测¶

# Python: 使用JASPAR数据库获取TF motif
from Bio import motifs
from Bio.motifs.jaspar.db import JASPAR5
import numpy as np

# 或者直接从JASPAR网站下载PFM文件
# https://jaspar.genereg.net/

# 读取JASPAR格式motif
with open("MA0106.1.jaspar") as f:
    m = motifs.read(f, "jaspar")

print(m.name)        # ETS1
print(m.matrix_id)   # MA0106.1
print(m.counts)      # PFM矩阵
print(m.pwm)         # PWM矩阵
print(m.pssm)        # PSSM (log-odds)

# 计算信息含量
print(m.pssm.mean())
print(m.pssm.std())

# 扫描序列
from Bio.Seq import Seq
test_seq = Seq("ATCGGAAGTCGATCGATCG")
for pos, score in m.pssm.search(test_seq, threshold=3.0):
    print(f"Position {pos}: score {score:.2f}")

# 安装MEME Suite
conda install -c bioconda meme

# 从ChIP-seq peak提取序列
# 先取peak summit ±150bp
awk '{OFS="\t"; mid=int(($2+$3)/2); print $1,mid-150,mid+150}' \
    peaks.narrowPeak > peak_regions.bed

bedtools getfasta -fi genome.fa -bed peak_regions.bed \
    -fo peak_sequences.fasta

# 方法1: MEME (经典，适合<1000条序列)
meme peak_sequences.fasta \
    -dna \
    -oc meme_output \
    -mod zoops \          # 每条序列0或1个motif实例
    -nmotifs 10 \         # 找10个motif
    -minw 6 -maxw 20 \   # motif长度范围
    -revcomp \            # 考虑反向互补
    -markov_order 1       # 1阶马尔可夫背景模型

# 方法2: STREME (快速，适合大数据集)
streme --dna \
    --oc streme_output \
    -p peak_sequences.fasta \
    --minw 6 --maxw 20 \
    --thresh 0.05

# 方法3: HOMER
findMotifsGenome.pl peaks.narrowPeak hg38 homer_output/ \
    -size 200 \
    -mask \
    -p 8

# HOMER优势：自动与已知motif比较，报告p-value和百分比

# 从JASPAR下载motif文件（MEME格式）
# https://jaspar.genereg.net/downloads/

# 单个motif扫描
fimo --oc fimo_output \
    --thresh 1e-4 \
    --bfile background_model.txt \
    MA0139.1.meme \      # CTCF motif
    peak_sequences.fasta

# 全基因组扫描（使用所有JASPAR motif）
fimo --oc fimo_genome \
    --thresh 1e-5 \
    --max-stored-scores 10000000 \
    JASPAR2024_CORE_vertebrates.meme \
    genome.fa

# 只扫描开放染色质区域（大幅减少假阳性）
bedtools getfasta -fi genome.fa \
    -bed atac_peaks.bed -fo open_chromatin.fasta

fimo --oc fimo_open_regions \
    --thresh 1e-4 \
    JASPAR2024_CORE_vertebrates.meme \
    open_chromatin.fasta

# 解析FIMO输出
# fimo.tsv包含: motif_id, sequence_name, start, stop, strand, score, p-value, q-value
head fimo_output/fimo.tsv

# 过滤结果
awk -F'\t' '$8 < 0.05' fimo_output/fimo.tsv > significant_sites.tsv

# AME (Analysis of Motif Enrichment)
ame --oc ame_output \
    --control --shuffle-- \   # 对照：打乱输入序列
    --method fisher \
    peak_sequences.fasta \
    JASPAR2024_CORE_vertebrates.meme

# 使用特定背景序列
ame --oc ame_output \
    --control background_sequences.fasta \
    --method ranksum \
    peak_sequences.fasta \
    JASPAR2024_CORE_vertebrates.meme

# CentriMo: 检测motif在peak中心的富集
centrimo --oc centrimo_output \
    --ethresh 10 \
    peak_sequences.fasta \
    JASPAR2024_CORE_vertebrates.meme

# 如果motif在peak中心显著富集，说明该TF可能是direct binder

# 方法1: TOBIAS (推荐，专门为ATAC-seq设计)
# 安装
pip install tobias

# 步骤1：Tn5偏好校正
TOBIAS ATACorrect \
    --bam atac_sorted.bam \
    --genome genome.fa \
    --peaks atac_peaks.bed \
    --outdir tobias_output/ \
    --cores 8

# 步骤2：计算footprint scores
TOBIAS FootprintScores \
    --signal tobias_output/*_corrected.bw \
    --regions atac_peaks.bed \
    --output tobias_output/footprint_scores.bw \
    --cores 8

# 步骤3：识别bound/unbound位点
TOBIAS BINDetect \
    --motifs JASPAR2024_CORE_vertebrates.meme \
    --signals tobias_output/condition1_footprints.bw \
              tobias_output/condition2_footprints.bw \
    --genome genome.fa \
    --peaks atac_peaks.bed \
    --outdir tobias_bindetect/ \
    --cores 8 \
    --cond_names condition1 condition2

# 步骤4：可视化聚合footprint
TOBIAS PlotAggregate \
    --TFBS tobias_bindetect/CTCF/beds/CTCF_condition1_bound.bed \
    --signals tobias_output/condition1_corrected.bw \
              tobias_output/condition2_corrected.bw \
    --output CTCF_footprint.pdf \
    --share_y both \
    --plot_boundaries

# 方法2: HINT-ATAC
rgt-hint footprinting \
    --atac-seq \
    --paired-end \
    --organism=hg38 \
    --output-location=hint_output/ \
    atac_sorted.bam atac_peaks.bed

# R: 整合分析
library(motifmatchr)
library(chromVAR)
library(GenomicRanges)
library(BSgenome.Hsapiens.UCSC.hg38)

# 使用chromVAR评估TF活性变异
# 读取peak计数矩阵
counts <- readRDS("peak_counts_matrix.rds")
peaks <- readRDS("peak_granges.rds")

# 获取JASPAR motif
library(JASPAR2024)
library(TFBSTools)
motifs <- getMatrixSet(JASPAR2024,
    opts=list(collection="CORE", tax_group="vertebrates"))

# 匹配motif到peaks
motif_ix <- matchMotifs(motifs, peaks,
    genome = BSgenome.Hsapiens.UCSC.hg38,
    out = "scores")

# chromVAR分析
library(chromVAR)
dev <- computeDeviations(object = counts, 
                         annotations = motif_ix)

# 提取TF活性z-scores
tf_scores <- deviationScores(dev)

# 差异TF活性
# 使用limma或t-test比较条件间TF活性差异
library(limma)
design <- model.matrix(~condition)
fit <- lmFit(tf_scores, design)
fit <- eBayes(fit)
results <- topTable(fit, number=Inf)

# === 完整TF结合位点分析流程 ===

# 1. 获取motif数据库
wget https://jaspar.genereg.net/download/data/2024/CORE/JASPAR2024_CORE_vertebrates_non-redundant_pfms_meme.txt

# 2. 从差异ATAC peaks提取序列
bedtools getfasta -fi hg38.fa -bed diff_peaks.bed -fo diff_peaks.fasta

# 3. De novo motif发现
meme diff_peaks.fasta -dna -oc meme_out -mod zoops -nmotifs 15 -minw 6 -maxw 20 -revcomp

# 4. 已知motif富集
ame --oc ame_out --control --shuffle-- diff_peaks.fasta JASPAR2024_CORE_vertebrates.meme

# 5. Motif扫描
fimo --oc fimo_out --thresh 1e-4 JASPAR2024_CORE_vertebrates.meme diff_peaks.fasta

# 6. TF Footprinting
TOBIAS ATACorrect --bam atac.bam --genome hg38.fa --peaks peaks.bed --outdir tobias/ --cores 16
TOBIAS FootprintScores --signal tobias/*_corrected.bw --regions peaks.bed --output tobias/scores.bw --cores 16
TOBIAS BINDetect --motifs JASPAR.meme --signals tobias/scores.bw --genome hg38.fa --peaks peaks.bed --outdir bindetect/ --cores 16

Level 1: 序列偏好（PWM匹配）→ 必要但不充分
Level 2: 染色质可及性（开放区域）→ 过滤大部分假阳性
Level 3: 辅因子共结合 → 确定细胞类型特异性
Level 4: 3D基因组结构 → 确定作用靶基因
Level 5: 表观遗传修饰 → H3K4me1增强子标记