哈佛刘小乐实验室出品的软件,可以跟MACS软件call到的peaks文件无缝连接,实现peaks的注释以及可视化分析
该软件的主页里面很清楚的介绍了它可以做的图,以及每个图的意义:http://liulab.dfci.harvard.edu/CEAS/usermanual.html
我这里简单讲一下该软件如何安装以及使用,我这里还是使用我们的CHIP-seq分析系列教程的测试数据。
## 首先安装软件,是一个python程序,非常好安装
## Download and install CEAScd ~/biosoftmkdir CEAS && cd CEAStar zxvf CEAS-Package-1.0.2.tar.gzcd CEAS-Package-1.0.2python setup.py install --user~/.local/bin/ceas --help
## 然后测试软件自带的数据
cd ~/biosoft/CEASmkdir testData && cd testData## Human CD4T+ H3K36me3 ChIP-Seq data#########################run CEAS in the default mode##########################$ ceas -g gdb -b bed -w wigwhere gdb, bed, and wig stands for a sqlite3 file with a gene annotation table and genome background annotation, a BED file with ChIP regions, and a WIG file with ChIP enrichment signal file, respectively.#########################example for test data : ####################
数据详情如下:
300M May 20 2009 H3K36me3.wig
68M May 20 2009 H3K36me3.wig.zip
110K May 20 2009 H3K36me3_MACS_pval1e-5_peaks.bed
43K May 20 2009 H3K36me3_MACS_pval1e-5_peaks.bed.zip
11M Dec 2 2010 hg18.refGene
用测试数据来测试我们的CEAS软件:
~/.local/bin/ceas --name=H3K36me3_ceas --pf-res=20 \--gn-group-names='Top 10%,Bottom 10%' -g hg18.refGene -b H3K36me3_MACS_pval1e-5_peaks.bed -w H3K36me3.wig
但是有个重点,如何获取wig文件: https://github.com/crazyhottommy/ChIP-seq-analysis
因为我这里只是用了bowtie比对了得到了bam文件,并没有用MACS软件
GSM1278641 Xu_MUT_rep1_BAF155_MUT SRR1042593
GSM1278642 Xu_MUT_rep1_Input SRR1042594
我用的是perl单行命令把bam文件转为wig格式,我这里就拿上面两个样本数据做例子:
samtools depth SRR1042593.sorted.bam | perl -ne 'BEGIN{ print "track type=print wiggle_0 name=SRR1042593 description=SRR1042593\n"}; ($c, $start, $depth) = split; if ($c ne $lastC) { print "variableStep chrom=$c span=10\n"; };$lastC=$c; next unless $. % 10 ==0;print "$start\t$depth\n" unless $depth<3' > SRR1042593.wigsamtools depth SRR1042594.sorted.bam | perl -ne 'BEGIN{ print "track type=print wiggle_0 name=SRR1042594 description=SRR1042594\n"}; ($c, $start, $depth) = split; if ($c ne $lastC) { print "variableStep chrom=$c span=10\n"; };$lastC=$c; next unless $. % 10 ==0;print "$start\t$depth\n" unless $depth<3' > SRR1042594.wig
通过上面的学习,我们学会了该软件的使用,就可以拿自己的数据来玩一玩了。
## 然后处理我们直接的数据mkdir annotation && cd annotationwget http://liulab.dfci.harvard.edu/CEAS/src/hg19.refGene.gz ; gunzip hg19.refGene.gz~/.local/bin/ceas --name=H3K36me3_ceas --pf-res=20 --gn-group-names='Top 10%,Bottom 10%' \-g hg19.refGene -b ../paper_results/GSM1278641_Xu_MUT_rep1_BAF155_MUT.peaks.bed -w ../rawData/SRR1042593.wig