ChIP-seq實踐（非轉錄因子，非組蛋白）

轉自 簡書 生信start_site

在實戰之前，首先對CHIP-seq分析做一些了解，以下是從各個地方copy過來綜合起來的，有些散亂，是我認為重要以及應該了解的知識點。
chip-seq 實戰就跟在轉錄組和外顯子上的實戰是沒有本質區別的，只是它多了一些分析，轉錄組只找表達量做差異分析，外顯子只找變異做變異注釋，那chip-seq 找的是peaks和motif 。
這里biobabble公眾號有比較完整的chip-seq的介紹：
CS1: ChIPseq簡介
CS2: BED文件
CS3: peak注釋
CS4：關于ChIPseq注釋的幾個問題
CS5: 吃著火鍋，唱著歌，還把分析給做了
CS6: ChIPseeker的可視化方法
CS7：Genomic coordination的富集性分析（1）
CS8：Genomic coordination的富集性分析（2）
CS9: GEO數據挖掘
CS番外1：如何定義下游的區間？

ChIP-seq 數據分析流程
1.得到原始序列文件后，第一步是執行標準的高通量數據質量控制。
2.原始的 reads 回比到研究物種的參考基因組上
3.特異性 ChIP-seq 的質量控制，檢查 ChIP-seq 樣品中的富集情況，并排除過度碎片
4.【核心】鑒定全基因上感興趣的因子富集的區域，也叫 Peak calling
5.一旦確定了峰值區域，還可以評估需要峰值位置的特定 ChIP-seq 的質量控制措施
6.在 Peak calling 之后，要確保預測的峰值準確地捕捉到binding pattern。
7.Peak calling后的分析取決于實驗設計與生物學問題。比如，如果在實驗中有重復和/或多個條件，那么下一步可能是樣本之間的比較（見第 8 章）。生物學重復將提供有關數據(噪聲)中的重現性和內在生物和技術變異性的信息。差異結合分析解決了哪一Peak區域在兩種情況下顯示出明顯不同的豐度（例如：與觀察到的相同條件的重復之間的變化相比，差異比預期的要大得多）。
8.Peak區間或差異Peak做下游分析，比如基因組注釋、GO分析、Pathway 分析、motif 查找、與其他基因組數據聯合分析。（參考文章：https://www.jianshu.com/p/19b178ea406e）

image.png

