1 ChIP-Seq技术
- 1.1 概念
- 1.2 ChIP-seq技术原理
2 ChIP-Seq数据分析
- 2.1 数据下载
- 2.2 质量控制（data_assess）
- 2.3 比对到参考基因组（mapping_analysis）
- 2.4 搜峰（Peak_calling）
  - MACS2
    - 2.4.1 MACS2 核心: callpeak 用法
    - 2.4.2 callpeak 结果文件说明
    - 2.4.3 bdg file → wig file
- 2.5 峰注释（Peak_anno）
  - ChIPseeker

ChIP-Seq仅仅是第一个表观遗传学领域比较成熟的技术而已，目前还有很多其他的技术，比如说

DNA修饰: DNA甲基化免疫共沉淀技术（MeDIP), 目标区域甲基化，全基因组甲基化（WGBS)，氧化-重亚硫酸盐测序（oxBS-Seq),
TET辅助重亚硫酸盐测序（TAB-Seq）

RNA修饰: RNA甲基化免疫共沉淀技术（MeRIP)

蛋白质与核酸相互作用： RIP-Seq, ChIP-Seq, CLIP-Seq

还有最近比较火的 ATAC-Seq ATAC-seq 能干啥？(
http://www.biotrainee.com/thread-1218-1-1.html
)

1 ChIP-Seq技术

1.1 概念

染色质免疫共沉淀技术 （Chromatin Immunoprecipitation， ChIP
）也称结合位点分析法，研究体内蛋白质与DNA相互作用的一种方法，通常用于转录因子结合位点或组蛋白特异性修饰位点的研究。
将 ChIP 与 第二代测序技术 相结合的 ChIP-seq技术
，能高效的在全基因组范围内检测与组蛋白、转录因子等互作的DNA片段。

1.2 ChIP-seq技术原理

在生理状态下，把细胞内的DNA与蛋白质交联（Crosslink）后裂解细胞，分离染色体，通过超声或酶处理将染色质随机切割；
利用抗原抗体的特异性识别反应，将与目的蛋白相结合的DNA片段沉淀下来；
再通过反交联（Reverse crosslink）释放结合蛋白的DNA片段；
纯化；
测序获得DNA片段的序列，最后将这些DNA片段比对到对应的参考基因组上。
![这里写图片描述](https://img-
blog.csdn.net/20180814104707287?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTc0NTg1OA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)

2 ChIP-Seq数据分析

2.1 数据下载

GSE98149 （包含H3K9me3的全部阶段，H3K4me3和H3K27me3的zygote、E6.5 Epi、E6.5 Exe、E7.5
Epi、E7.5 Exe、E8.5 embryo、Esc）
Title：Reprogramming of H3K9me3-dependent heterochromatin during mammalian
early embryo development [ChIP-seq]
Organism：Mus musculus

    for ((i=594;i<=670;i++));dowget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP105/SRP105176/SRR5479$i/SRR5479$i.sra;done &
