怎么使用SICER进行peak calling,针对这个问题,这篇文章详细介绍了相对应的分析和解答,希望可以帮助更多想解决这个问题的小伙伴找到更简单易行的方法。
chip_seq数据中peak的长度范围跨度较大,既有覆盖几个核小体的几百bp的peak, 也有包含多个基因长度在上千kb的peak。比如H3K4me2和H3K4me3这两种组蛋白修饰中peak在几百bp左右, 而H3K27me3中则为长度在几十到几百kb之间。组蛋白修饰中peak长度跨度大,弱信号分散都特点,使得基于转录因子TF结合位点的peak calling软件在分析这类数据时准确度较差。
SICER是一款专门针对组蛋白修饰的chip数据进行peak calling的软件,核心思想也是基于滑动窗口和局部泊松分布的方式来识别富集区域,下图所示为该软件用默认参数识别到的H3K27me3的peak区域
黑色区域为ENCODE分析得到的peak区域,红色区域为SICER分析得到的peak区域。该软件官网如下
https://home.gwu.edu/~wpeng/Software.htm
为例方便使用,有人对该软件进行了分装,使用起来更加方便,源代码托管在github上,网址如下
https://github.com/dariober/SICERpy
基本用法如下
python SICERpy \
-c input.bam \
-w 200 \
-g 3 \
-t ip.bam \
> peak.bed
-w
参数表示滑动窗口的大小,默认值为200。数值越小, 识别到的peak区间长度相对越短且越分散;数值越大,会造成过渡拟合,识别到的peak区间过长,丢失掉真实的信息,示意如下
对于转录因子,官方推荐滑动窗口设置为50-100bp, 对于组蛋白修饰,推荐设置为200bp。
-g
参数代表gap的大小,默认值为3。和windows size类似,该参数也直接影响peak区间的定义,示意如下
对于转录因子,官方推荐该数值和滑动窗口数值保持相同;对于组蛋白修饰,推荐值为3。
输出文件为bed格式,共8列,每列含义如下
chrom
start
end
chip read count
input read count
pvalue
fold_enrichment
fdr
可以最后一列的fdr值,来筛选得到高可信度的peak信息,用法如下
awk '$8 < 0.01' peaks.bed > peaks.01.bed
关于怎么使用SICER进行peak calling问题的解答就分享到这里了,希望以上内容可以对大家有一定的帮助,如果你还有很多疑惑没有解开,可以关注亿速云行业资讯频道了解更多相关知识。
免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。