怎么使用SICER进行peak calling

发布时间：2021-11-10 10:17:03 作者：柒染
来源：亿速云阅读：198

# 怎么使用SICER进行peak calling

## 一、SICER简介

SICER（Spatial Clustering for Identification of ChIP-Enriched Regions）是一种专门用于分析ChIP-seq数据的peak calling算法，由Zang等人于2009年开发。与常规peak caller（如MACS2）不同，SICER通过结合**空间聚类**和**统计建模**，特别适合识别**弥散型修饰**（如H3K27me3、H3K9me2等）的广泛富集区域。

### 核心优势
- 能识别宽peak和弱信号区域
- 通过窗口聚类减少假阳性
- 支持重复样本的联合分析

---

## 二、安装与环境配置

### 1. 系统要求
- Linux/Unix系统（推荐Ubuntu/CentOS）
- Python 2.7（SICER原生版本依赖）
- R语言（用于统计分析和可视化）

### 2. 安装步骤
```bash
# 从官网下载SICER
wget http://home.gwu.edu/~wpeng/Software/SICER_V1.1.tar.gz
tar -xzvf SICER_V1.1.tar.gz
cd SICER_V1.1

# 安装依赖
pip install numpy scipy

注意：若使用Python 3，需手动修改SICER和SICER-rb.sh中的语法（如print语句）。

三、数据预处理

1. 输入文件格式

实验组（Treatment）和对照组（Control）的BAM/BED文件
需转换为SICER要求的BED格式（chr、start、end三列）

# 使用bedtools转换BAM到BED
bedtools bamtobed -i treatment.bam > treatment.bed

2. 去除冗余reads

# 使用SICER自带的去重脚本
python SICER-df.sh treatment.bed output_dir hg38 1 200 0.74 600 0.01

参数说明： - 1：物种对应的基因组版本ID（1=hg38） - 200：读段延伸长度（需与实验片段大小一致） - 0.74：冗余读段阈值

四、运行Peak Calling

基础命令

python SICER.sh /path/to/bed_files treatment.bed control.bed output_dir hg38 1 200 150 0.74 600 0.01

关键参数解析

参数	说明
200	窗口大小（bp）
150	间隙允许距离（bp）
0.01	FDR阈值
600	有效基因组长度比例

输出文件

treatment-W200-G600-islands-summary：peak统计信息
treatment-W200-G600-graph：可视化数据
treatment-W200-G600-islands.bed：最终peak文件

五、结果解读与优化

1. 评估peak质量

FRiP（Fraction of Reads in Peaks）：计算落在peak区域的reads比例
```
bedtools intersect -a treatment.bed -b output_dir/*-islands.bed | wc -l
```

Peak宽度分布：使用R绘制直方图

peaks <- read.table("treatment-W200-G600-islands.bed")
hist(peaks$V3-peaks$V2, main="Peak Width Distribution")

2. 参数调优建议

宽peak数据：增大窗口（如500bp）和间隙（如300bp）
高噪声数据：提高FDR阈值（如0.05）
小样本量：降低有效基因组比例（如400）

六、高级应用

1. 差异Peak分析

使用SICER-df.sh比较两组ChIP-seq数据：

python SICER-df.sh group1.bed group2.bed output_dir hg38 1 200 0.74 600 0.01

2. 与MACS2联用

先用MACS2检测sharp peak，再用SICER检测宽peak：

macs2 callpeak -t treatment.bam -c control.bam -f BAM -g hs -n sharp_peaks
python SICER.sh treatment.bed control.bed output_dir hg38 1 500 300 0.74 600 0.05

七、常见问题解决

1. 报错”Invalid BED file”

检查BED文件是否含header
确保染色体格式一致（chr1 vs 1）

2. 运行速度慢

预筛选高表达区域：bedtools genomecov -bg -i treatment.bed | awk '$4>10' > filtered.bed
使用split命令分割大文件并行处理

3. 结果空文件

检查FDR是否过于严格
确认对照组是否覆盖度不足

八、参考文献

Zang C, et al. (2009). Bioinformatics 25(15):1952-8
SICER官方文档：http://home.gwu.edu/~wpeng/Software.htm

作者注：本文基于SICER 1.1版本撰写，实际使用时请参考最新版本更新日志。 “`

这篇文章共计约1150字，包含了安装指南、参数详解、结果分析和故障排查等完整流程，采用Markdown格式便于阅读和编辑。需要调整细节可进一步补充具体案例或截图说明。