[/code]```code# sratookit： .sra 文件 → fastq文件ls *sra |while read id;do/home/chen/sratoolkit.2.8.2-ubuntu64/bin/fastq-dump --gzip --split-3 $id;done &

    # 下载小鼠参考基因组的 indexwget -c "ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip" &# 解压unzip mm10.zip &
[/code]##  2.2 质量控制（data_assess）```code# Fastqc 进行质控ls *fq | while read id; do fastqc -t 4 $id; done &# multiqc：质控结果批量查看multiqc *fastqc.zip --export &
[/code]```code## trimmomatic # 安装 trimmomaticwget -c http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip &unzip Trimmomatic-0.38.zip# 数据清理# -threads 设置多线程运行java -jar "/data/chen/biosoft/Trimmomatic-0.38/trimmomatic-0.38.jar" PE -threads 2 -phred33 \# 2个输入文件${name}_1.fq.gz ${name}_2.trim.fq.gz \# 4个输出文件${name}_R1.clean.fq.gz ${name}_R1.unpaired.fq.gz \${name}_R2.clean.fq.gz ${name}_R2.unpaired.fq.gz \# ILLUMINACLIP:去接头# "$adapter"/Exome.fa ：adapter 序列的 fasta 文件# 2：16 个碱基长度的种子序列中可以有 2 个错配# 30：采用回文模式时匹配得分至少为30 (约50个碱基)# 10：采用简单模式时匹配得分至少为10 (约17个碱基)ILLUMINACLIP:"$adapter"/Exome.fa:2:30:10 \# LEADING:3，从序列的开头开始去掉质量值小于 3 的碱基；# TRAILING:3，从序列的末尾开始去掉质量值小于 3 的碱基；# SLIDINGWINDOW:4:15，从 5' 端开始以 4 bp 的窗口计算碱基平均质量，# 如果此平均值低于 15，则从这个位置截断 read；# HEADCROP:<length> 在reads的首端切除指定的长度；# MINLEN:36， 如果 reads 长度小于 36 bp 则扔掉整条 read。LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 HEADCROP:10 MINLEN:36
[/code]![这里写图片描述](https://img-
blog.csdn.net/20180904092004409?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTc0NTg1OA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)
![这里写图片描述](https://img-
blog.csdn.net/20180904092011835?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTc0NTg1OA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)##  2.3 比对到参考基因组（mapping_analysis）Bowtie2 或 BWA```code# bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r>} -S [<hit>]# -p/--threads NTHREADS 设置线程数. Default: 1# -q reads 是 fastq 格式的# -x <bt2-idx> index 路径# -1 <m1> 双末端测序的 _1.fastq 路径。可以为多个文件，并用逗号分开；多个文件必须和 -2 <m2> 中制定的文件一一对应。# -2 <m2> 双末端测序的 _2.fastq 路径.# -U <r> 非双末端测序的 fastq 路径。可以为多个文件，并用逗号分开。# -S <hit> 输出 Sam 格式文件。# -3/--trim3 <int> 剪掉3'端<int>长度的碱基，再用于比对。(default: 0).# 用fastqc看了看数据质量，发现3端质量有点问题，我就用了-3 5 --local参数，# --local 如果fq文件是没有经过 trim 的，可以用局部比对执行 soft-clipping，加上参数--local 。该模式下对read进行局部比对, 从而, read 两端的一些碱基不比对，从而使比对得分满足要求.bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479594_1.fastq -2 SRR5479594_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/MII_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479595_1.fastq -2 SRR5479595_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/MII_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479597_1.fastq -2 SRR5479597_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/sperm_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479598_1.fastq -2 SRR5479598_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/sperm_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479596_1.fastq -2 SRR5479596_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/sperm_input.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479599_1.fastq -2 SRR5479599_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/zygote_input.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479605_1.fastq -2 SRR5479605_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/zygote_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479606_1.fastq -2 SRR5479606_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/zygote_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479607_1.fastq -2 SRR5479607_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/2cell_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479608_1.fastq -2 SRR5479608_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/2cell_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479609_1.fastq -2 SRR5479609_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/4cell_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479610_1.fastq -2 SRR5479610_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/4cell_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479611_1.fastq -2 SRR5479611_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/8cell_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479612_1.fastq -2 SRR5479612_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/8cell_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479613_1.fastq -2 SRR5479613_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/morula_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479614_1.fastq -2 SRR5479614_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/morula_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479615_1.fastq -2 SRR5479615_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ICM_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479616_1.fastq -2 SRR5479616_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ICM_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479617_1.fastq -2 SRR5479617_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TE_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479618_1.fastq -2 SRR5479618_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TE_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479623_1.fastq -2 SRR5479623_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ESC_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479624_1.fastq -2 SRR5479624_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/ESC_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479625_1.fastq -2 SRR5479625_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TSC_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479626_1.fastq -2 SRR5479626_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TSC_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479627_1.fastq -2 SRR5479627_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/TSC_rep3.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479628_1.fastq -2 SRR5479628_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Epi_input.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479634_1.fastq -2 SRR5479634_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Epi_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479635_1.fastq -2 SRR5479635_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Epi_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479636_1.fastq -2 SRR5479636_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Exe_input.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479642_1.fastq -2 SRR5479642_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Exe_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479643_1.fastq -2 SRR5479643_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E65Exe_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479644_1.fastq -2 SRR5479644_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_input.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479650_1.fastq -2 SRR5479650_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479651_1.fastq -2 SRR5479651_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479652_1.fastq -2 SRR5479652_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep3.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479653_1.fastq -2 SRR5479653_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Epi_rep4.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479654_1.fastq -2 SRR5479654_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Exe_input.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479661_1.fastq -2 SRR5479661_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Exe_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479662_1.fastq -2 SRR5479662_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E75Exe_rep2.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479663_1.fastq -2 SRR5479663_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E85Epi_input.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479669_1.fastq -2 SRR5479669_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E85Epi_rep1.sam &bowtie2 -p 6 -q -x /data/chen/Data/reference/index/mm10/mm10 -1 SRR5479670_1.fastq -2 SRR5479670_2.fastq -S /data/chen/Data/GSE98149/chip-seq/h3k9me3/alignment/E85Epi_rep2.sam &
[/code]![这里写图片描述](https://img-
blog.csdn.net/2018090411251642?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTc0NTg1OA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)
![这里写图片描述](https://img-
blog.csdn.net/20180904112404199?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTc0NTg1OA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)
![这里写图片描述](https://img-
blog.csdn.net/20180904161257837?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTc0NTg1OA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)
![这里写图片描述](https://img-
blog.csdn.net/2018090416131067?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MTc0NTg1OA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)##  2.4 搜峰（Peak_calling）###  **MACS2**peaks calling：寻找可能的结合位点，即基因组中大量reads富集的区域。####  2.4.1 MACS2 核心: callpeak 用法```code# Example for regular peak callingmacs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01# Example for broad peak callingmacs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
[/code]```code# 批量callpeakmacs2 callpeak -c IgGold.bam -t suz12.bam -q 0.05 -f BAM -g mm -n suz12 2> suz12.macs2.log &macs2 callpeak -c IgGold.bam -t ring1B.bam -q 0.05 -f BAM -g mm -n ring1B 2> ring1B.macs2.log &macs2 callpeak -c IgGold.bam -t cbx7.bam -q 0.05 -f BAM -g mm -n cbx7 2> cbx7.macs2.log &macs2 callpeak -c IgGold.bam -t RYBP.bam -q 0.01 -f BAM -g mm -n RYBP 2>RYBP.macs2.log &
[/code]**-t** /–treatment FILENAME——处理组输入
This is the only REQUIRED parameter for MACS. File can be in any supported
format specified by –format option. Check –format for detail. If you have more
than one alignment files, you can specify them as ` -t A B C ` . MACS will
pool up all these files together.**-c** /–control——对照组输入
The control or mock data file. Please follow the same direction as for
-t/–treatment.**-f** /–format FORMAT——-t和-c提供文件的格式，目前MACS能够识别的格式有 “ELAND”, “BED”,
“ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (双端测序), “SAM”, “BAM”, “BOWTIE”,
“BAMPE”, “BEDPE”. 除”BAMPE”, “BEDPE”需要特别声明外，其他格式都可以用 AUTO自动检测。如果不提供这项，就是自动检测选择。
Format of tag file, can be “ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”,
“ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” or
“BEDPE”. Default is “AUTO” which will allow MACS to decide the format
automatically. “AUTO” is also usefule when you combine different formats of
files. Note that MACS can’t detect “BAMPE” or “BEDPE” format with “AUTO”, and
you have to implicitly specify the format for “BAMPE” and “BEDPE”.**-g** /–gsize——基因组大小，默认提供了hs, mm, ce,
dm选项，不在其中的话，比如说拟南芥，就需要自己提供了（拟南芥根据NCBI显示是119,667,750，也就是1.2e8）。
PLEASE assign this parameter to fit your needs!
It’s the mappable genome size or effective genome size which is defined as the
genome size which can be sequenced. Because of the repetitive features on the
chromsomes, the actual mappable genome size will be smaller than the original
size, about 90% or 70% of the genome size. The default hs – 2.7e9 is
recommended for UCSC human hg18 assembly. Here are all precompiled parameters
for effective genome size:
hs: 2.7e9 （人类是2.7e9，也就是2.7G）
mm: 1.87e9
ce: 9e7
dm: 1.2e8**-n** /–name——输出文件的前缀名。表示实验的名字, 请取一个有意义的名字。
The name string of the experiment. MACS will use this string NAME to create
output files like ‘NAME_peaks.xls’, ‘NAME_negative_peaks.xls’,
‘NAME_peaks.bed’ , ‘NAME_summits.bed’, ‘NAME_model.r’ and so on. So please
avoid any confliction between these filenames and your existing files.**-B** /–bdg
会保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and
-log10qvalue scores。
If this flag is on, MACS will store the fragment pileup, control lambda,
-log10pvalue and -log10qvalue scores in bedGraph files. The bedGraph files
will be stored in current directory named NAME+’_treat_pileup.bdg’ for
treatment data, NAME+’_control_lambda.bdg’ for local lambda values from
control, NAME+’_treat_pvalue.bdg’ for Poisson pvalue scores (in -log10(pvalue)
form), and NAME+’_treat_qvalue.bdg’ for q-value scores from
Benjamini–Hochberg–Yekutieli procedure [
http://en.wikipedia.org/wiki/False_discovery_rate#Dependent_tests
](http://en.wikipedia.org/wiki/False_discovery_rate#Dependent_tests)**-q** /–qvalue——q值，也就是最小的PDR阈值， 默认是0.05。q值是根据p值利用BH计算，也就是多重试验矫正后的结果。
The qvalue (minimum FDR) cutoff to call significant regions. Default is 0.01.
For broad marks, you can try 0.05 as cutoff. Q-values are calculated from
p-values using Benjamini-Hochberg procedure.**-p** /–pvalue——p值，指定 p 值后 MACS2 就不会用 q 值了。
The pvalue cutoff. If -p is specified, MACS2 will use pvalue instead of
qvalue.**-m**
/–mfold——和MFOLD有关，而MFOLD和MACS预构建模型有关，默认是5：50，MACS会先寻找100多个peak区构建模型，一般不用改，因为你不懂。
This parameter is used to select the regions within MFOLD range of high-
confidence enrichment ratio against background to build model. The regions
must be lower than upper limit, and higher than the lower limit of fold
enrichment. DEFAULT:5,50 means using all regions not too low (>5) and not too
high (<50) to build paired-peaks model. If MACS can not find more than 100
regions to build model, it will use the –extsize parameter to continue the
peak detection ONLY if –fix-bimodal is set.####  2.4.2 callpeak 结果文件说明```code# （在当前目录下）统计 *bed 的行数（peak数）wc -l *bed2384 cbx7_summits.bed8342 ring1B_summits.bed0 RYBP_summits.bed7619 suz12_summits.bed# 在文件a中统计 hello 出现的行数:# grep hello a | wc -l# wc（word count）#-c 统计字节数。#-l 统计行数。line#-m 统计字符数。这个标志不能与 -c 标志一起使用。#-w 统计字数。一个字被定义为由空白、跳格或换行字符分隔的字符串。#-L 打印最长行的长度。
[/code]callpeak会得到如下结果文件：**NAME_summits.bed** ：Browser Extensible Data，记录每个peak的peak
summits，话句话说就是记录极值点的位置。MACS建议用该文件寻找结合位点的motif。能够直接载入UCSC
browser，用其他软件分析时需要去掉第一行。**NAME_peaks.xls** ：以表格形式存放peak信息，虽然后缀是xls，但其实能用文本编辑器打开，和bed格式类似，但是 **以1为基**
，而bed文件是以0为基.也就是说xls的坐标都要减一才是bed文件的坐标。**NAME_peaks.narrowPeak** ， **NAME_peaks.broadPeak** 类似。后面4列表示为， integer score
for display， fold-change，-log10pvalue，-log10qvalue，relative summit position to
peak start。内容和NAME_peaks.xls基本一致，适合用于导入R进行分析。**NAME_model.r** ：能通过 ` $ Rscript NAME_model.r ` 作图，得到是基于你提供数据的peak模型。**.bdg** ：能够用 UCSC genome browser 转换成更小的 bigWig 文件。####  2.4.3 bdg file → wig file> 为了方便在IGV上查看ChIP-seq的结果和后期的可视化展示，需要把macs2的结果转化为bw提供给IGV。一共分为三步：第一步：使用 bdgcmp 得到 **FE** 或者 **logLR 转化后的文件** (Run MACS2 bdgcmp to generate
fold-enrichment and logLR track)```codemacs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_FE.bdg -m FEmacs2 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，给予一个很小的值
[/code]第二步： **预处理** 文件，下载对应 **参考基因组染色体长度** 文件使用conda安装以下三个处理软件：
bedGraphToBigWig
bedClip
bedtools下载染色体长度文件： [
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz
](http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz)
并解压（针对human，其余物种的皆可以按照类似网址下载）写一个小小的sh脚本方便一步转化name.sh：```code#!/bin/bash# check commands: slopBed, bedGraphToBigWig and bedClipwhich 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 checkingif [ $# -lt 2 ];thenecho "Need 2 parameters! <bedgraph> <chrom info>"exitfiF=$1G=$2bedtools slop -i ${F} -g ${G} -b 0 | bedClip stdin ${G} ${F}.clipLC_COLLATE=C sort -k1,1 -k2,2n ${F}.clip > ${F}.sort.clipbedGraphToBigWig ${F}.sort.clip ${G} ${F/bdg/bw}rm -f ${F}.clip ${F}.sort.clip
[/code]chmod +x name.sh第三步：生成 **bw** 文件```code./name.sh H3K36me1_EE_rep1_FE.bdg hg19.len./name.sh H3K36me1_EE_rep1_logLR.bdg hg19.len
[/code]最后得到产物，至于的使用哪一个作为输入文件大家就根据需要来吧H3K36me1EErep1_FE.bw
H3K36me1EErep1_logLR.bw##  2.5 峰注释（Peak_anno）###  ChIPseeker> ChIPseeker的功能分为三类：
>  注释：提取peak附近最近的基因，注释peak所在区域。
>  比较：估计ChIP peak数据集中重叠部分的显著性；整合GEO数据集，以便于将当前结果和已知结果比较。
>  可视化：peak的覆盖情况；TSS区域结合的peak的平均表达谱和热图；基因组注释；TSS距离；peak和基因的重叠。
```code# 加载ChIPseeker、基因组注释包和bed数据biocLite("ChIPseeker")biocLite("org.Mm.eg.db")biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")library("ChIPseeker")# 下载source ("https://bioconductor.org/biocLite.R")biocLite("ChIPseeker")biocLite("org.Mm.eg.db")biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")biocLite("clusterProfiler")biocLite("ReactomePA")biocLite("DOSE")#载入library("ChIPseeker")library("org.Mm.eg.db")library("TxDb.Mmusculus.UCSC.mm10.knownGene")txdb <- TxDb.Mmusculus.UCSC.mm10.knownGenelibrary("clusterProfiler")# 读入bed文件ring1B <- readPeakFile("F:/Data/ChIP/ring1B_peaks.narrowPeak")cbx7 <- readPeakFile("F:/Data/ChIP/cbx7_peaks.narrowPeak")RYBP <- readPeakFile("F:/Data/ChIP/RYBP2_summits.bed")suz12 <- readPeakFile("f:/Data/ChIP/suz12_peaks.narrowPeak")
[/code]Chip peaks coverage plot查看peak在全基因组的位置covplot(ring1B,chrs=c(“chr17”, “chr18”)) #specific chrring1Bsuz12cbx7RYBPring1B(chr17&18)● Heatmap of ChIP binding to TSS regionspromoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(ring1B, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color=”red”)Average Profile of ChIP peaks binding to TSS region● Confidence interval estimated by bootstrap methodplotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)peak的注释peak的注释用annotatePeak**，TSS (transcription start site) region
可以自己设定，默认是（-3000，3000），TxDb 是指某个物种的基因组，例如TxDb.Hsapiens.UCSC.hg38.knownGene,
TxDb.Hsapiens.UCSC.hg19.knownGene for human genome hg38 and hg19,
TxDb.Mmusculus.UCSC.mm10.knownGene and TxDb.Mmusculus.UCSC.mm9.knownGene for
mouse mm10 and mm9.peakAnno <- annotatePeak(ring1B, tssRegion=c(-3000, 3000),
TxDb=txdb, annoDb=”org.Mm.eg.db”)可视化 Pie and Bar plotplotAnnoBar(peakAnno)
vennpie(peakAnno)
upsetplot(peakAnno)饼图:条形图：upsetplot: upset技术适用于多于5个集合的表示情况。可视化TSS区域的TF binding lociplotDistToTSS(peakAnno,
title=”Distribution of transcription factor-binding loci\nrelative to TSS”)多个peak的比较多个peak set注释时，先构建list,然后用lapply.
list(name1=bed_file1,name2=bed_file2)RYBP的数据有问题，这里加上去，会一直报错。peaks <- list(cbx7=cbx7,ring1B=ring1B,suz12=suz12)
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf=0.95,resample=500,
facet=”row”)
tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=NULL)ChIP peak annotation comparisionpeakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb,
tssRegion=c(-3000, 3000), verbose=FALSE)
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)Overlap of peaks and annotated genesgenes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)![在这里插入图片描述](https://img-blog.csdnimg.cn/20210608151750993.gif)

ChIP-seq 数据分析1 ChIP-Seq技术2 ChIP-Seq数据分析相关推荐

学习全基因组测序数据分析1：测序技术
本文转载自微信公众号解螺旋的矿工,作者为黄树嘉,已获得授权.黄树嘉写了WGS系列的文章,堪称教科书级别的生物信息学习材料.虽然本平台只关注宏基因组领域,但此系列文章知识体系完善.干货满满,是值得每位专 ...
Interview：算法岗位面试—10.11下午—上海某公司算法岗位(偏数据分析，证券金融行业)技术面试考点之sqlserver语言相关考察点复习
Interview:算法岗位面试-10.11下午-上海某公司算法岗位(偏数据分析,证券金融行业)技术面试考点之sqlserver语言相关考察点复习导读:其实,考察的知识点,博主都做过, 还包括sql ...
Interview：算法岗位面试—上海某公司算法岗位(偏数据分析，互联网行业)技术面试考点之特征工程考察点
ML岗位面试:上海某公司算法岗位(偏数据分析,互联网行业)技术面试考点之特征工程考察点 Interview:算法岗位面试-上海某公司算法岗位(偏数据分析,互联网行业)技术面试考点之特征工程考察点导读 ...
中国临床数据分析市场趋势报告、技术动态创新及市场预测
出版商:贝哲斯咨询获取报告样本: 企业竞争态势临床数据分析市场报告涉及的主要国际市场参与者有Athenahealth.Cerner.McKesson.Xerox等.这些参与者的市场份额.收入.公司 ...
大数据分析处理及挖掘技术
数据处理是对纷繁复杂的海量数据价值的提炼,而其中最有价值的地方在于预测性分析,即可以通过数据可视化.统计模式识别.数据描述等数据挖掘形式帮助数据科学家更好的理解数据,根据数据挖掘的结果得出预测性决策. ...
python基于爬虫技术的海量电影数据分析源码,数据处理分析可视化，GUI界面展示
基于爬虫技术的海量电影数据分析介绍一个基于爬虫技术的海量电影数据分析系统系统架构本系统主要分为四个部分,分别为后端爬虫抓取.数据处理分析可视化.GUI界面展示.启动运行,分别对应getData ...
郑大研究生计算机科学与技术,21郑大考研计算机科学与技术、软件工程考研数据分析...
原标题:21郑大考研计算机科学与技术.软件工程考研数据分析 [天任郑大考研网]21郑大考研计算机科学与技术.软件工程(408计算机科专业基础综合)考研数据分析 1.专业代码及专业名称: 081200计 ...
基于大数据分析的安全管理平台技术研究及应用
http://www.venustech.com.cn/NewsInfo/531/25566.Html [内容摘要]本文首先通过介绍大数据的起因,给出了大数据的定义和特征描述,并简要说明了当前大数据的 ...
数据分析的方法与技术
数据分析的方法与技术数据分析是指采用准确适宜的分析方法和工具来分析经过处理的数据,提取有价值的信息,从而形成有效的结论并通过可视化技术展现出来的过程. 数据分析的方法有: 基本分析方法:主要以基础的 ...
高性能计算系统——大数据/快速数据分析中的高性能技术
大数据/快速数据分析中的高性能技术高性能计算的目的是为了数据密集型以及处理密集型的工作实现少费而多用的目标.计算机.存储设备和网络解决方案也相应变得高性能和可扩展. 高通量计算(HTC)同高性能计算 ...

ChIP-seq 数据分析1 ChIP-Seq技术2 ChIP-Seq数据分析