怎樣識別富集區域？(https://www.jianshu.com/p/df33bcb68548)
從DNA測序的角度來說，因為測序都是5'端的reads，對于一個DNA序列來說（有正負鏈的），它mapping的位置正負鏈都有的（也就是紅色和藍色的reads都有），對這些reads位置進行統計畫圖可以看到一個紅色的peak，一個藍色的peak。這兩個peak說明的是一個事情，就是這個地方有富集。最后對這兩個peak進行merge，最后變成了一個富集區域?；疑膒eak!

image

ChIP-seq的常見峰型
不同類型的蛋白或者組蛋白修飾會得到不同的峰形：
sharp binding sites, CTCF (red);
a mixture of shapes, RNA Polymerase II (orange);
medium size, H3K36me3 (green);
large domains, H3K27me3 (blue)

image

富集倍數
一般來說富集倍數要在5以上才算是顯著，如何計算富集比率呢？看圖~

image

Chip-Seq實踐

需要用的軟件：
fastqc, sra-tools, bowtie2, samtools, MACS2, deeptools
上述軟件之前在RNA-seq實踐中已經安裝了大部分，需要安裝MACS2和deeptools。
MACS2比較特殊，因為它需要Python2.7的環境才能安裝，而我的電腦安裝的是python3環境。所以首先創建一個虛擬環境，在python2.7的虛擬環境中使用（https://mp.weixin.qq.com/s/mEwl7-f2WTNfmK-FkNkP0A）：

conda create -n py27 python=2
source activate py27
pip install MACS2

本次實踐的文章是：
RYBP and Cbx7 define specific biological functions of polycomb complexes in mouse embryonic stem cells
（探索PRC1，PCR2蛋白復合物，不是轉錄因子或者組蛋白的CHIP-seq）
本文主要講了PRC1(Polycomb repressive complex 1)在小鼠的胚胎干細胞中有兩類亞型Cbx7-PRC1和RYBP-PRC1，但是Cbx7和RYBP這兩種是不能共存的，盡管兩者的全基因組定位在某些gene上交叉，也就是說在基因組上結合在了同一個位置。在分子水平上，Cbx7用于招募Ring1B結合到染色質上，然而RYBP可以增強PRC1的酶活性。RYBP結合的基因，其有著了較低水平的Ring1B和H2AK119ub，但其相比結合Cbx7有較高的表達量。在功能上，RYBP結合的基因主要涉及代謝調節和細胞周期，Cbx7結合的基因主要涉及early-lineage commitment of ESCs。（這段來源：http://www./archives/358）

1.下載原始數據：
可以寫一個腳本，批量下載：

#!/bin/bash
for ((i=204;i<=209;i++))
do wget https://sra-download.ncbi.nlm./traces/sra6/SRR/000605/SRR620$i
done

2.下載之后提取fastq文件：

#!/bin/bash
for i in SRR62020*
do
        echo $i
        fastq-dump --gzip --split-files $i
done

3.fastqc檢查測序質量
fastqc發現Ring1B、Suz12,IgG_old，IgG樣本的3’端有大約5個bp長度的堿基質量不是太好，可以在之后的bowtie2比對的時候設置參數，把質量差的幾個堿基去掉。

4.bowtie2比對
首先下載bowtie2的mm10基因組索引：

 wget ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip

還需要mm10 refseq注釋bed文件，進行網址下載http://genome./cgi-bin/hgTables，其中track和table選擇RefSeq，然后輸出格式為bed即可，最后下載（可以先在自己電腦上下載再傳到服務器上，也可以像下面直接在服務器上下載）（參考http://www./archives/358）

curl 'http://genome./cgi-bin/hgTables?hgsid=646311755_P0RcOBvAQnWZSzQz2fQfBiPPSBen&boolshad.hgta_printCustomTrackHeaders=0&hgta_ctName=tb_ncbiRefSeq&hgta_ctDesc=table+browser+query+on+ncbiRefSeq&hgta_ctVis=pack&hgta_ctUrl=&fbQual=whole&fbUpBases=200&fbExonBases=0&fbIntronBases=0&fbDownBases=200&hgta_doGetBed=get+BED' >mm10.refseq.bed

用bowtie2進行比對（用腳本進行）：

#!/bin/bash

bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620204_1.fastq.gz | samtools sort -O bam -o ring1B.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620205_1.fastq.gz | samtools sort -O bam -o cbx7.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620206_1.fastq.gz | samtools sort -O bam -o Suz12.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620207_1.fastq.gz | samtools sort -O bam -o RYBP.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620208_1.fastq.gz | samtools sort -O bam -o IgGold.bam
bowtie2 -p 6 -3 5 --local -x /media/yanfang/FYWD/CHIPSEQ/bowtie2mm10/mm10 -U /media/yanfang/FYWD/CHIPSEQ/fqfile/SRR620209_1.fastq.gz | samtools sort -O bam -o IgG.bam
~                                                   

但是貌似對比率不高：

19687027 reads; of these:
  19687027 (100.00%) were unpaired; of these:
    7778437 (39.51%) aligned 0 times
    7766495 (39.45%) aligned exactly 1 time
    4142095 (21.04%) aligned >1 times
60.49% overall alignment rate
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
20710168 reads; of these:
  20710168 (100.00%) were unpaired; of these:
    2584870 (12.48%) aligned 0 times
    10914615 (52.70%) aligned exactly 1 time
    7210683 (34.82%) aligned >1 times
87.52% overall alignment rate
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
21551927 reads; of these:
  21551927 (100.00%) were unpaired; of these:
    7109852 (32.99%) aligned 0 times
    9455003 (43.87%) aligned exactly 1 time
    4987072 (23.14%) aligned >1 times
67.01% overall alignment rate
[bam_sort_core] merging from 4 files and 1 in-memory blocks...
39870646 reads; of these:
  39870646 (100.00%) were unpaired; of these:
    6312903 (15.83%) aligned 0 times
    20702502 (51.92%) aligned exactly 1 time
    12855241 (32.24%) aligned >1 times
84.17% overall alignment rate
[bam_sort_core] merging from 8 files and 1 in-memory blocks...
13291429 reads; of these:
  13291429 (100.00%) were unpaired; of these:
    5608796 (42.20%) aligned 0 times
    4663584 (35.09%) aligned exactly 1 time
    3019049 (22.71%) aligned >1 times
57.80% overall alignment rate
[bam_sort_core] merging from 2 files and 1 in-memory blocks...
31218866 reads; of these:
  31218866 (100.00%) were unpaired; of these:
    5370101 (17.20%) aligned 0 times
    15180536 (48.63%) aligned exactly 1 time
    10668229 (34.17%) aligned >1 times
82.80% overall alignment rate
[bam_sort_core] merging from 6 files and 1 in-memory blocks...

5.用MACS2進行call peak:
首先要進入python2.7的環境：

source activate py27

然后寫一個腳本：

#!/bin/bash
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/Suz12.bam -m 10 30 -p 1e-5 -f BAM -g mm -n Suz12 2>Suz12.masc2.log
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/ring1B.bam -m 10 30 -p 1e-5 -f BAM -g mm -n ring1B 2>ring1B.masc2.log
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/cbx7.bam -m 10 30 -p 1e-5 -f BAM -g mm -n cxb7 2>cxb7.masc2.log
macs2 callpeak -c /media/yanfang/FYWD/CHIPSEQ/sambam/IgGold.bam -t /media/yanfang/FYWD/CHIPSEQ/sambam/RYBP.bam -m 10 30 -p 1e-5 -f BAM -g mm -n RYBP 2>RYBP.masc2.log

-t/–treatment 輸入文件，支持很多格式，搭配-f/–format使用
-f/–format 設定輸入文件的格式，這里選擇bam格式
-c/–control 對照組（或者模擬的）數據，需要跟-t的輸出文件在相同目錄下
-n/–name 輸出文件（有很多個文件）的前綴
–outdir 輸出文件所在目錄（不設定的話就默認當前目錄了）
-g/–gsize 提供基因組的大小，程序有默認的幾個物種可以選hs,mm,ce,dm
-q/–qvalue 設定FDR的閾值，默認是0.05
-p/–pvalue 設定pvalue的閾值，如果參數設定了-p，則其會覆蓋參數-q的作用
-B/–bdg If this flag is on, MACS will store the fragment pileup, control lambda, -log10pvalue and -log10qvalue scores in bedGraph files.

參考文章（http://www./archives/363）

運行后得到的文件不只是上述提到的一類，還有如下格式：
1.NAME_peaks.xls: 以表格形式存放peak信息，雖然后綴是xls，但其實能用文本編輯器打開，和bed格式類似，但是以1為基，而bed文件是以0為基.也就是說xls的坐標都要減一才是bed文件的坐標。
2.NAME_peaks.narrowPeak： NAME_peaks.broadPeak 類似。后面4列表示為， integer score for display， fold-change，-log10 pvalue，-log10 qvalue，relative summit position to peak start。內容和NAME_peaks.xls基本一致，適合用于導入R進行分析。
3.NAME_summits.bed：記錄每個peak的peak summits，話句話說就是記錄極值點的位置。MACS建議用該文件尋找結合位點的motif。
4.NAME_model.r： 能通過$ Rscript NAME_model.r作圖，得到是基于你提供數據的peak模型。

運行之后，看一下每個樣本找到多少個peak（bed的格式，儲存了每個peaks的位置信息）:

$ wc -l *.bed     2754 cxb7_summits.bed     7193 ring1B_summits.bed      473 RYBP_summits.bed     6696 Suz12_summits.bed

本文作者在GEO給的PEAKS個數如下：
2754 GSE42466_Cbx7_peaks_10.txt
6982 GSE42466_Ring1b_peaks_10.txt
6872 GSE42466_RYBP_peaks_5.txt
8054 GSE42466_Suz12_peaks_10.txt
然而我查了幾篇對該文獻重復的文章，call peak的數目和文章也有較大差距。（https://cloud.tencent.com/developer/article/1055109
http://www./archives/363
https://mp.weixin.qq.com/s/_A0rHldzEgVk7bgwt457qQ?
https://mp.weixin.qq.com/s/mEwl7-f2WTNfmK-FkNkP0A）

6.結果注釋與可視化
可視化這里我查閱了很多文章，分成兩種，有些人用的是chpseeker包進行peak的可視化和注釋；而還有一些人用了deeptool軟件進行分析。我兩種都試了試。

首先我用的是Chipseeker包。
ChIPseeker的功能分為三類(https://www.jianshu.com/p/c83a38915afc):
(1) 注釋：提取peak附近最近的基因， 注釋peak所在區域
(2)比較：估計ChIP peak數據集中重疊部分的顯著性；整合GEO數據集，以便于將當前結果和已知結果比較
(3)可視化： peak的覆蓋情況；TSS區域結合的peak的平均表達譜和熱圖；基因組注釋；TSS距離；peak和基因的重疊。

Chip peaks coverage plot

#調用要用的包
library("ChIPseeker")
library("org.Mm.eg.db")
library("TxDb.Mmusculus.UCSC.mm10.knownGene")
library("clusterProfiler")
#定位文件位置
setwd("/media/yanfang/FYWD/CHIPSEQ/bed")
#讀入bed文件
cbx7 <- readPeakFile("cxb7_peaks.narrowPeak.bed")
ring1B <- readPeakFile("ring1B_peaks.narrowPeak.bed")
Suz12 <- readPeakFile("Suz12_peaks.narrowPeak.bed")
#查看不同的蛋白的peak在染色體上的分布，因為RYBP的peak幾乎沒有，所以沒有用這個蛋白做后續的畫圖處理
covplot(cbx7, weightCol="V5")
covplot(ring1B, weightCol="V5")
covplot(Suz12, weightCol="V5")

例：cbx7的peak在染色體上的分布

#查看在指定的染色體上的分布
covplot(cbx7,chrs=c("chr17", "chr18"))

cbx7在第17和18號染色體上的分布

熱圖（Heatmap of ChIP binding to TSS regions）

# 定義TSS上下游的距離
> promoter <- promoters(TxDb.Mmusculus.UCSC.mm10.knownGene, upstream = 3000, downstream = 3000)> list_peak <- list(cbx7, ring1B, Suz12)> tagMatrix <- lapply(list_peak, getTagMatrix, window=promoter)> tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

image.png

Average Profile of ChIP peaks binding to TSS region
用譜圖的形式來展示不同蛋白結合的強度：

> library(org.Mm.eg.db)
> library(TxDb.Mmusculus.UCSC.mm10.knownGene)
> library(VennDiagram)
> setwd("/media/yanfang/FYWD/CHIPSEQ/bed")
> cbx7 <- readPeakFile("cxb7_peaks.narrowPeak.bed")> ring1B <- readPeakFile("ring1B_peaks.narrowPeak.bed")> Suz12 <- readPeakFile("Suz12_peaks.narrowPeak.bed")> peaks <- list(cbx7=cbx7,ring1B=ring1B,Suz12=Suz12)
> txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene> promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)> tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)> plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))

