2022-09-13 calling peaks问题

1. 利用macs2 callpeak进行 calling peaks,其中有个参数--keep-dup,这是如果输入为过滤掉重复reads的bam文件,这--keep-dup all。如果为未过滤重复reads的bam文件,--keep-dup 1。此参数是指定输入bam文件中重复reads有无的。

macs2 callpeak 的基本命令

-t

IP.bam

-c

input.bam

-g

genome.size

-B

输出bgd文件,下游bigwig文件生成所需

-f

双端测序使用BAMPE,单端的话不需要加参数(或 -f BAM),默认是auto识别。除”BAMPE”, “BEDPE”需要特别声明外,其他格式都可以用 AUTO自动检测。

标签文件的格式,可以是“ELAND”,“BED”,“ELANDMULTI”,“ELANDEXPORT”,“ELANDMULTIPET”(用于对端标签),“SAM”,“BAM”,“BOWTIE”,“BAMPE”或“BEDPE”。默认为“AUTO”,这将允许MACS自动决定格式。当您使用“AUTO”时也会使用

结合不同格式的文件。请注意,MACS无法检测到“BAMPE”或“BEDPE”格式带有“AUTO”,你必须隐含指定“BAMPE”和“BEDPE”的格式。

格式指定'BAMPE'或'BEDPE'时将触发特殊模式。这样,MACS2将处理BAM或BED文件作为配对结束数据。而不是建立双峰分布正负链读数预测片段大小,MACS2会使用读取对的实际插入大小来构建片段积累。(所以,当你的数据是双端测序数据时,你应该用BAMPE或者BEDPE参数。当你设置成双端参数的时候,MACS2就会跳过建模计算d的那一步,而是直接用片段的insert size来建立堆积。)

BAMPE格式只是包含配对末端对齐的BAM格式信息,例如来自BWA或BOWTIE的信息。

BEDPE格式是一种简化且更灵活的BED格式只包含定义染色体名称的前三列,来自Paired-end的片段的左右位置测序。请注意,这与BEDTOOLS使用的格式不同,BEDTO的BEDTOOLS版本实际上不在标准BED中格式。

-q

设置FDR阈值

-p

设置pvalue阈值

--nomodel

MACS 不构建模型。这个参数和extsize、shift是配套使用的,有这个参数才可以设置extsize和shift。

--extsize

当设置了nomodel时,MACS会用--extsize这个参数从5'->3'方向扩展reads修复fragments。比如说你的转录因子结合范围200bp,就设置这个参数是200。

--shift

当设置了--nomodel,MACS用这个参数从5' 端移动剪切,然后用--extsize延伸,如果--shift是负值表示从3'端方向移动。建议ChIP-seq数据集这个值保持默认值为0,对于检测富集剪切位点如DNAsel数据集设置为EXTSIZE的一半。

--SPMR

需要-B被设置,不影响FDR和pvalue

--outdir

输出文件的路径

--broad

peak有narrow peak和broad peak, 设置时可以call broad peak 的结果文件。

--broad-cutoff

和pvalue、以及qvalue相似

其实,这里面讨论最多的是--nomodel --shift -100 --extsize 200这些参数如何选择,下面的图很形象的展示了参数的作用。当然,我也是查阅了很多资料与文献,

一般默认在ATAC-seq,DNase-seq,FAIRE-seq的时候将shift设置为extsize的一半,且参数固定为:--nomodel --shift -100 --extsize 200 (猪项目中为shift -75 --extsize 150)。

而在MNase-seq的时候,参数固定为:--nomodel --shift 37 --extsize 73

在ChiP-seq的时候不用移峰,所以只使用-nomodel!!当做组蛋白修饰的时候,由于peak并不典型,所以使用--nomodel --extsize 73参数。


对人细胞系ATAC-seq 数据call peak的参数设置如下:

macs2 callpeak -t sample.final.bam -n sample --nomodel --shift -100 --extsize 200 -B --SPMR -g hs --outdir Macs2_out 2 --keep-dup all --call-summits > sample.macs2.log (单端read时,bam已经去过重复则--keep-dup all,不去的话--keep-dup设置为1)

macs2 callpeak -t sample.final.bam -n sample -f BAMPE -B --SPMR -g hs --outdir Macs2_out 2 --keep-dup all --call-summits > sample.macs2.log(双端reads)


思考:-f 设置为BAMPE时,好像与--shift -100 --extsize 200 --nomode冲突,因为加不加这三个参数,结果是一样的。看-f里面,BAMPE的解释好像也是这个意思。

