當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

ChIP-seq 数据分析

發布時間：2023/12/14 编程问答 31 豆豆

生活随笔收集整理的這篇文章主要介紹了 ChIP-seq 数据分析小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

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片段比對到對應的參考基因組上。

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++)); do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP105/SRP105176/SRR5479$i/SRR5479$i.sra; done & # sratookit： .sra 文件 → fastq文件 ls *sra |while read id; do /home/chen/sratoolkit.2.8.2-ubuntu64/bin/fastq-dump --gzip --split-3 $id; done & # 下載小鼠參考基因組的 index wget -c "ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip" &# 解壓 unzip mm10.zip &

2.2 質量控制（data_assess）

# Fastqc 進行質控 ls *fq | while read id; do fastqc -t 4 $id; done &# multiqc：質控結果批量查看 multiqc *fastqc.zip --export & ## trimmomatic # 安裝 trimmomatic wget -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

2.3 比對到參考基因組（mapping_analysis）

Bowtie2 或 BWA

# 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 &

2.4 搜峰（Peak_calling）

MACS2

peaks calling：尋找可能的結合位點，即基因組中大量reads富集的區域。

2.4.1 MACS2 核心: callpeak 用法

# 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 # 批量callpeak macs2 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 &

-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

-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 結果文件說明

# （在當前目錄下）統計 *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 打印最長行的長度。

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)

macs2 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，給予一個很小的值

第二步：預處理文件，下載對應參考基因組染色體長度文件

使用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

第三步：生成 bw 文件

./name.sh H3K36me1_EE_rep1_FE.bdg hg19.len ./name.sh H3K36me1_EE_rep1_logLR.bdg hg19.len

最后得到產物，至于的使用哪一個作為輸入文件大家就根據需要來吧

H3K36me1EErep1_FE.bw
H3K36me1EErep1_logLR.bw

2.5 峰注釋（Peak_anno）

ChIPseeker

ChIPseeker的功能分為三類：
注釋：提取peak附近最近的基因，注釋peak所在區域。
比較：估計ChIP peak數據集中重疊部分的顯著性；整合GEO數據集，以便于將當前結果和已知結果比較。
可視化：peak的覆蓋情況；TSS區域結合的peak的平均表達譜和熱圖；基因組注釋；TSS距離；peak和基因的重疊。

# 加載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.knownGene library("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")

Chip peaks coverage plot

查看peak在全基因組的位置

covplot(ring1B,chrs=c(“chr17”, “chr18”)) #specific chr

ring1B

suz12

cbx7

RYBP

ring1B(chr17&18)

● Heatmap of ChIP binding to TSS regions

promoter <- 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 method

plotAvgProf(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 plot

plotAnnoBar(peakAnno)
vennpie(peakAnno)
upsetplot(peakAnno)

餅圖:

條形圖：

upsetplot: upset技術適用于多于5個集合的表示情況。

可視化TSS區域的TF binding loci

plotDistToTSS(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 comparision

peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb,
tssRegion=c(-3000, 3000), verbose=FALSE)
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)

Overlap of peaks and annotated genes

genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)

總結

以上是生活随笔為你收集整理的ChIP-seq 数据分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇：常用地理信息数据下载平台
下一篇：工资软件测试白盒测试报告,白盒测试测试报