image.png

還支持使用bootstrap估計置信區間：

plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)

image.png

注釋peak
ChIPseeker負責注釋的就是annotatePeak

#對于非向量的循環，推薦構建list，然后用lapply> peaks <- list(ring1B, cbx7, Suz12)#注釋。啟動子區域是沒有明確的定義的，所以你可能需要指定tssRegion，把基因起始轉錄位點的上下游區域來做為啟動子區域。> peakAnnoList <- lapply(peaks, annotatePeak, tssRegion=c(-2500,2500), TxDb=TxDb.Mmusculus.UCSC.mm10.knownGene, 
                       addFlankGeneInfo=TRUE, flankDistance=5000)#結果可以用as.GRanges或者as.data.frame查看#Granges是R語言中的數據結構，它主要用來存儲基因組中的坐標區間。> as.GRanges(peakAnnoList[[1]])> as.data.frame(peakAnnoList[[1]])#畫圖，不同蛋白的靶位點的注釋情況> plotDistToTSS(peakAnnoList)

可視化TSS區域的TF binding loci

韋恩圖

genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)vennplot(genes[1:2:])

韋恩圖

注釋的可視化：

#先注釋
peakAnno <- annotatePeak(cbx7, tssRegion=c(-3000, 3000),
                         TxDb=txdb, annoDb="org.Mm.eg.db")
