安装默认使用Conda
macs2的基本命令
macs2 callpeak -t treatment.bam -c control.bam -g hs -B -f BAM -n prefix -q 0.00001
macs2 callpeak -t treatment.bam -c control.bam -g hs -B -f BAMPE -n prefix -q 0.00001
参数解读：
-g 基因组的选择
-B 输出bgd文件，下游bigwig文件生成所需
-f 双端测序使用BAMPE，单端的话不需要加参数，默认是auto识别，但是识别不了BAMPE
-q 使用阈值，根据需要选择
输出文件包括
prefix_model.r
prefix_peaks.narrowPeak
prefix_peaks.xls
prefix_summits.bed
prefix_treat_pileup.bdg
prefxi_ _treat_pvalue.bdg
prefix_model.r 为模型建立示意图，使用Rscript prefix_model.r 命令可以在Linux环境下作图
prefix_peaks.xls 为包含了peak信息的表格，可以用excel打开，其中包括：
染色体
peak起始位点
peak终止位点
peak区域的长度
peak定点的绝对位置
堆积信号值及其-log10(pvalue)
fold enrichment及其-log10(pvalue)
#注意：XLS文件里面的起始位点和GFF文件相似从1起始，bed文件从0开始
prefix_peaks.narrowPeak 为包含了peak位置及summit位置，p值，q值的文件，可以直接上传到UCSC browser，其中特定列的信息为：
5th：整数信号值
7th：fold-change
8th：-log10 pvalue
9th：-log10 qvalue
10th：summit距离peak起始位点的距离
prefix_summits.bed 包含每个peak峰顶summit的位置，等同于从narrowPeak文件里面剥离出来的，方便用来寻找motif的binding，同样可以载入UCSC
prefix_peaks.broadPeak &prefix_peaks.gappedPeak 文件都是为选择 -broad 参数后得到的文件，本质上和narrowPeak文件一致，不过增加了broadregion，当然summit的定义也有区别，需要自己去指定
prefix_treat_pileup.bdg 文件为输入文件treat组的bedGraph file， prefxi_ _treat_pvalue.bdg为对照组的bedGraph file
为了方便在IGV上查看ChIP-seq的结果和后期的可视化展示，所以我们需要把macs2的结果转化为bw提供给IGV(bdg file to wig file transformation)
一共分为三步
1.首先使用bdgcmp得到FE或者logLR转化后的文件(Run MACS2 bdgcmp to generate fold-enrichment and logLR track)
macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_FE.bdg -m FE
macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_logLR.bdg -m logLR -p 0.00001
# 参数解释
-m FE 计算富集倍数降低噪音
-p 为了避免log的时候input值为0时发生error，给予一个很小的值
2. 预处理文件与对应参考基因组染色体长度文件下载
使用conda安装以下三个处理软件：bedGraphToBigWig/bedClip/bedtools
下载染色体长度文件：http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz 并解压（针对human，其余物种的皆可以按照类似网址下载）
写一个小小的sh脚本方便一步转化name.sh：
#!/bin/bash
# check commands: slopBed, bedGraphToBigWig and bedClip
which bedtools &>/dev/null || { echo "bedtools not found! Download bedTools: <http://code.google.com/p/bedtools/>"; exit 1; }
which bedGraphToBigWig &>/dev/null || { echo "bedGraphToBigWig not found! Download: <http://hgdownload.cse.ucsc.edu/admin/exe/>"; exit 1; }
which bedClip &>/dev/null || { echo "bedClip not found! Download: <http://hgdownload.cse.ucsc.edu/admin/exe/>"; exit 1; }
# end of checking
if [ $# -lt 2 ];then
echo "Need 2 parameters! <bedgraph> <chrom info>"
exit
fi
F=$1
G=$2
bedtools slop -i ${F} -g ${G} -b 0 | bedClip stdin ${G} ${F}.clip
LC_COLLATE=C sort -k1,1 -k2,2n ${F}.clip > ${F}.sort.clip
bedGraphToBigWig ${F}.sort.clip ${G} ${F/bdg/bw}
rm -f ${F}.clip ${F}.sort.clip
chmod +x name.sh
3. 最后就是生成bw文件
./name.sh H3K36me1_EE_rep1_FE.bdg hg19.len
./name.sh H3K36me1_EE_rep1_logLR.bdg hg19.len
最后得到产物，至于的使用哪一个作为输入文件大家就根据需要来吧
H3K36me1_EE_rep1_FE.bw
H3K36me1_EE_rep1_logLR.bw
参考文献：
网友评论