deeptools？还可以这么用！

    说到deeptools，大家很容易想到ChIP。deeptools集美貌与才华 是表观分析的深海利器。关于deeptools的应用，其中最为经典的就是plotHeatmap和plotProfiles（如下图）。几乎每个经典的表观文章都能找到他们的身影。

经典的用法主要是研究某种组蛋白或者转录因子是不是在TSS区域富集（如下图）。

但殊不知deeptools不仅仅是查看组蛋白或者转录因子在基因区域的修饰，还能为三维基因组理论奠定基础添砖加瓦，例如EP_loop（enhancer_promoter构成的loop）的CTCF以及蛋白修饰情况（如下图）：

例如，TAD边界以及ABcompartment两端的组蛋白修饰情况（如下图）：

说到这里，小编的手机响了，收到 一条信息～～

Hi～～fan，原来deeptools还可以用作HiC多组学分析呀，那我想拿TAD边界的组蛋白修饰和AB compartment组蛋白修饰练练手，我具体该怎么做呢？

哈哈，别急，容我一一道来～～

首先看图～～ 

TAD边界这幅图呢，主要是以一个点为中心，而ABcompart这幅图则是有两个中心点border1，border2（箭头所指）

这就涉及deeptools的computeMatrix两种模式（如下图）：

deeptools  computeMatrix模块可以计算每个基因区域的结合得分，生成中间文件用以给plotHeatmap和plotProfiles作图。

reference-point mode则是给定一个bed file，以某个点为中心开始统计信号（如TSS/TES/center）。

来举个'🌰'：

首先，整理一下TAD边界的bed文件boundaries.bed  

其次，第二步通过以下命令将bam转换成bigwig：

    bamCoverage --bam  CTCF.bam -o CTCF.bw

    --binSize 10

    --numberOfProcessors 4 

   --normalizeUsingRPKM

    --extendReads

之后，第三步利用reference-point 计算以这个点为中心上下游500k范围内以10k为窗口计算每个窗口CHIP的信号强度：

computeMatrix reference-point  --referencePoint center 

       -b 500000 -a 500000 

       -R boundaries.bed

        -S CTCF.bw H3K27ac.bw H3K27me3.bw H3K4me3.bw --skipZeros

        -o  boundaries_modification.gz 

       --binSize 10000

然后，第四步画线图：

plotProfile -m boundaries_modification.gz 

    -out boundaries_modification.pdf 

      --plotType se 

      --plotFileFormat 'pdf' 

    --refPointLabel Boundary --numPlotsPerRow 4 

     --perGroup 

    --plotHeight 8 --plotWidth 8 -

    --legendLocation lower-center 

    --yMin 0 --plotTitle ""

     --yAxisLabel "Normalized signal"

经过以上四步就可以收工了，快来看结果，怎么都感觉有点小清新呢～

scale-regions mode简单来说会将给定bed file范围内的结合信号做一个统计（指的是一段长度），并将基因长度统一scale到设定regionBdoyLength的长度，加上统计基因上游和下游（通过beforeRegionStartLength参数和afterRegionStartLength参数设定）n bp的信号。

同样的针对ABcompartment先要识别连续的ABcomaprtment位点形成chr，start，end三列bed文件(A_com.bed ,B_com.bed)。

之后scale-regions  计算当前以及延伸区域的信号值～

computeMatrix scale-regions 

    --regionsFileName A_com.bed  B_com.bed

     --scoreFileName  H3K27ac.bw   H3K27me3.bw  H3K36me3.bw   H3K4me3.bw              H3K9me3.bw  

     --skipZeros 

    --regionBodyLength 6000000 

    --beforeRegionStartLength 2000000 

    --afterRegionStartLength 2000000 

    --startLabel border --endLabel border 

    -outFileName  AB_modification.gz  

     --binSize 100000 

统计完信号之后，就可以画图了：

plotProfile  -m   AB_modification.gz  

    -out  AB_mod_lineProfile.pdf 

    --plotType  se 

    --startLabel border 

    --endLabel border 

    --plotFileFormat  pdf 

   --plotHeight 8 --plotWidth 12

    --yAxisLabel "Reads Density"

     --perGroup

咚咚咚～～，通过这样简单的步骤，图就画好了，惊不惊喜？意不意外？

艾玛，就扯到这儿吧～～，朕累了～～