#柱形圖
plotAnnoBar(peakAnno)
#vennpie
vennpie(peakAnno)
#餅圖
plotAnnoPie(peakAnno)

柱形圖

vennpie

餅圖

# 展示注釋信息的交集$ upsetplot(peakAnno)

image.png

7.deeptools的使用
BAM文件是SAM的二進制轉換版，應該都知道。那么bigWig格式是什么？bigWig是wig或bedGraph的二進制版，存放區間的坐標軸信息和相關計分(score)，主要用于在基因組瀏覽器上查看數據的連續密度圖，可用wigToBigWig從wiggle進行轉換。為什么要用bigWig? 主要是因為BAM文件比較大，直接用于展示時對服務器要求較大。因此在GEO上僅會提供bw,即bigWig下載，便于下載和查看。
（一）將bam文件轉換成bw(bigwig)文件
這一步需要用到的是deeptool里的bamCoverage工具：

#!/bin/bash
bamCoverage -b cbx7.bam -o cbx7.bw -p 10 --normalizeUsing RPKM
bamCoverage -b ring1B.bam -o ring1B.bw -p 10 --normalizeUsing RPKM
bamCoverage -b Suz12.bam -o Suz12.bw -p 10 --normalizeUsing RPKM
bamCoverage -b IgGold.bam -o IgGold.bw -p 10 --normalizeUsing RPKM

