怎么使用SICER進(jìn)行peak calling

怎么使用SICER進(jìn)行peak calling，針對(duì)這個(gè)問題，這篇文章詳細(xì)介紹了相對(duì)應(yīng)的分析和解答，希望可以幫助更多想解決這個(gè)問題的小伙伴找到更簡(jiǎn)單易行的方法。

chip_seq數(shù)據(jù)中peak的長(zhǎng)度范圍跨度較大，既有覆蓋幾個(gè)核小體的幾百bp的peak, 也有包含多個(gè)基因長(zhǎng)度在上千kb的peak。比如H3K4me2和H3K4me3這兩種組蛋白修飾中peak在幾百bp左右, 而H3K27me3中則為長(zhǎng)度在幾十到幾百kb之間。組蛋白修飾中peak長(zhǎng)度跨度大，弱信號(hào)分散都特點(diǎn)，使得基于轉(zhuǎn)錄因子TF結(jié)合位點(diǎn)的peak calling軟件在分析這類數(shù)據(jù)時(shí)準(zhǔn)確度較差。

SICER是一款專門針對(duì)組蛋白修飾的chip數(shù)據(jù)進(jìn)行peak calling的軟件，核心思想也是基于滑動(dòng)窗口和局部泊松分布的方式來識(shí)別富集區(qū)域，下圖所示為該軟件用默認(rèn)參數(shù)識(shí)別到的H3K27me3的peak區(qū)域

黑色區(qū)域?yàn)镋NCODE分析得到的peak區(qū)域，紅色區(qū)域?yàn)镾ICER分析得到的peak區(qū)域。該軟件官網(wǎng)如下

https://home.gwu.edu/~wpeng/Software.htm

為例方便使用，有人對(duì)該軟件進(jìn)行了分裝，使用起來更加方便，源代碼托管在github上，網(wǎng)址如下

https://github.com/dariober/SICERpy

基本用法如下

python SICERpy \
-c input.bam \
-w 200 \
-g 3 \
-t ip.bam \
> peak.bed

-w參數(shù)表示滑動(dòng)窗口的大小，默認(rèn)值為200。數(shù)值越小， 識(shí)別到的peak區(qū)間長(zhǎng)度相對(duì)越短且越分散；數(shù)值越大，會(huì)造成過渡擬合，識(shí)別到的peak區(qū)間過長(zhǎng)，丟失掉真實(shí)的信息，示意如下

對(duì)于轉(zhuǎn)錄因子，官方推薦滑動(dòng)窗口設(shè)置為50-100bp， 對(duì)于組蛋白修飾，推薦設(shè)置為200bp。

-g參數(shù)代表gap的大小，默認(rèn)值為3。和windows size類似，該參數(shù)也直接影響peak區(qū)間的定義，示意如下

對(duì)于轉(zhuǎn)錄因子，官方推薦該數(shù)值和滑動(dòng)窗口數(shù)值保持相同；對(duì)于組蛋白修飾，推薦值為3。

輸出文件為bed格式，共8列，每列含義如下

1. chrom

2. start

3. end

4. chip read count

5. input read count

6. pvalue

7. fold_enrichment

8. fdr

   
   

可以最后一列的fdr值，來篩選得到高可信度的peak信息，用法如下

awk '$8 < 0.01' peaks.bed > peaks.01.bed

關(guān)于怎么使用SICER進(jìn)行peak calling問題的解答就分享到這里了，希望以上內(nèi)容可以對(duì)大家有一定的幫助，如果你還有很多疑惑沒有解開，可以關(guān)注億速云行業(yè)資訊頻道了解更多相關(guān)知識(shí)。