做甲基化分析时,生成了bedGraph
文件,想将其转化为bigwig
文件格式,以固定窗口的形式展示其窗口里的甲基化程度,最开始想使用USCS的工具bedGraphToBigWig
转换格式,但是貌似还是不能设置滑动窗口,后面在论坛查到了方法,链接为:https://www.biostars.org/p/226051/
示例
我是用bismark
进行甲基化的提取,最后生成的文件为GSM1032058_1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz
,将其解压。
- step1:
用sort-bed
进行排序,其为软件BEDOPS
的套件
awk '{ print $1"\t"$2"\t"$3"\tid-"NR"\t"$4; }' GSM1032058_1_val_1_bismark_bt2_pe.deduplicated.bedGraph | sort-bed - > input.bed
- step2:
fetchChromSizes
是UCSC
工具套件
fetchChromSizes hg38 > hg38.bounds.unsorted.txt
awk '{ print $1"\t0\t"$2; }' hg38.bounds.unsorted.txt | sort-bed - > hg38.bounds.bed
- step3:
bedops --chop 1000 --stagger 100 hg38.bounds.bed | bedmap --faster --echo --mean --delim "\t" --skip-unmapped - input.bed > final.bed
最后生成的final.bed文件内容如下
chr1 10200 11200 0.000000
chr1 10300 11300 0.000000
chr1 10400 11400 0.000000
chr1 10500 11500 0.000000
chr1 25800 26800 50.000000
chr1 25900 26900 78.571429
chr1 26000 27000 78.571429
chr1 26100 27100 78.571429
chr1 26200 27200 78.571429
chr1 26300 27300 78.571429
chr1 26400 27400 78.571429
chr1 26500 27500 78.571429
chr1 26600 27600 78.571429
chr1 26700 27700 78.571429
chr1 26800 27800 90.000000
chr1 38800 39800 0.000000
chr1 38900 39900 0.000000
chr1 39000 40000 0.000000
chr1 39100 40100 0.000000
chr1 39200 40200 0.000000
chr1 39300 40300 0.000000
chr1 39400 40400 0.000000
chr1 39500 40500 0.000000
chr1 39600 40600 0.000000
chr1 39700 40700 0.000000
chr1 87400 88400 33.333333
chr1 87500 88500 33.333333
chr1 87600 88600 33.333333
chr1 87700 88700 50.000000
chr1 87800 88800 40.000000
chr1 87900 88900 28.571429
chr1 88000 89000 28.571429
最后将其转化为了窗口的1000的形式,值为窗口内的平均值。
导入IGV
看看
导入没有问题,这些为0的也显示了,可以在final.bed文件中踢除这些为0的数值,但是为什么bed文件中的1000窗口不是连续地?如果是未检测到,为啥会有为0的数值,或者为0的数值是因为四舍五入才出现的,其实是有值地?
===========================================================================================
更正
后面想将生成的final.bed
文件转换为bw格式,使用命令
$ bedGraphToBigWig final_sorted.bed hg38.bounds.unsorted.txt final.bw
Error - overlapping regions in bedGraph line 2 of final_sorted.bed
报错,bed文件有overlapping,细看final.bed文件,确实每隔100bp有overlap,是前面转换时参数的设置问题,更改step3中的命令如下:
bedops --chop 1000 hg38.bounds.bed | bedmap --faster --echo --mean --delim "\t" --skip-unmapped - input.bed > final_no_stagger.bed
生成的final_no_stagger.bed
文件为
chr1 10000 11000 0.000000
chr1 26000 27000 78.571429
chr1 39000 40000 0.000000
chr1 88000 89000 28.571429
chr1 121000 122000 46.428571
chr1 128000 129000 0.000000
chr1 135000 136000 80.000000
chr1 136000 137000 83.333333
chr1 147000 148000 50.000000
chr1 151000 152000 37.500000
现在没有overlap了,然后再转换为bw文件
bedGraphToBigWig final_no_stagger.bed hg38.bounds.unsorted.txt final.bw
将 final_no_stagger.bed和 final.bw导入IGV查看
image.png
还是bw文件高高低低地看着舒服些,而且值为0的区域bw文件也不会显示
网友评论