-b:指要進行操作的bam文件，后面跟的是bam文件名。
-o:指輸出的文件，后面跟的是要輸出的文件名。
-p: --numberOfProcessors，進程的數量。
--normalizeUsing: 標準化的方法。這個工具有RPKM, CPM, BPM, RPGC, None幾個不同的選擇，這里我選了RPKM。

（二）peak分布可視化
也是deeptool里的工具。為了統計全基因組范圍的peak在基因特征的分布情況，需要用到computeMatrix計算，用plotHeatmap以熱圖的方式對覆蓋進行可視化，用plotProfile以折線圖的方式展示覆蓋情況。
computeMatrix有兩種模式，scale-regions mode和reference-point mode（參考文章：https://www.jianshu.com/p/2b386dd437d3）
reference-point(relative to a point): 計算某個點的信號豐度。給定一個bed file，以某個點為中心開始統計信號（TSS/TES/center）。
scale-regions(over a set of regions): 把所有基因組區段縮放至同樣大小，然后計算其信號豐度。簡單來說會將給定bed file范圍內的結合信號做一個統計（指的是一段長度），并將基因長度統一scale到設定regionBdoyLength的長度，加上統計基因上游和下游Xbp的信號（beforeRegionStartLength參數和afterRegionStartLength參數）
不管是哪一種，有兩個參數是必須的，-S提供bigwig文件，-R提供基因的注釋信息。

image

如果想得到單獨的文件，就是說每個樣品都有單獨的matrix文件，可以用下面的代碼:

#!/bin/bash
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/cbx7.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o cbx7_matrix_TSS.gz
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/ring1B.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o ring1B_matrix_TSS.gz
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/Suz12.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o Suz12_matrix_TSS.gz
computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/IgGold.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o IgGold_matrix_TSS.gz

reference-point: 選擇模式，這里我選擇了第二種。
-p 10:指10個進程。我的電腦比較渣，所以沒有加這個參數。
--referencePoint TSS：選擇參考點，這里選擇轉錄起始位點。
-a 3000:下游3000bp
-b 3000:上游3000bp
-S:指bw文件，后面跟的是文件的位置。
-R:參考基因組的注釋文件，這里要求是bed格式的。后面跟的是文件的位置。
--skipZeros：Whether regions with only scores of zero should be included or not. Default is to include them.
-o:輸出文件，后面跟的是輸出文件的名字

如果想把四個樣品合在一起，后續畫圖也是merge在一起的，可以用下面的代碼：

$ computeMatrix reference-point --referencePoint TSS -a 2500 -b 2500 -S /media/yanfang/FYWD/CHIPSEQ/bw/*.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o merge_matrix_TSS.gz

畫基因的TSS附近的profile和heatmap圖

image.png

#heatmap圖$ plotHeatmap -m merge_matrix_TSS.gz -out merge_heatmap.png#將heatmap圖輸出為pdf格式$ plotHeatmap --dpi 720 -m merge_matrix_TSS.gz -out merge_heatmap.pdf --plotFileFormat pdf#輪廓圖plotProfile -m merge_matrix_TSS.gz -out merge_profile.png#將輪廓圖輸出為pdf格式$ plotProfile --dpi 720 -m merge_matrix_TSS.gz -out merge_profile.pdf --plotFileFormat pdf#如果想把每個樣品的輪廓圖merge上在一個圖片里，在代碼后面加上--perGroup$ plotProfile -m merge_matrix_TSS.gz -out merge_profile.png --perGroup
$ plotProfile --dpi 720 -m merge_matrix_TSS.gz -out merge_profile.pdf --plotFileFormat pdf --perGroup

熱圖

把幾個樣品的輪廓圖merge到同一張圖里

也可以整合所有的chipseq的bw文件，畫基因的genebody附近的profile和heatmap圖（圖就不放上來了）:

$ computeMatrix scale-regions -a 3000 -b 3000 -m 5000 -S /media/yanfang/FYWD/CHIPSEQ/bw/*.bw -R /media/yanfang/FYWD/CHIPSEQ/bw/mm10.bed --skipZeros -o merge_scaleregion_matrix_TSS.gz
$ plotHeatmap -m merge_scaleregion_matrix_TSS.gz -out merge_scaleregion_matrix_TSS.png
$ plotProfile -m merge_scaleregion_profile_TSS.gz -out merge_scaleregion_profile_TSS.png