ChIP-seq笔记
文章目錄
- ChIP-seq學(xué)習(xí)
- 1 數(shù)據(jù)下載
- 1.1 數(shù)據(jù)主要分為三個(gè)部分
- 1.2 從NCBI上下載數(shù)據(jù)
- 2 質(zhì)量控制
- 2.1 軟件安裝
- 2.2 轉(zhuǎn)化數(shù)據(jù)格式 sra ----- fastq
- 2.3 下載小鼠參考基因組的index
- 2.4 質(zhì)量檢測
- 2.5 質(zhì)控結(jié)果批量查看
- 2.6 數(shù)據(jù)清理
- 3 對比到參考基因組
- 3.1 bowtie2
- 3.1.1 首先介紹一下bowtie2的一些參數(shù)
- 3.2 index必須與待匹配的read.fastq文件放在同一個(gè)目錄下，方可開始匹配。
- 4 搜峰
- 4.1 MACS2
- 4.2MACS2 核心: callpeak 用法
- 4.3 實(shí)際操作語句
- 4.4 參數(shù)介紹
- 4.5 callpeak 結(jié)果文件說明
- 4.5.1
- 5 峰注釋
- 6 作圖分析
- 6.1 ChIP peaks結(jié)合TSS 區(qū)域的情況
- 6.1.1 Heatmap of ChIP binding to TSS regions (chip結(jié)合到TSS區(qū)域的熱圖)
- 6.1.2 Average Profile of ChIP peaks binding to TSS region (chip峰值對TSS區(qū)域的平均配置)
- 6.2 peaks注釋
- 7 組合可視化分析
ChIP-seq學(xué)習(xí)
這是我第一次做ChIP-seq，將所有的步驟以及代碼全部記錄下來，如有錯(cuò)誤歡迎大家指正。
chip-seq主要有四個(gè)步驟
Cross-linking（DNA和蛋白質(zhì)交聯(lián)）
Sonication（超聲將染色體切割）
IP（利用抗原抗體的特異性識別）
Sequencing（測序）
（Linux操作系統(tǒng)CentOS）
流程圖
1 數(shù)據(jù)下載
1.1 數(shù)據(jù)主要分為三個(gè)部分
（1）ivf GSE112546
（2）scnt GSE112546
（3）zyy GSE73952
1.2 從NCBI上下載數(shù)據(jù)
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR714/SRR7145717/SRR7145717.sra #Morula.IVF.H3K27me3.Rep1 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR714/SRR7145718/SRR7145718.sra #Morula.IVF.H3K27me3.Rep2 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR714/SRR7145722/SRR7145722.sra #Morula.IVF.Inputftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR714/SRR7145719/SRR7145719.sra #Morula.SCNT.H3K27me3.Rep1 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR714/SRR7145720/SRR7145720.sra #Morula.SCNT.H3K27me3.Rep2 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR714/SRR7145721/SRR7145721.sra #Morula.SCNT.Inputftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208744/SRR3208744.sra #MII Oocyte Input ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208749/SRR3208749.sra #MII Oocyte H3K27me3 rep1 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208750/SRR3208750.sra #MII Oocyte H3K27me3 rep2 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208751/SRR3208751.sra #MII Oocyte H3K27me3 rep3ESC ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208788/SRR3208788.sra #ESC input ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208791/SRR3208791.sra #ESC H3K27me3 rep1 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208792/SRR3208792.sra #ESC H3K27me3 rep2 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208793/SRR3208793.sra #ESC H3K27me3 rep3 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208794/SRR3208794.sra #ESC H3K27me3 rep4TSC ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208795/SRR3208795.sra #TSC input ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208799/SRR3208799.sra #TSC H3K27me3 rep1 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208800/SRR3208800.sra #TSC H3K27me3 rep2 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR320/SRR3208801/SRR3208801.sra #TSC H3K27me3 rep32 質(zhì)量控制
2.1 軟件安裝
#安裝sratoolkit
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.2/sratoolkit.2.9.2-centos_linux64.tar.gz#解壓sratoolkit
tar -zxvf sratoolkit.2.9.2-centos_linux64.tar.gz2.2 轉(zhuǎn)化數(shù)據(jù)格式 sra ----- fastq
for i in *sra do echo $i /data/sunyu/h3k27/sratoolkit/sratoolkit.2.9.2-centos_linux64/bin/fastq-dump --split-3 $i; done2.3 下載小鼠參考基因組的index
wget -c "ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip" &``` unzip mm10.zip2.4 質(zhì)量檢測
for i in *fastq do fastqc -t 4 $i done2.5 質(zhì)控結(jié)果批量查看
multiqc *fastqc.zip --export如果沒有添加環(huán)境變量可采用一下方法
export PATH=$PATH:/data/sunyu/anaconda2/bin multiqc *fastqc.zip --export##trimmomatic
安裝 trimmomatic
wget -c http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip & unzip Trimmomatic-0.38.zip2.6 數(shù)據(jù)清理
-threads 設(shè)置多線程運(yùn)行
java -jar /data/sunyu/h3k27/trimmomatic/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 10 -phred33 SRR7145717.fastq output_SRR7145717.fastq ILLUMINACLIP:/data/sunyu/h3k27/trimmomatic/Trimmomatic-0.38/adapters/TruSeq3-SE.fa:2:30:10 TRAILING:4 HEADCROP:15 MINLEN:36在數(shù)據(jù)清理的過程中存在兩個(gè)模式
單末端 SE
和雙末端 PE
java -jar /path/Trimmomatic/trimmomatic-0.36.jar PE -phred33 -trimlog logfile reads_1.fq.gz reads_2.fq.gz out.read_1.fq.gz out.trim.read_1.fq.gz out.read_2.fq.gz out.trim.read_2.fq.gz ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50如果大家想要了解更多數(shù)據(jù)質(zhì)控給大家推薦個(gè)連接
https://zhuanlan.zhihu.com/p/28802083
| ILLUMINACLIP | 去接頭 |
| adapters | Exome.fa ：adapter 序列的 fasta 文件 |
| 2 | 16 個(gè)堿基長度的種子序列中可以有 2 個(gè)錯(cuò)配 |
| 30 | 采用回文模式時(shí)匹配得分至少為30 (約50個(gè)堿基) |
| 10 | 采用簡單模式時(shí)匹配得分至少為10 (約17個(gè)堿基) |
| ILLUMINACLIP | “$adapter”/Exome.fa:2:30:10 \ |
| LEADING:3 | 從序列的開頭開始去掉質(zhì)量值小于 3 的堿基； |
| TRAILING:3 | 從序列的末尾開始去掉質(zhì)量值小于 3 的堿基； |
| SLIDINGWINDOW:4:15 | 從 5’ 端開始以 4 bp 的窗口計(jì)算堿基平均質(zhì)量， |
| 如果此平均值低于 15 | 則從這個(gè)位置截?cái)?read； |
| HEADCROP | 在reads的首端切除指定的長度； |
| MINLEN:36 | 如果 reads 長度小于 36 bp 則扔掉整條 read。 |
3 對比到參考基因組
3.1 bowtie2
3.1.1 首先介紹一下bowtie2的一些參數(shù)
#bowtie2 [options]* -x {-1 -2 | -U } -S []
#-p/–threads NTHREADS 設(shè)置線程數(shù). Default: 1
#-q reads 是 fastq 格式的
#-x index 路徑
#-1 雙末端測序的 _1.fastq 路徑。可以為多個(gè)文件，并用逗號分開；多個(gè)文件必須和 -2 中制定的文件一一對應(yīng)。
#-2 雙末端測序的 _2.fastq 路徑.
#-U 非雙末端測序的 fastq 路徑。可以為多個(gè)文件，并用逗號分開。
#-S 輸出 Sam 格式文件。
#-3/–trim3 剪掉3’端長度的堿基，再用于比對。(default: 0).
#用fastqc看了看數(shù)據(jù)質(zhì)量，發(fā)現(xiàn)3端質(zhì)量有點(diǎn)問題，我就用了-3 5 --local參數(shù)，–local 如果fq文件是沒有經(jīng)過 trim 的，可以用局部比對執(zhí)行 soft-clipping，加上參數(shù)–local 。該模式下對read進(jìn)行局部比對, 從而, read 兩端的一些堿基不比對，從而使比對得分滿足要求.
3.2 index必須與待匹配的read.fastq文件放在同一個(gè)目錄下，方可開始匹配。
bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/ivf/SRR7145717.fastq -S /data/sunyu/h3k27/ivf/ivf_1/MII_rep1.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/ivf/SRR7145718.fastq -S /data/sunyu/h3k27/ivf/ivf_1/MII_rep2.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/ivf/SRR7145722.fastq -S /data/sunyu/h3k27/ivf/ivf_1/MII_rep3.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/scnt/SRR7145719.fastq -S /data/sunyu/h3k27/scnt/scnt_1/MII_rep4.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/scnt/SRR7145720.fastq -S /data/sunyu/h3k27/scnt/scnt_1/MII_rep5.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/scnt/SRR7145721.fastq -S /data/sunyu/h3k27/scnt/scnt_1/MII_rep6.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208744_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208744_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep7.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208749_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208749_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep8.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208750_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208750_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep9.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208751_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208751_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep10.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/zyy/SRR3208791.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep11.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/zyy/SRR3208792.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep12.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/zyy/SRR3208793.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep13.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -U /data/sunyu/h3k27/zyy/SRR3208794.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep14.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208795_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208795_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep15.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208799_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208799_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep16.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208800_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208880_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep17.sam bowtie2 -p 6 -3 5 --local -x /data/sunyu/h3k27/bowtie/bowtie2-2.3.4.3-linux-x86_64/index/mm10 -1 /data/sunyu/h3k27/zyy/SRR3208801_1.fastq -2 /data/sunyu/h3k27/zyy/SRR3208881_2.fastq -S /data/sunyu/h3k27/zyy/zyy_1/MII_rep18.sam
這里只展示三張圖片其余的模式和這三張截屏一樣。
在程序運(yùn)行時(shí)可采用腳本形式
4 搜峰
4.1 MACS2
peaks calling：尋找可能的結(jié)合位點(diǎn)，即基因組中大量reads富集的區(qū)域。
MACS2介紹 https://www.jianshu.com/p/e83a7e10ea2e
4.2MACS2 核心: callpeak 用法
Example for regular peak calling
macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01Example for broad peak calling
macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.14.3 實(shí)際操作語句
macs2 callpeak -t /data/sunyu/h3k27/ivf/ivf_1/MII_rep1.sam -c /data/sunyu/h3k27/ivf/ivf_1/MII_rep3.sam -f SAM -g mm -n ivf -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/ivf/ivf_1/MII_rep2.sam -c /data/sunyu/h3k27/ivf/ivf_1/MII_rep3.sam -f SAM -g mm -n ivf1 -B -q 0.01macs2 callpeak -t /data/sunyu/h3k27/scnt/scnt_1/MII_rep4.sam -c /data/sunyu/h3k27/scnt/scnt_1/MII_rep6.sam -f SAM -g mm -n scnt1 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/scnt/scnt_1/MII_rep5.sam -c /data/sunyu/h3k27/scnt/scnt_1/MII_rep6.sam -f SAM -g mm -n scnt2 -B -q 0.01macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep8.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep7.sam -f SAM -g mm -n zyy1 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep9.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep7.sam -f SAM -g mm -n zyy2 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep10.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep7.sam -f SAM -g mm -n zyy3 -B -q 0.01macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep16.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep15.sam -f SAM -g mm -n zyy4 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep17.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep15.sam -f SAM -g mm -n zyy5 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep18.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep15.sam -f SAM -g mm -n zyy6 -B -q 0.01macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep11.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep19.sam -f SAM -g mm -n zyy7 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep12.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep19.sam -f SAM -g mm -n zyy8 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep13.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep19.sam -f SAM -g mm -n zyy9 -B -q 0.01 macs2 callpeak -t /data/sunyu/h3k27/zyy/zyy_1/MII_rep14.sam -c /data/sunyu/h3k27/zyy/zyy_1/MII_rep19.sam -f SAM -g mm -n zyy10 -B -q 0.014.4 參數(shù)介紹
-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自動(dòng)檢測。如果不提供這項(xiàng)，就是自動(dòng)檢測選擇。
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——基因組大小，默認(rèn)提供了hs, mm, ce, dm選項(xiàng)，不在其中的話，比如說擬南芥，就需要自己提供了（擬南芥根據(jù)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——輸出文件的前綴名。表示實(shí)驗(yàn)的名字, 請取一個(gè)有意義的名字。
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閾值， 默認(rèn)是0.05。q值是根據(jù)p值利用BH計(jì)算，也就是多重試驗(yàn)矯正后的結(jié)果。
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有關(guān)，而MFOLD和MACS預(yù)構(gòu)建模型有關(guān)，默認(rèn)是5：50，MACS會先尋找100多個(gè)peak區(qū)構(gòu)建模型，一般不用改，因?yàn)槟悴欢?br /> 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.
4.5 callpeak 結(jié)果文件說明
callpeak會得到如下結(jié)果文件（先用ivf的結(jié)果舉例）：
4.5.1
ivf_summits.bed：Browser Extensible Data，記錄每個(gè)peak的peak summits，換句話說就是記錄極值點(diǎn)的位置。MACS建議用該文件尋找結(jié)合位點(diǎn)的motif。能夠直接載入U(xiǎn)CSC browser，用其他軟件分析時(shí)需要去掉第一行。bed包含有3個(gè)必須的字段和9個(gè)可選字段。
三個(gè)字段包括
| chrom | 染色體名字 |
| chromStart | 染色體起始位點(diǎn) |
| chromEND | 染色體終止位點(diǎn) |
這里的chromStart是起始于0，而不是1。
可選的9個(gè)字段包括：
| name | 名字 |
| score | 分值（0-1000），用于genome browser展示時(shí)上色。 |
| strand | 正負(fù)鏈，對于chip-seq數(shù)據(jù)來說，一般沒有正負(fù)鏈信息 |
| thickstart | 畫矩形圖的起點(diǎn) |
| thickEnd | 畫矩形的終點(diǎn) |
| itemRgd | RGD值 |
| blockCount | 子元件（比如外顯子）的數(shù)目 |
| blockSizes | 子元件的大小 |
| blockStarts | 子元件的起始位點(diǎn) |
一般情況下，我們只用到前面的5個(gè)字段，這也是做peak calling的MACS輸出的字段。
其中第5個(gè)字段，MACS的解釋是這個(gè)樣子的：
The 5th column in this file is the summit height of fragment pileup.
5 峰注釋
主要會用到R語言中的ChIPseeker包
ChIPseeker的功能分為三類：
注釋：提取peak附近最近的基因，注釋peak所在區(qū)域。
比較：估計(jì)ChIP peak數(shù)據(jù)集中重疊部分的顯著性；整合GEO數(shù)據(jù)集，以便于將當(dāng)前結(jié)果和已知結(jié)果比較。
可視化：peak的覆蓋情況；TSS區(qū)域結(jié)合的peak的平均表達(dá)譜和熱圖；基因組注釋；TSS距離；peak和基因的重疊。
安裝r語言包
source ("https://bioconductor.org/biocLite.R") biocLite("DOSE") biocLite("ChIPseeker") biocLite("org.Mm.eg.db") biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene") biocLite("clusterProfiler") biocLite("ReactomePA")簡單出個(gè)圖（數(shù)據(jù)zyy1查看峰值在全基因組上的位置）
library("ChIPseeker") library("org.Mm.eg.db") library("TxDb.Mmusculus.UCSC.mm10.knownGene") txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene library("clusterProfiler") setwd("E:/data/h3k27/zyy") suz12<-readPeakFile("zyy1_peaks.narrowPeak") covplot(suz12,weightCol=5)covplot函數(shù)可以計(jì)算peak 在染色體上的覆蓋區(qū)域，并可視化。
covplot(peak, weightCol = NULL, xlab = “Chromosome Size (bp)”, ylab = “”,
title = “ChIP Peaks over Chromosomes”, chrs = NULL, xlim = NULL,
lower = 1)
Arguments
| peak | 峰值文件 或者 GRanges 對象 |
| weightCol | weight column of peak（峰值柱的寬度） |
| xlab | xlab |
| ylab | ylab |
| title | title |
| chrs | 選擇的染色體在默認(rèn)情況下，所有的染色體 |
| xlim | 繪圖的范圍，默認(rèn)是整個(gè)染色體 |
| lower | lower cutoff of coverage signal（覆蓋信號的最下端節(jié)點(diǎn)） |
(先采用一個(gè)數(shù)據(jù)把所有的圖都做出來）
在安裝R語言包是真的是比較復(fù)雜的，安了好多次，有的是因?yàn)榫W(wǎng)絡(luò)的問題鏡像連接不上，有的是因?yàn)槠渌妮o助包的問題，希望大家在做的時(shí)候耐心一些。
這張圖更難看 哈哈 沒法展示
6 作圖分析
接下來將利用r語言來做出更加細(xì)致的圖
6.1 ChIP peaks結(jié)合TSS 區(qū)域的情況
TSS:transcription start site(轉(zhuǎn)錄起始位點(diǎn))
首先，計(jì)算ChIP peaks結(jié)合在TSS區(qū)域的情況。這就需要準(zhǔn)備TSS區(qū)域，這一般定義在TSS位點(diǎn)的側(cè)翼序列（默認(rèn)-3000~+3000）。然后比對 map到這些區(qū)域的peak，并生成tagMatrix。
getPromoters(TxDb = NULL, upstream = 1000, downstream = 1000, by = “gene”)
| TxDb | TxDb |
| upstream | upstream from TSS site （轉(zhuǎn)錄起始位點(diǎn)上游） |
| downstream | downstream from TSS site （轉(zhuǎn)錄起始位點(diǎn)下游） |
| by | one of gene or transcript （一個(gè)基因或者是轉(zhuǎn)錄組） |
getTagMatrix(peak, weightCol = NULL, windows)
peak 峰值文件 或者 GRanges 對象
weightCol column name of weight, default is NULL（柱形圖的列名，默認(rèn)值為空）
windows a collection of region with equal size, eg. promoter region.（一組大小相同的區(qū)域，例如。啟動(dòng)子區(qū)域。）
6.1.1 Heatmap of ChIP binding to TSS regions (chip結(jié)合到TSS區(qū)域的熱圖)
為了計(jì)算結(jié)合TSS區(qū)域的ChIP峰的分布，我們應(yīng)該準(zhǔn)備TSS區(qū)域，其被定義為TSS位點(diǎn)的側(cè)翼序列。然后對齊映射到這些區(qū)域的峰，并生成tagMatrix。
tagHeatmap(tagMatrix, xlim, xlab = “”, ylab = “”, title = NULL, color = “red”)
| tagMatrix | tagMatrix or a list of tagMatrix(標(biāo)記矩陣或標(biāo)記矩陣列表) |
| xlim | xlim |
| xlab | xlab |
| ylab | ylab |
| title | title |
| color | color |
6.1.2 Average Profile of ChIP peaks binding to TSS region (chip峰值對TSS區(qū)域的平均配置)
plotAvgProf2(peak, weightCol = NULL, TxDb = NULL, upstream = 1000, downstream = 1000, xlab = “Genomic Region (5’->3’)”, ylab = “Peak Count Frequency”, conf, facet = “none”, free_y = TRUE, verbose = TRUE, …)
| peak | peak file or GRanges object |
| weightCol | column name of weight |
| TxDb | TxDb object |
| upstream | upstream position |
| downstream | downstream position |
| xlab | xlab |
| ylab | ylab |
| conf | confidence interval |
| facet | one of ‘none’, ‘row’ and ‘column’ |
| free_y | if TRUE, y will be scaled by AvgProf |
| verbose | print message or not |
| … | additional parameter |
6.2 peaks注釋
annotatePeak函數(shù)進(jìn)行peaks注釋，可以定義TSS（轉(zhuǎn)錄起始位點(diǎn)）區(qū)域，默認(rèn)情況下TSS定義為-3kb到+ 3kb。
annotatePeak的輸出是csAnno格式。 ChIPseeker中國的as.GRanges函數(shù)將csAnno轉(zhuǎn)換為GRanges格式，as.data.frame將csAnno轉(zhuǎn)換為data.frame，然后通過write.table將其導(dǎo)出到文件。
TxDb.Hsapiens.UCSC.hg38.knownGene TxDb.Hsapiens.UCSC.hg19.knownGene分別對應(yīng)人類基因組hg38和hg19，TxDb.Mmusculus。
UCSC.mm10.knownGene和TxDb.Mmusculus.UCSC.mm9.knownGene則對應(yīng)小鼠基因組mm10和mm9。
用戶還可以通過RmakeTxDbFromBiomart和makeTxDbFromUCSC從UCSC Genome Bioinformatics和BioMart數(shù)據(jù)庫檢索準(zhǔn)備自己的TxDb對象。然后進(jìn)行峰值注釋。
所有的峰值信息都會保存在輸出文件中。其中包含peak最近的gene的位置和鏈的信息，從peak到最近的gene的TSS的距離等。鑒于某些信息可能的重疊，ChIPseeker采取以下優(yōu)先級別
annotatePeak(peak, tssRegion = c(-3000, 3000), TxDb = NULL,
level = “transcript”, assignGenomicAnnotation = TRUE,
genomicAnnotationPriority = c(“Promoter”, “5UTR”, “3UTR”, “Exon”, “Intron”,
“Downstream”, “Intergenic”), annoDb = NULL, addFlankGeneInfo = FALSE,
flankDistance = 5000, sameStrand = FALSE, ignoreOverlap = FALSE,
ignoreUpstream = FALSE, ignoreDownstream = FALSE, overlap = “TSS”,
verbose = TRUE)
| peak | peak file or GRanges object |
| tssRegion | Region Range of TSS |
| TxDb | TxDb object |
| level | one of transcript and gene |
| assignGenomicAnnotation | logical, assign peak genomic annotation or not |
| genomicAnnotationPriority | genomic annotation priority |
| annoDb | annotation package |
| addFlankGeneInfo | logical, add flanking gene information from the peaks |
| flankDistance | distance of flanking sequence |
| sameStrand | logical, whether find nearest/overlap gene in the same strand |
| ignoreOverlap | logical, whether ignore overlap of TSS with peak |
| ignoreUpstream | logical, if True only annotate gene at the 3’ of the peak. |
| ignoreDownstream | logical, if True only annotate gene at the 5’ of the peak. |
| overlap | one of ‘TSS’ or ‘a(chǎn)ll’, if overlap=“all”, then gene overlap with peak will be reported as nearest gene, no matter the overlap is at TSS region or not. |
| verbose | print message or not |
有以下幾種可視化方式(本方法根據(jù)簡書https://www.jianshu.com/p/2b8e2ea26665)
要根據(jù)基因組特征注釋給定峰的位置，annotatePeak在輸出的“注釋”列中為基因組注釋指定峰，其中包括峰是否在TSS，外顯子，5’UTR，3’UTR，Intronic或基因間。許多研究人員對這些注釋非常感興趣。TSS區(qū)域可以由用戶定義并annotatePeak輸出哪些基因的外顯子/內(nèi)含子的細(xì)節(jié)，如前一部分所示。
A：pie plot
啟動(dòng)
5’UTR
3’UTR
外顯子
內(nèi)含子
下游
基因間
B: bar plot
C: venn plot(重疊）
D: upsetplot（重疊）
6.3.2 Visualize distribution of TF-binding loci relative to TSS
peak（TF結(jié)合位點(diǎn)）到最近的gene的TSS之間的距離可以有annotatePeak函數(shù)進(jìn)行計(jì)算。作者提供了plotDistToTSS函數(shù)計(jì)算最近基因的TSS上游和下游的結(jié)合位點(diǎn)的百分比，并可視化這種分布。
7 組合可視化分析
上述的圖都是單個(gè)一個(gè)數(shù)據(jù)的圖，接下來我將ivf，scnt，ESC，TSC，oocyte細(xì)胞數(shù)據(jù)整合到一起做成venn圖。
這篇博客是在這各路大神的基礎(chǔ)上完成的，是我的一個(gè)學(xué)習(xí)筆記，歡迎大家參考借鑒。如有問題可以相互溝通交流。
總結(jié)
以上是生活随笔為你收集整理的ChIP-seq笔记的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： PASCAL VOC 2012 and
- 下一篇： 地理空间数据和大数据平台Spark结合能