bdg2bw(macs2的peak要修正坐标）
Run MACS2 combining replicates.
You don’t need to use “samtools merge”, just

$ macs2 callpeak -t H3K36me1_EE_rep1.bam H3K36me1_EE_rep2.bam -c Input_EE_rep1.bam Input_EE_rep2.bam -B --nomodel --extsize 147 --SPMR -g ce -n H3K36me1_EE
https://github.com/macs3-project/MACS/wiki/Build-Signal-Track#Fix_the_bedGraph_and_convert_them_to_bigWig_files
bigwig归一化方式详解
https://blog.csdn.net/weixin_43569478/article/details/108079478?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161495063116780255276217%2522%252C%2522scm%2522%253A%252220140713.130102334…%2522%257D&request_id=161495063116780255276217&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2_allsobaiduend~default-2-108079478.first_rank_v2_pc_rank_v29&utm_term=bamCoverage

for j in DM GM
do
ls jMyoD∗.bam∣sed′s/MyoDChIPseqRep[0−9]trimmed.bam//g′∣sort−u∣whilereadid;dosamtoolsmerge./mergeBam/{j}_MyoD*.bam|sed 's/_MyoD_ChIPseq_Rep[0-9]_trimmed.bam//g' |sort -u |while read id;do samtools merge ./mergeBam/jMyoD∗.bam∣sed′s/MyoDChIPseqRep[0−9]trimmed.bam//g′∣sort−u∣whilereadid;dosamtoolsmerge./mergeBam/id.merge.bam $id.bam ;done

done

2.4.3 bdg file → wig file
https://blog.csdn.net/weixin_41745858/article/details/82628775

为了方便在IGV上查看ChIP-seq的结果和后期的可视化展示，需要把macs2的结果转化为bw提供给IGV。一共分为三步：

第一步：使用 bdgcmp 得到 FE 或者 logLR 转化后的文件 (Run MACS2 bdgcmp to generate fold-enrichment and logLR track)
macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_FE.bdg -m FE
macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_logLR.bdg -m logLR -p 0.00001

参数解释

-m FE 计算富集倍数降低噪音

-p 为了避免log的时候input值为0时发生error，给予一个很小的值

第二步：预处理文件，下载对应参考基因组染色体长度文件

使用conda安装以下三个处理软件：
bedGraphToBigWig
bedClip
bedtools

下载染色体长度文件：http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz 并解压（针对human，其余物种的皆可以按照类似网址下载）

第三步：生成 bw 文件

#!/bin/bash

#$ -cwd
#$ -S /bin/bash
#$ -l p=10

export PATH=$PATH:~/miniconda3/envs/rna/bin
source activate rna

cd /data/zhangyong/yyp/yypold/chips/lh11_3/data/X101SC20070420-Z01-J003/allcleandata/sam/chenalign/16HD

macs2 callpeak -t DM24h-16HD-1_FKDL202607369.sort.bam -c DM24h-16HD-Input_FKDL202607368.sort.bam -B --nomodel --extsize 147 --SPMR -g mm -n H3K36me1_EE_rep1

macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_FE.bdg -m FE
macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_logLR.bdg -m logLR -p 0.00001

#######在macs2中 bdg2bw######
Fix the bedGraph and convert them to bigWig files.

You need UCSC tools from (if linux):

http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bedGraphToBigWig

http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bedClip

And bedTools:

http://code.google.com/p/bedtools/

Since sometimes, chromosome names in alignment and bedGraph files are not consistent with UCSC naming convention, you may need to fix the chromosomenames using awk or perl.

Then you need to convert bedGraph to bigWig using this simple script:https://gist.github.com/2469050

You have to download chromosome length file from UCSC too.

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz for human ( I renamed it as hg19.len )

http://hgdownload.cse.ucsc.edu/goldenPath/dm3/database/chromInfo.txt.gz for fly ( I renamed it as dm3.len )
http://hgdownload.cse.ucsc.edu/goldenPath/mm10/database/chromInfo.txt.gz
Then:

$ bdg2bw H3K36me1_EE_rep1_FE.bdg hg19.len $ bdg2bw H3K36me1_EE_rep1_logLR.bdg hg19.len
And you will have these bigwig files:

H3K36me1_EE_rep1_FE.bw H3K36me1_EE_rep1_logLR.bw
Calculate correlation between replicates.
You need UCSC tool:

http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigCorrelate

Run it:

$ wigCorrelate H3K36me1_EE_rep1_FE.bw H3K36me1_EE_rep2_FE.bw
If the correlation is good enough, you can combine replicates and run MACS2 again.

Run MACS2 combining replicates.
You don’t need to use “samtools merge”, just

$ macs2 callpeak -t H3K36me1_EE_rep1.bam H3K36me1_EE_rep2.bam -c Input_EE_rep1.bam Input_EE_rep2.bam -B --nomodel --extsize 147 --SPMR -g ce -n H3K36me1_EE
Then use the similar approaches as step 2 and 3 to generate signal tracks in bigWig format.

You can easily write a script to build a pipeline from step1 to step5.

MACS2学习笔记

为了方便在IGV上查看ChIP-seq的结果和后期的可视化展示，所以我们需要把macs2的结果转化为bw提供给IGV(bdg file to wig file transformation)
一共分为三步
1.首先使用bdgcmp得到FE或者logLR转化后的文件(Run MACS2 bdgcmp to generate fold-enrichment and logLR track)

macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_FE.bdg -m FE

macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_logLR.bdg -m logLR -p 0.00001

参数解释

-m FE 计算富集倍数降低噪音

-p 为了避免log的时候input值为0时发生error，给予一个很小的值
2. 预处理文件与对应参考基因组染色体长度文件下载

使用conda安装以下三个处理软件：bedGraphToBigWig/bedClip/bedtools下载染色体长度文件：http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz 并解压（针对human，其余物种的皆可以按照类似网址下载）写一个小小的sh脚本方便一步转化name.sh：

#!/bin/bash

check commands: slopBed, bedGraphToBigWig and bedClip

which bedtools &>/dev/null || { echo “bedtools not found! Download bedTools: http://code.google.com/p/bedtools/”; exit 1; }

which bedGraphToBigWig &>/dev/null || { echo “bedGraphToBigWig not found! Download: http://hgdownload.cse.ucsc.edu/admin/exe/”; exit 1; }

which bedClip &>/dev/null || { echo “bedClip not found! Download: http://hgdownload.cse.ucsc.edu/admin/exe/”; exit 1; }

end of checking

if [ $# -lt 2 ];then

echo "Need 2 parameters! "

exit

F=$1

G=$2

bedtools slop -i ${F} -g ${G} -b 0 | bedClip stdin ${G} ${F}.clip

LC_COLLATE=C sort -k1,1 -k2,2n ${F}.clip > ${F}.sort.clip

bedGraphToBigWig ${F}.sort.clip ${G} ${F/bdg/bw}

rm -f ${F}.clip ${F}.sort.clip

chmod +x name.sh
3. 最后就是生成bw文件

./name.sh H3K36me1_EE_rep1_FE.bdg hg19.len

./name.sh H3K36me1_EE_rep1_logLR.bdg hg19.len
最后得到产物，至于的使用哪一个作为输入文件大家就根据需要来吧

H3K36me1_EE_rep1_FE.bw

H3K36me1_EE_rep1_logLR.bw

作者：面面的徐爷
链接：https://www.jianshu.com/p/ceace75a3dcf
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

bedGraphToBigWig
之前ChIP-seq里面讲过deeptools来转格式bam2bw，现在学的是bdg2bw，是一回事吗？
这一步使用到的工具：bedSort bedClip bedGraphToBigWig。下载：http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bedSort
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bedClip
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bedGraphToBigWig
下载之后，赋予权限，再转移到一个PATH包含的目录下，
mv bedSort bedClip bedGraphToBigWig ~/local/bin/

#求染色体序列大小
samtools view -H SRR5874657.final.bam | perl -ne ‘if(/SN:(\S+)\s+LN:(\d+)/){print “$1\t$2\n”}’ > chr.size

for i in SRR58746{57…62}
do
#对treat_pileup.bdg文件的前三列排序，原文件的染色体顺序不定
~/local/bin/bedSort ${i}_treat_pileup.bdg ${i}_treat_pileup.sort.bdg
#根据染色体的最大坐标，修正部分记录的第三列（原文件的这些记录的第三列大于max_chr_pos，原因是前面在call peaks的时候reads有偏移）；
#这样一来，这些记录的第二列就比第三列大了，需要过滤掉
~/local/bin/bedClip -truncate itreatpileup.sort.bdg/ATACseq/align/chr.sizestdout∣perl−ane′printif({i}_treat_pileup.sort.bdg ~/ATAC_seq/align/chr.size stdout | perl -ane 'print if(itreatpileup.sort.bdg /ATACseq/align/chr.sizestdout∣perl−ane′printif(F[1] < $F[2])’ > ${i}_treat_pileup.bedgrafh
#以上两步的目的是生成一个“可读性”良好的bedgrafh文件，与原bedgrafh文件有区别
#bedgraph to bigwig
~/local/bin/bedGraphToBigWig ${i}_treat_pileup.bedgrafh ~/ATAC_seq/align/chr.size ${i}_treat_pileup.bw
done
8. 可视化
8.1 TSS Enrichment
TSS富集结果依赖于基因注释结果。这个图在学习ChIP-seq的时候也画过，方法一样

cd ~/ATAC_seq/visual
for i in SRR58746{57…62}
do
~/miniconda3/bin/computeMatrix reference-point -S ~/ATAC_seq/peak/${i}_treat_pileup.bw -R ~/Ref/Arabidopsis/Arabidopsis_thaliana.TAIR10.44.gtf -a 3000 -b 3000 -p 1 -o ${i}.matrix.TSS.3k.gz
~/miniconda3/bin/plotHeatmap -m ${i}.matrix.TSS.3k.gz -o ${i}.TSS.3k.png --colorMap Reds
done
8.2 IGV中查看数据
导入bigwig文件和narrowPeak文件到IGV中（可以显示peaks的名称），即可查看ATAC-Seq数据的结果。这一步之前也做过

多样品数据可重复评估
9.1 用bedtools计算overlap
wc -l SRR587465{7…9}_peaks.narrowPeak
30123 SRR5874657_peaks.narrowPeak
30644 SRR5874658_peaks.narrowPeak
30806 SRR5874659_peaks.narrowPeak
bedtools intersect -a SRR5874657_peaks.narrowPeak -b SRR5874658_peaks.narrowPeak -wa > tmp
bedtools intersect -a tmp -b SRR5874659_peaks.narrowPeak -wa | wc -l && rm -f tmp
29310

wc -l SRR587466{0…2}_peaks.narrowPeak
23965 SRR5874660_peaks.narrowPeak
26038 SRR5874661_peaks.narrowPeak
24812 SRR5874662_peaks.narrowPeak
bedtools intersect -a SRR5874660_peaks.narrowPeak -b SRR5874661_peaks.narrowPeak -wa > tmp
bedtools intersect -a tmp -b SRR5874662_peaks.narrowPeak -wa | wc -l && rm -f tmp
22773
由以上可知，重复性不错。

9.2 IDR评估
同时考虑peaks间的overlap和富集倍数的一致性。
使用和结果解读参考：https://www.jianshu.com/p/d8a7056b4294

#安装
~/miniconda3/bin/conda install idr
#先根据-log10(p-value)进行排序
for i in SRR58746{57…62}
do
sort -k8,8nr ~/ATAC_seq/peak/ipeaks.narrowPeak>/ATACseq/idr/{i}_peaks.narrowPeak > ~/ATAC_seq/idr/ipeaks.narrowPeak> /ATACseq/idr/{i}_peaks_sort.narrowPeak
done
#使用示例
~/miniconda3/bin/idr --samples SRR5874657_peaks_sort.narrowPeak SRR5874658_peaks_sort.narrowPeak -o SRR5874657_SRR5874658_idr --plot --input-file-type narrowP
eak --rank p.value &

作者：TOP生物信息
链接：https://www.jianshu.com/p/b7bd87315b6e
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

(rna) [zhangyong@cluster 16HD]$ ./bedClip
./bedClip: /lib64/libz.so.1: version ZLIB_1.2.3.3' not found (required by ./bedClip) ./bedClip: /lib64/libc.so.6: versionGLIBC_2.17’ not found (required by ./bedClip)
./bedClip: /lib64/libc.so.6: version `GLIBC_2.14’ not found (required by ./bedClip)
(rna) [zhangyong@cluster 16HD]$ Connection to 222.28.163.113 closed by remote host.

其他问题
如果还是继续报错，请注意报的so已经不一样了，把每一个so的问题都解决。
如果你依赖的gcc是4.8+，请用4.6版本，否则你要拥有root权限。

建议
这个东西根本原因就是系统版本库太老，如果更新glibc将会引起系统不稳定，所以建议还是升级系统吧。如果是centos，建议用6u3及以上版本。

作者：Kinva
链接：https://www.jianshu.com/p/308a4e803c81
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

H3K36me1_EE_rep1_FE.bdg
H3K36me1_EE_rep1_logLR.bdg

Method1

qrsh
cd /data/zhangyong/yyp/yypold/chips/lh11_3/data/X101SC20070420-Z01-J003/allcleandata/sam/chenalign/16HD
source activate rna

F=H3K36me1_EE_rep1_FE.bdg

G=mm10.len

bedtools slop -i ${F} -g ${G} -b 0 | ./bedClip stdin ${G} ${F}.clip

LC_COLLATE=C sort -k1,1 -k2,2n ${F}.clip > ${F}.sort.clip

./bedGraphToBigWig ${F}.sort.clip ${G} ${F/bdg/bw}

rm -f ${F}.clip ${F}.sort.clip

(rna) yyp@DESKTOP-LRUCLJJ:~$ F=H3K36me1_EE_rep1_FE.bdg
(rna) yyp@DESKTOP-LRUCLJJ:~$ G=mm10.len
(rna) yyp@DESKTOP-LRUCLJJ:~$ bedtools slop -i ${F} -g ${G} -b 0 | ./bedClip stdin ${G} F.clip(rna)yyp@DESKTOP−LRUCLJJ:{F}.clip (rna) yyp@DESKTOP-LRUCLJJ:~F.clip(rna)yyp@DESKTOP−LRUCLJJ: LC_COLLATE=C sort -k1,1 -k2,2n ${F}.clip > F.sort.clip(rna)yyp@DESKTOP−LRUCLJJ:{F}.sort.clip (rna) yyp@DESKTOP-LRUCLJJ:~F.sort.clip(rna)yyp@DESKTOP−LRUCLJJ: ./bedGraphToBigWig ${F}.sort.clip ${G} F/bdg/bw(rna)yyp@DESKTOP−LRUCLJJ:{F/bdg/bw} (rna) yyp@DESKTOP-LRUCLJJ:~F/bdg/bw(rna)yyp@DESKTOP−LRUCLJJ: ls
H3K36me1_EE_rep1_FE.bdg H3K36me1_EE_rep1_FE.bdg.sort.clip bedClip ee miniconda3 project
H3K36me1_EE_rep1_FE.bdg.clip H3K36me1_EE_rep1_FE.bw bedGraphToBigWig libz.so.1 mm10.len ref
(rna) yyp@DESKTOP-LRUCLJJ:~$

Method2
H3K36me1_EE_rep1_FE.bdg
H3K36me1_EE_rep1_logLR.bdg
DM24h-16HD-1_FKDL202607369.sort.bam

samtools view -H DM24h-16HD-1_FKDL202607369.sort.bam | perl -ne ‘if(/SN:(\S+)\s+LN:(\d+)/){print "$1\tKaTeX parse error: Undefined control sequence: \n at position 2: 2\̲n̲"}' > chr.size#…F[1] < $F[2])’ > H3K36me1_EE_rep1_FE.sort.bedgrafh
./bedGraphToBigWig H3K36me1_EE_rep1_FE.sort.bedgrafh mm10.len H3K36me1_EE_rep1_FE.1.bw

(rna) yyp@DESKTOP-LRUCLJJ:~$ ./bedSort H3K36me1_EE_rep1_FE.bdg H3K36me1_EE_rep1_FE.sort.bdg1
(rna) yyp@DESKTOP-LRUCLJJ:~$ ls

H3K36me1_EE_rep1_FE.bdg H3K36me1_EE_rep1_FE.bw bedGraphToBigWig libz.so.1 project
H3K36me1_EE_rep1_FE.bdg.clip H3K36me1_EE_rep1_FE.sort.bdg1 bedSort miniconda3 ref
H3K36me1_EE_rep1_FE.bdg.sort.clip bedClip ee mm10.len
(rna) yyp@DESKTOP-LRUCLJJ:~$ ./bedClip -truncate H3K36me1_EE_rep1_FE.sort.bdg1 mm10.len stdout | perl -ane 'print if($F[1] < F[2])′>H3K36me1EErep1FE.sort.bedgrafh(rna)yyp@DESKTOP−LRUCLJJ:F[2])' > H3K36me1_EE_rep1_FE.sort.bedgrafh (rna) yyp@DESKTOP-LRUCLJJ:~F[2])′>H3K36me1EErep1FE.sort.bedgrafh(rna)yyp@DESKTOP−LRUCLJJ: ./bedGraphToBigWig H3K36me1_EE_rep1_FE.sort.bedgrafh mm10.len H3K36me1_EE_rep1_FE.1.bw
(rna) yyp@DESKTOP-LRUCLJJ:~$

#求染色体序列大小
samtools view -H SRR5874657.final.bam | perl -ne ‘if(/SN:(\S+)\s+LN:(\d+)/){print “$1\t$2\n”}’ > chr.size ##等同下载的

bdg2bw(macs2的peak要修正坐标）相关推荐

川崎机器人here指令_川崎机器人定点修正坐标设置指导书.pdf
川崎机器人定点修正坐标设置指导书先按示教器上的[菜单]键然后打开[键盘]页然后把机器人移动到需要定点的位置上输入:here 变量这里的pick_t1 是位姿变量(示例) 然后确认输入请翻阅 ...
android 获取定位坐标，在百度地图出现大偏移的处理办法
问题现象: 通过android原生api获取到gps坐标,调用百度的地图api在手机上显示时,发现偏移距离超过1KM. android获取定位坐标的代码,网上很多不再粘贴. 问题分析: android ...
关于Delphi通用涵数
DELPHI程序注册码设计(转载) 思路是这样的:程序运行时先检测注册表,如果找到注册项,则表明已经注册,如果没有找到注册项,则提示要求注册. <注册例程> ...
转载：使用 OpenCV 识别 QRCode
原文链接:http://coolshell.cn/articles/10590.html#jtss-tsina 识别二维码的项目数不胜数,每次都是开箱即用,方便得很. 这次想用 OpenCV 从零识别 ...
pytorch 三维点分类_用于RGBD语义分割的三维图神经网络(2017ICCV,已开源)
3D Graph Neural Networks for RGBD Semantic Segmentation(2017ICCV, citation:78) 开源地址:https://github.c ...
如何采集指定年份的poi_房价关键影响因素分析：从数据采集到建模全过程
决定房价高低的关键性因素分析房价数据从采集到建模的详细过程以宁波为例作者:数字567作者简介:作者本人是CDA数据科学家,同时也是宁波校区的副校长本文共计6000字左右,如果不喜欢技 ...
javascript的拖放（第1部分）
既然说是入门,本文给出的函数都尽可能短小,但不失强大,并兼容所有浏览器.要想拖动页面上的一个元素,说白了就是让它实现位移.CSS中最能体现这种思想的是绝对定位,因为绝对定位能使元素脱离原来的文档流,但 ...
autojs 悬浮框演示代码
很多人在写autojs的时候,都是没有界面的,所以就少了很多智能化操作,这一节给大家分享一段代码,如何使用悬浮框. var 悬浮控制 = function(window, windowid, ar) ...
【目标检测】YOLO系列——YOLOv1详解
本篇是关于目标检测算法YOLOv1的学习笔记.网上的博客大多是摘抄翻译论文中的重点内容加上自己的理解,或者是梳理论文的大致思路,在没看过原文的情况下可能总会有些看不懂的地方.所以强烈建议博客搭配原文食 ...

bdg2bw(macs2的peak要修正坐标）

参数解释

-m FE 计算富集倍数降低噪音

-p 为了避免log的时候input值为0时发生error，给予一个很小的值

参数解释

check commands: slopBed, bedGraphToBigWig and bedClip

end of checking

bdg2bw(macs2的peak要修正坐标）相关推荐

最新文章

热门文章