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]
## 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]
## 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)