所以对于双末端数据可以只设置-f BAMPE与--nomode --shift -100 --extsize 200 参数二选一???



参考1:组蛋白ChIP分析要注意的2个要点(基迪奥):https://www.genedenovo.com/news/333.html

–shiftsize已经被 –extsize所替代;(https://www.imooc.com/article/270403)


参考2: https://www.jianshu.com/p/e83a7e10ea2e?tdsourcetag=s_pcqq_aiomsg

参考3: https://www.jianshu.com/p/9aa719faa4b5


-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

---------------------------------------------------------最近macs3说明书建议分析方法---------------------------------------------------------------------------------------

Examples:

1. Peak calling for regular TF ChIP-seq:

    $ macs3 callpeak -t ChIP.bam -c Control.bam -f BAM  -g hs -n test -B -q 0.01

2. Broad peak calling on Histone Mark ChIP-seq:

    $ macs3 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1

3. Peak calling on ATAC-seq (paired-end mode):

    $ macs3 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01

4. Peak calling on ATAC-seq ( focusing on insertion sites, and using single-end mode):

    $ macs3 callpeak -f BAM -t ATAC.bam -g hs -n test -B -q 0.01 --shift -50 --extension 100


-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

---------------------------------------------------------------综上感觉对的应该是-------------------------------------------------------------------------------------------------------

1.分析普通转录因子时:

macs3 callpeak -t ChIP.bam -c Control.bam -f BAMPE -g hs -n sample -B -q 0.01 (双末端测序,--keep-dup all必须根据实际情况申明)

2.分析组蛋白ChiP时

鉴定方法1:(建议优选这个吧???

Narrow peak: macs3 callpeak -t ChIP.bam -c Control.bam -n sample -f BAMPE -B -g hs -q 0.05 (双末端测序,--keep-dup all必须根据实际情况申明)

Broad peak:macs3 callpeak -t ChIP.bam -c Control.bam -n sample -f BAMPE -B -g hs -q 0.05 --broad --broad-cutoff 0.1 (双末端测序,--keep-dup all必须根据实际情况申明)

注:DNA methylation underpins the epigenomic landscape regulating genome transcription in Arabidopsis (GB论文),就用了这个方法鉴定peaks,但是-q和--broad-cutoff参数(-q 0.00001 --broad-cutoff 0.00001)更为严格,可以根据实际需求进行调整。

鉴定方法2(如基迪奥推荐方法,用--nomodel,--extsize 指定拓展值)

Narrow peak:macs2 callpeak -t ChIP.bam -c Control.bam -n sample -g hs --keep-dup all --nomodel --shift 0 --extsize 200  -q 0.01

Broad peak:macs2 callpeak -t ChIP.bam -c Control.bam -n sample -g hs --keep-dup all --nomodel --shift 0 --extsize 200 -q 0.01 --broad --broad-cutoff 0.1

注:A chromatin integration labelling method enables epigenomic profiling with lower input (nature cell biology),就用了这个方法鉴定peaks,但是基迪奥推荐方法推荐--extsize 73,《猪基因组顺势元件鉴定及功能SNP注释》,建议--extsize n,n由计算获得。

3. 分析ATAC-seq时

鉴定方法1:(建议优选这个

macs2 callpeak -t ATAC.bam -n sample --nomodel --shift -100 --extsize 200 -B --SPMR -g hs --outdir Macs2_out 2 --keep-dup all --call-summits

注:但是--shift -100 --extsize 200这个所用论文极多,但是也有其他数值,如《猪基因组顺势元件鉴定及功能SNP注释》用的是shift -75 --extsize 150;macs3又推荐shift -50 --extsize 100;最新的Genomic innovation and regulatory rewiring during evolution of the cotton genus Gossypium (NG棉花论文)用的是—extsize 38 - shift -15(UMI-ATAC-seq),总体来说五花八门。

鉴定方法2:

macs2 callpeak -f BAMPE -t ATAC.bam -g hs -n sample -B -q 0.01


20240602更新:

在基因课最新的课程中

TF和组蛋白修饰:

进行双末端测序时,直接-f BAMPE;单末端测序,直接-f BAM,如果不能有效估计d值,在手动设置--nomodel --shift --extsize等参数

在ATAC-seq时:

建议-f BAM --nomodel --shift -100 --extsize 200

基本和我的建议一样。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,172评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,346评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 159,788评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,299评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,409评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,467评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,476评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,262评论 0 269
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,699评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,994评论 2 328
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,167评论 1 343
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,827评论 4 337
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,499评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,149评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,387评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,028评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,055评论 2 352

推荐阅读更多精彩内容