滑动窗口展示甲基化程度

作者: scdzzdw | 来源:发表于2019-11-05 15:57 被阅读0次

    做甲基化分析时,生成了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:
      fetchChromSizesUCSC工具套件
    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看看

    image.png
    导入没有问题,这些为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文件也不会显示

    相关文章

      网友评论

        本文标题:滑动窗口展示甲基化程度

        本文链接:https://www.haomeiwen.com/subject/iqtebctx.html