本次給學徒講解的文章是 : Brookes, E. et al. Polycomb associates genome-wide with a specific RNA polymerase II variant, and regulates metabolic genes in ESCs. Cell Stem Cell 10, 157–170 (2012).

查看文章發現數據是： Polycomb associates genome-wide with a specific RNA polymerase II variant, and regulates metabolic genes in ES cells (ChIP-Seq) https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE34518 總共是9個樣本。

但是很多樣本都分開在多個lane測序的，所以每個樣本其實是有多個sra文件，多個fastq文件。

在SRA數據庫可以下載 ：https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP009883 包括：Examination of 4 different RNAPII modifications (S5p, S7p, 8WG16, S2p), and the histone modifications H2Aub1 and H3K36me3 in mouse ES cells 這里需要警覺了，參考基因組應該是鼠。

我這里 ： Use prefetch to download them all, then transform those SRA files to fastq files by **sra-toolkits **, then align them to mm10, and call peaks.

作者并沒有給peaks文件，要想利用這個數據，只能自己重新處理，這就是為什么需要學會ChIP-seq數據處理的原因。不過作者給了bw文件，所以可以勉強跟自己的結果相互驗證。

這里作者使用的是 Illumina Genome Analyzer II 測序儀，有點過時了，測序策略是 se50。

從文章找到數據的ID： https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP009883 把下面的內容保存到文件，命名為 srr.list 就可以使用prefetch這個函數來下載。

SRR391032
SRR391033
SRR391034
SRR391035
SRR391036
SRR391037
SRR391038
SRR391039
SRR391040
SRR391041
SRR391042
SRR391043
SRR391044
SRR391045
SRR391046
SRR391047
SRR391048
SRR391049
SRR391050
安裝必備軟件：
#!/bin/bash

wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

conda create -n chipseq python=2 bwa
conda info --envs
source activate chipseq

可以用search先進行檢索

conda search trim_galore

保證所有的軟件都是安裝在 wes 這個環境下面

conda install -y sra-tools
conda install -y trim-galore samtools
conda install -y deeptools homer meme
conda install -y macs2 bowtie bowtie2

R

還需要安裝必備R包：

options(“repos” = c(CRAN=“https://mirrors.tuna.tsinghua.edu.cn/CRAN/”))
install.packages(“devtools”,
repos=“https://mirrors.tuna.tsinghua.edu.cn/CRAN/”)
library(devtools)
source(“https://bioconductor.org/biocLite.R”)
options(BioC_mirror=“https://mirrors.ustc.edu.cn/bioc/”)
BiocInstaller::biocLite(c(‘airway’,‘DESeq2’,‘edgeR’,‘limma’))
BiocInstaller::biocLite(c(‘ChIPpeakAnno’,‘ChIPseeker’))
BiocInstaller::biocLite(‘TxDb.Hsapiens.UCSC.hg19.knownGene’,
ask=F,suppressUpdates=T)
BiocInstaller::biocLite(‘TxDb.Hsapiens.UCSC.hg38.knownGene’,
ask=F,suppressUpdates=T)
BiocInstaller::biocLite(‘TxDb.Mmusculus.UCSC.mm10.knownGene’,
ask=F,suppressUpdates=T)

值得注意的是Y叔的包檢查會有版本的問題，包括 ChIPseeker

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(ChIPpeakAnno)
library(ChIPseeker)
下載sra并且轉換為fastq
使用下面的代碼，指定好自己的 prefetch 軟件命令即可。

prefetch=/home/jianmingzeng/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/prefetch
source activate chipseq
prefetch=prefetch

cat srr.list |while read id;do (nohup $prefetch $id -X 100G & );done

mkdir -p ~/project/epi/
cd ~/project/epi/
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd sra

vim 或者cat命令創建 srr.list 文件。

cat srr.list |while read id;do ( nohup $prefetch $id & );done

默認下載目錄：~/ncbi/public/sra/

ls -lh ~/ncbi/public/sra/

下載耗時，自行解決，學員使用現成數據：/public/project/epi/Chipseq-OS25_Esc/OS25_Esc/sra

假如提前下載好了數據。

cd ~/project/epi/
ln -s /public/project/epi/Chipseq-OS25_Esc/OS25_Esc/sra sra
第一步需要制作配置文件，代碼是：

直接用excel制作config文件，或者寫代碼

cut -f 4,7 sra.table |cut -d":" -f 2 |sed ‘s/ChIPSeq//g’ | sed ‘s/MockIP//g’|sed ‘s/^ //’ |tr ’ ’ ‘_’ |perl -alne ‘{KaTeX parse error: Expected '}', got 'EOF' at end of input: h{F[0]}++ if exists KaTeX parse error: Expected '}', got 'EOF' at end of input: h{F[0]}; KaTeX parse error: Expected '}', got 'EOF' at end of input: h{F[0]}=1 unless exists KaTeX parse error: Expected '}', got 'EOF' at end of input: h{F[0]};print “$F[0]$h{KaTeX parse error: Expected 'EOF', got '}' at position 5: F[0]}?\tF[1]”}’ > config
得到內容如下：

RNAPII_S5P_1 SRR391032
RNAPII_S5P_2 SRR391033
RNAPII_S2P_1 SRR391034
RNAPII_S7P_1 SRR391035
RNAPII_8WG16_1 SRR391036
RNAPII_8WG16_2 SRR391037
RNAPII_S2P_2 SRR391038
RNAPII_S2P_3 SRR391039
RNAPII_S7P_2 SRR391040
H2Aub1_1 SRR391041
H2Aub1_2 SRR391042
H3K36me3_1 SRR391043
H3K36me3_2 SRR391044
Control_1 SRR391045
Control_2 SRR391046
Ring1B_1 SRR391047
Ring1B_2 SRR391048
Ring1B_3 SRR391049
RNAPII_S5PRepeat_1 SRR391050
有了上面的配置文件就可以批量sra轉fq文件：

下面需要用循環

cd ~/project/epi/
source activate chipseq
dump=’/home/jianmingzeng/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump’
dump=fastq-dump
analysis_dir=raw

下面用到的 config 文件，就是上面自行制作的。

cat config|while read id;
do echo $idarr=($id)
srr=$arr[1]sample=${arr[0]}

單端測序數據的sra轉fasq

nohup $dump -A $sample -O $analysi{s}_{d}ir??gzip??split?3sra/$srr.sra &
done
下載的sra文件如下：

-rw-rw-r-- 1 jianmingzeng jianmingzeng 474M Mar 23 14:29 SRR391032.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 473M Mar 23 14:30 SRR391033.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 406M Mar 23 14:30 SRR391034.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 346M Mar 23 14:31 SRR391035.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 324M Mar 23 14:31 SRR391036.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 393M Mar 23 14:32 SRR391037.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 191M Mar 23 14:32 SRR391038.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 304M Mar 23 14:34 SRR391039.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 234M Mar 23 14:34 SRR391040.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 322M Mar 23 14:35 SRR391041.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 548M Mar 23 14:36 SRR391042.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 597M Mar 23 14:37 SRR391043.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 537M Mar 23 14:37 SRR391044.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 255M Mar 23 14:38 SRR391045.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 342M Mar 23 14:38 SRR391046.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 159M Mar 23 14:39 SRR391047.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 164M Mar 23 14:39 SRR391048.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 438M Mar 23 14:40 SRR391049.sra
-rw-rw-r-- 1 jianmingzeng jianmingzeng 165M Mar 23 14:40 SRR391050.sra
得到的fastq測序數據如下：

427M Jul 28 15:25 Control_1.fastq.gz
527M Jul 28 15:26 Control_2.fastq.gz
507M Jul 28 15:26 H2Aub1_1.fastq.gz
843M Jul 28 15:28 H2Aub1_2.fastq.gz
881M Jul 28 15:28 H3K36me3_1.fastq.gz
858M Jul 28 15:29 RNAPII_S2P_1.fastq.gz
326M Jul 28 15:25 RNAPII_S2P_2.fastq.gz
489M Jul 28 15:26 RNAPII_S2P_3.fastq.gz
283M Jul 28 15:25 RNAPII_S5PRepeat_1.fastq.gz
745M Jul 28 15:27 RNAPII_S5P_2.fastq.gz
533M Jul 28 15:26 RNAPII_S7P_1.fastq.gz
393M Jul 28 15:25 RNAPII_S7P_2.fastq.gz
266M Jul 28 15:25 Ring1B_1.fastq.gz
274M Jul 28 15:25 Ring1B_2.fastq.gz

使用trim_galore軟件進行質控
這個時候選擇trim_galore軟件進行過濾，單端測序數據的代碼如下；

cd ~/project/epi/clean
analysis_dir=/home/jmzeng/project/epi
bin_trim_galore=“trim_galore”
ls …/raw/*gz | while read fq1;
do
nohup $bin_trim_galore -q 25 --phred33 --length 25 -e 0.1 --stringency 4 -o $analysis_dir/clean $fq1 &
done
過濾后的文件如下：

356M Jul 28 15:46 Control_1_trimmed.fq.gz
438M Jul 28 15:46 Control_2_trimmed.fq.gz
456M Jul 28 15:46 H2Aub1_1_trimmed.fq.gz
693M Jul 28 15:47 H2Aub1_2_trimmed.fq.gz
773M Jul 28 15:47 H3K36me3_1_trimmed.fq.gz
824M Jul 28 15:48 RNAPII_S2P_1_trimmed.fq.gz
282M Jul 28 15:45 RNAPII_S2P_2_trimmed.fq.gz
458M Jul 28 15:46 RNAPII_S2P_3_trimmed.fq.gz
210M Jul 28 15:45 RNAPII_S5PRepeat_1_trimmed.fq.gz
626M Jul 28 15:47 RNAPII_S5P_2_trimmed.fq.gz
405M Jul 28 15:46 RNAPII_S7P_1_trimmed.fq.gz
306M Jul 28 15:46 RNAPII_S7P_2_trimmed.fq.gz
218M Jul 28 15:45 Ring1B_1_trimmed.fq.gz
218M Jul 28 15:45 Ring1B_2_trimmed.fq.gz
很明顯，QC應該走一波

cd ~/project/epi/qc

相對目錄需要理解

ls …/raw/*gz|xargs fastqc -t 10 -o ./
ls …/clean/*gz|xargs fastqc -t 10 -o ./
使用bowtie2進行比對
然后直接用bowtie2進行比對和統計比對率, 需要提前下載參考基因組然后使用命令構建索引，或者直接就下載索引文件：

下載小鼠參考基因組的索引和注釋文件, 這里用常用的mm10

索引大小為3.2GB， 不建議自己下載基因組構建，可以直接下載索引文件，代碼如下：

mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
單端測序數據的比對代碼如下：

cd ~/project/epi/align

相對目錄需要理解

bin_bowtie2=’/home/jianmingzeng/biosoft/bowtie/bowtie2-2.2.9/bowtie2’
bin_bowtie2=bowtie2
bowtie2_index="/home/jianmingzeng/reference/index/bowtie/mm10"
bowtie2_index=/public/reference/index/bowtie/mm10

一定要搞清楚自己的bowtie2軟件安裝在哪里，以及自己的索引文件在什么地方！！！

ls …/clean/gz |while read id;
do
file=$(basename $id)sample=${file%%.}
echo $file $sample

比對過程3分鐘一個樣本

$bin_bowtie2 -p 5 -x $bowtie2_index -U $id | samtools sort -O bam -@ 5 -o - > ${sample}.bam
done
得到的文件如下：

374M Jul 28 16:50 Control_1_trimmed.bam
469M Jul 28 16:50 Control_2_trimmed.bam
502M Jul 28 16:50 H2Aub1_1_trimmed.bam
767M Jul 28 16:50 H2Aub1_2_trimmed.bam
834M Jul 28 16:50 H3K36me3_1_trimmed.bam
731M Jul 28 16:28 RNAPII_S2P_1_trimmed.bam
302M Jul 28 16:29 RNAPII_S2P_2_trimmed.bam
483M Jul 28 16:32 RNAPII_S2P_3_trimmed.bam
218M Jul 28 16:33 RNAPII_S5PRepeat_1_trimmed.bam
609M Jul 28 16:36 RNAPII_S5P_2_trimmed.bam
416M Jul 28 16:38 RNAPII_S7P_1_trimmed.bam
309M Jul 28 16:39 RNAPII_S7P_2_trimmed.bam
238M Jul 28 16:40 Ring1B_1_trimmed.bam
239M Jul 28 16:41 Ring1B_2_trimmed.bam
對bam文件進行QC

cd ~/project/epi/align
ls *.bam |xargs -i samtools index {}
ls *.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id “.bam”).stat & );done
比對成功率都挺好的：

Control_1_trimmed.stat:7438540 + 0 mapped (88.03% : N/A)
Control_2_trimmed.stat:7221781 + 0 mapped (86.40% : N/A)
H2Aub1_1_trimmed.stat:8969578 + 0 mapped (97.40% : N/A)
H2Aub1_2_trimmed.stat:13229916 + 0 mapped (97.53% : N/A)
H3K36me3_1_trimmed.stat:11737310 + 0 mapped (98.89% : N/A)
Ring1B_1_trimmed.stat:4634240 + 0 mapped (93.59% : N/A)
Ring1B_2_trimmed.stat:4646919 + 0 mapped (93.85% : N/A)
RNAPII_S2P_1_trimmed.stat:25018794 + 0 mapped (97.26% : N/A)
RNAPII_S2P_2_trimmed.stat:6112834 + 0 mapped (95.00% : N/A)
RNAPII_S2P_3_trimmed.stat:8675514 + 0 mapped (96.99% : N/A)
RNAPII_S5P_2_trimmed.stat:12182274 + 0 mapped (98.17% : N/A)
RNAPII_S5PRepeat_1_trimmed.stat:4163763 + 0 mapped (82.81% : N/A)
RNAPII_S7P_1_trimmed.stat:6386269 + 0 mapped (80.90% : N/A)
RNAPII_S7P_2_trimmed.stat:5971178 + 0 mapped (82.66% : N/A)
合并bam文件
因為一個樣品分成了多個lane進行測序，所以在進行peaks calling的時候，需要把bam進行合并。

如果不用循環：

samtools merge control.merge.bam Control_1_trimmed.bam Control_2_trimmed.bam

通常我們用批處理。

cd ~/project/epi/
mkdir mergeBam
source activate chipseq
cd ~/project/epi/align
ls .bam|sed ‘s/_[0-9]_trimmed.bam//g’ |sort -u |while read id;do samtools merge …/mergeBam/$id.merge.bam $id.bam ;done
得到全新的bam文件如下：

847M Jul 28 17:05 Control.merge.bam
1.3G Jul 28 17:06 H2Aub1.merge.bam
834M Jul 28 17:06 H3K36me3.merge.bam
1.5G Jul 28 17:08 RNAPII_S2P.merge.bam
831M Jul 28 17:09 RNAPII_S5P.merge.bam
218M Jul 28 17:09 RNAPII_S5PRepeat.merge.bam
722M Jul 28 17:09 RNAPII_S7P.merge.bam
472M Jul 28 17:10 Ring1B.merge.bam
14個fq測序數據只剩下8個樣本啦。（我下載的時候漏掉了2個sra文件，也就是漏掉了一個樣本。）

假如需要去除PCR重復
cd ~/project/epi/mergeBam
source activate chipseq
ls *merge.bam | while read id ;do (nohup samtools markdup -r $id $(basename $id “.bam”).rmdup.bam & );done
ls *.rmdup.bam |xargs -i samtools index {}
ls *.rmdup.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id “.bam”).stat & );done
去除PCR重復前后比較：

847M Jul 28 17:05 Control.merge.bam
753M Jul 28 17:11 Control.merge.rmdup.bam
1.3G Jul 28 17:06 H2Aub1.merge.bam
1.1G Jul 28 17:12 H2Aub1.merge.rmdup.bam
834M Jul 28 17:06 H3K36me3.merge.bam
793M Jul 28 17:11 H3K36me3.merge.rmdup.bam
1.5G Jul 28 17:08 RNAPII_S2P.merge.bam
1.2G Jul 28 17:12 RNAPII_S2P.merge.rmdup.bam
831M Jul 28 17:09 RNAPII_S5P.merge.bam
568M Jul 28 17:11 RNAPII_S5P.merge.rmdup.bam
218M Jul 28 17:09 RNAPII_S5PRepeat.merge.bam
212M Jul 28 17:11 RNAPII_S5PRepeat.merge.rmdup.bam
722M Jul 28 17:09 RNAPII_S7P.merge.bam
618M Jul 28 17:11 RNAPII_S7P.merge.rmdup.bam
472M Jul 28 17:10 Ring1B.merge.bam
427M Jul 28 17:11 Ring1B.merge.rmdup.bam
使用macs2進行找peaks
macs2包含一系列的子命令，其中最主要的就是callpeak， 官方提供了使用實例

macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
一般而言，我們照葫蘆畫瓢，按照這個實例替換對應部分就行了，介紹一下各個參數的意義

-t: 實驗組的輸出結果

-c: 對照組的輸出結果

-f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一個。如果不提供這項，就是自動檢測選擇。

-g: 基因組大小， 默認提供了hs, mm, ce, dm選項， 不在其中的話，比如說擬南芥，就需要自己提供了。

-n: 輸出文件的前綴名

-B: 會保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores

-q: q值，也就是最小的PDR閾值， 默認是0.05。q值是根據p值利用BH計算，也就是多重試驗矯正后的結果。

-p： 這個是p值，指定p值后MACS2就不會用q值了。

-m: 和MFOLD有關，而MFOLD和MACS預構建模型有關，默認是5：50，MACS會先尋找100多個peak區構建模型，一般不用改，因為你很大概率上不會懂。

所以我這里給學徒講解的實戰代碼是：

cd ~/project/epi/mergeBam
source activate chipseq
ls *merge.bam |cut -d"." -f 1 |while read id;
do
if [ ! -s ${id}_summits.bed ];
then
echo $id
nohup macs2 callpeak -c Control.merge.bam -t $id.merge.bam -f BAM -B -g mm -n $id --outdir …/peaks 2> $id.log &
fi
done

mkdir dup
mv rmdup dup/
cd dup/

ls *.merge.rmdup.bam |cut -d"." -f 1 |while read id;
do
if [ ! -s ${id}_rmdup_summits.bed ];
then
echo $id
nohup macs2 callpeak -c Control.merge.rmdup.bam -t $id.merge.rmdup.bam -f BAM -B -g mm -n ${id}_rmdup --outdir …/peaks 2> $id.log &
fi
done
其實上面的-B 參數意義也不大，得到的bedgraph文件沒啥用。

得到的bed格式的peaks文件的行數如下：

0 Control_summits.bed 1102 H2Aub1_summits.bed

89739 H3K36me3_summits.bed
27705 Ring1B_summits.bed
20043 RNAPII_S2P_summits.bed
38643 RNAPII_S5PRepeat_summits.bed
61805 RNAPII_S5P_summits.bed
72498 RNAPII_S7P_summits.bed

0 Control_rmdup_summits.bed 1102 H2Aub1_rmdup_summits.bed

89739 H3K36me3_rmdup_summits.bed
27705 Ring1B_rmdup_summits.bed
20043 RNAPII_S2P_rmdup_summits.bed
38643 RNAPII_S5PRepeat_rmdup_summits.bed
61805 RNAPII_S5P_rmdup_summits.bed
72326 RNAPII_S7P_rmdup_summits.bed
因為MockIP是control，所以它自己跟自己比較，肯定是沒有peaks的。

值得注意的是S5P并不是一個樣本多個lane，而是本身樣本有重復，其實是需要分開的。

而且可以看到是否去除PCR重復，對找到的peaks數量沒有影響。

而且很有趣的是我前幾個月處理這個數據集的時候使用的過濾低質量reads參數是短于 35bp的全部丟棄，現在是短于25bp的全部拋棄，導致了得到的peaks從數量上千差別不小。

使用deeptool是進行可視化
下面的文字摘抄自生信技能樹論壇：https://vip.biotrainee.com/d/226 不過代碼純粹是我自己手打。

deeptools提供bamCoverage和bamCompare進行格式轉換，為了能夠比較不同的樣本，需要對先將基因組分成等寬分箱(bin)，統計每個分箱的read數，最后得到描述性統計值。對于兩個樣本，描述性統計值可以是兩個樣本的比率，或是比率的log2值，或者是差值。如果是單個樣本，可以用SES方法進行標準化。

bamCoverage的基本用法

source activate chipseq
bamCoverage -e 170 -bs 10 -b ap2_chip_rep1_2_sorted.bam -o ap2_chip_rep1_2.bw

ap2_chip_rep1_2_sorted.bam是前期比對得到的BAM文件

得到的bw文件就可以送去IGV/Jbrowse進行可視化。 這里的參數僅使用了-e/–extendReads和-bs/–binSize即拓展了原來的read長度，且設置分箱的大小。其他參數還有

–filterRNAstrand {forward, reverse}: 僅統計指定正鏈或負鏈

–region/-r CHR:START:END: 選取某個區域統計

–smoothLength: 通過使用分箱附近的read對分箱進行平滑化

如果為了其他結果進行比較，還需要進行標準化，deeptools提供了如下參數：

–scaleFactor: 縮放系數

–normalizeUsingRPKMReads: Per Kilobase per Million mapped reads (RPKM)標準化

–normalizeTo1x: 按照1x測序深度(reads per genome coverage, RPGC)進行標準化

–ignoreForNormalization： 指定那些染色體不需要經過標準化

如果需要以100為分箱，并且標準化到1x，且僅統計某一條染色體區域的正鏈，輸出格式為bedgraph,那么命令行可以這樣寫

bamCoverage -e 170 -bs 100 -of bedgraph -r Chr4:12985884:12997458 --normalizeTo1x 100000000 -b 02-read-alignment/ap2_chip_rep1_1_sorted.bam -o chip.bedgraph
bamCompare和bamCoverage類似，只不過需要提供兩個樣本，并且采用SES方法進行標準化，于是多了–ratio參數。

首先把bam文件轉為bw文件，詳情：http://www.bio-info-trainee.com/1815.html

cd ~/project/epi/mergeBam
source activate chipseq
ls *.bam |xargs -i samtools index {}
ls .bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.}.bw &
done

cd dup
ls *.bam |xargs -i samtools index {}
ls .bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.}.rm.bw &
done
查看TSS附件信號強度：

首先對單一樣本繪圖：

both -R and -S can accept multiple files

mkdir -p ~/project/epi/tss
cd ~/project/epi/tss
computeMatrix reference-point --referencePoint TSS -p 15
-b 10000 -a 10000
-R /public/annotation/CHIPseq/mm10/ucsc.refseq.bed
-S /home/jmzeng/project/epi/mergeBam/H2Aub1.bw
–skipZeros -o matrix1_test_TSS.gz
–outFileSortedRegions regions1_test_genes.bed

both plotHeatmap and plotProfile will use the output from computeMatrix

plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.png
plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_test_TSS.gz -out test_Profile.png
plotProfile -m matrix1_test_TSS.gz -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720

如果要批處理 ，需要學習好linux命令。

首先畫10K附近

bed=/public/annotation/CHIPseq/mm10/ucsc.refseq.bed
for id in /home/jmzeng/project/epi/mergeBam/bw ;
do
echo $idfile=$(basename $id)sample=${file%%.}
echo $sample

computeMatrix reference-point --referencePoint TSS -p 15
-b 10000 -a 10000
-R $bed
-S KaTeX parse error: Undefined control sequence: \ at position 5: id \? ?--skipZeros -o…{sample}TSS_10K.gz
–outFileSortedRegions regions1${sample}_TSS_10K.bed

輸出的gz為文件用于plotHeatmap, plotProfile

both plotHeatmap and plotProfile will use the output from computeMatrix

plotHeatmap -m matrix1_${sample}_TSS_10K.gz -out KaTeX parse error: Expected group after '_' at position 48: …tmap -m matrix1_?{sample}_TSS_10K.gz -out KaTeX parse error: Expected group after '_' at position 82: …file -m matrix1_?{sample}_TSS_10K.gz -out KaTeX parse error: Expected group after '_' at position 48: …file -m matrix1_?{sample}_TSS_10K.gz -out ${sample}_Profile_10K.pdf --plotFileFormat pdf --perGroup --dpi 720

done
使用命令批量提交：nohup bash 10k.sh 1>10k.log &

然后畫2K的

bed=/public/annotation/CHIPseq/mm10/ucsc.refseq.bed
for id in /home/jmzeng/project/epi/mergeBam/bw ;
do
echo $idfile=$(basename $id)sample=${file%%.}
echo $sample

computeMatrix reference-point --referencePoint TSS -p 15
-b 2000 -a 2000
-R $bed
-S KaTeX parse error: Undefined control sequence: \ at position 5: id \? ?--skipZeros -o…{sample}TSS_2K.gz
–outFileSortedRegions regions1${sample}_TSS_2K.bed

both plotHeatmap and plotProfile will use the output from computeMatrix

plotHeatmap -m matrix1_${sample}_TSS_2K.gz -out KaTeX parse error: Expected group after '_' at position 47: …tmap -m matrix1_?{sample}_TSS_2K.gz -out KaTeX parse error: Expected group after '_' at position 81: …file -m matrix1_?{sample}_TSS_2K.gz -out KaTeX parse error: Expected group after '_' at position 47: …file -m matrix1_?{sample}_TSS_2K.gz -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720

done
使用命令批量提交：nohup bash 2k.sh 1>2k.log &

還可以給多個bed文件來繪圖，還可以畫genebody的圖，因為原理一樣，我就不做過多介紹啦。

上面的批量代碼其實就是為了統計全基因組范圍的peak在基因特征的分布情況，也就是需要用到computeMatrix計算，用plotHeatmap以熱圖的方式對覆蓋進行可視化，用plotProfile以折線圖的方式展示覆蓋情況。

computeMatrix具有兩個模式:scale-region和reference-point。前者用來信號在一個區域內分布，后者查看信號相對于某一個點的分布情況。無論是那個模式，都有有兩個參數是必須的，-S是提供bigwig文件，-R是提供基因的注釋信息。還有更多個性化的可視化選項。

使用R包對找到的peaks文件進行注釋
bedPeaksFile = ‘8WG16_summits.bed’;
bedPeaksFile

loading packages

require(ChIPseeker)
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
require(clusterProfiler)
peak <- readPeakFile( bedPeaksFile )
keepChr= !grepl(’_’,seqlevels(peak))
seqlevels(peak, pruning.mode=“coarse”) <- seqlevels(peak)[keepChr]
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000),
TxDb=txdb, annoDb=“org.Mm.eg.db”)
peakAnno_df <- as.data.frame(peakAnno)
可以載入IGV看看效果，檢測軟件找到的peaks是否真的合理，還可以配合rmarkdown來出自動化報告。

也可以使用其它代碼進行下游分析； https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq

peaks相關基因集的注釋
都是得到感興趣基因集，然后注釋，分析方法等同于GEO數據挖掘課程或者轉錄組下游分析： https://github.com/jmzeng1314/GEO （有配套視頻，就不多說了這里）

homer軟件來尋找motif
這個軟件安裝當初特別麻煩： https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step8-Homer-findMotif.sh

但是現在有了conda，一句話搞定：conda install -c bioconda homer , 找到自己安裝的homer，然后使用其附帶的配置腳本來下載數據庫咯。

perl ~/miniconda3/envs/chipseq/share/homer-4.9.1-5/configureHomer.pl -install mm10
ls -lh ~/miniconda3/envs/chipseq/share/homer-4.9.1-5/data/

我們上游分析是基于mm10找到的peaks文件

數據庫下載取決于網速咯

下載成功后會多出 ~/miniconda3/envs/chipseq/share/homer-4.9.1-5/data/genomes/mm9/ 文件夾, 共 4.9G

這個文件夾取決于你把homer這個軟件安裝到了什么地方。

或者用下面代碼安裝:

cd ~/biosoft
mkdir homer && cd homer
wget http://homer.salk.edu/homer/configureHomer.pl
perl configureHomer.pl -install
perl configureHomer.pl -install hg19
homer軟件找motif整合了兩個方法，包括依賴于數據庫的查詢，和de novo的推斷,都是讀取ChIP-seq數據上游分析得到的bed格式的peaks文件。

運行homer軟件

但是使用起來很簡單：http://homer.ucsd.edu/homer/ngs/peakMotifs.html

cd ~/project/epi/motif
for id in /home/jmzeng/project/epi/peaks/.bed;
do
echo $idfile=$(basename $id)sample=${file%%.}
echo $sample
awk ‘{print $4"\t"$1"\t"$2"\t"$3"\t+"}’ $id >homer_peaks.tmp
findMotifsGenome.pl homer_peaks.tmp mm10 ${sample}_{m}otifDir?len8,10,12annotatePeaks.plhome{r}_{p}eaks.tmpmm101>${sample}.peakAnn.xls 2>${sample}.annLog.txt
done
把上面的代碼保存為腳本runMotif.sh，然后運行：nohup bash runMotif.sh 1>motif.log &

不僅僅找了motif，還順便把peaks注釋了一下。得到的后綴為peakAnn.xls 的文件就可以看到和使用R包注釋的結果是差不多的。

還可以使用meme來找motif，需要通過bed格式的peaks的坐標來獲取fasta序列。MEME，鏈接：http://meme-suite.org/

其它高級分析
比如可以 比較不同的peaks文件，代碼見：https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step6-ChIPpeakAnno-Venn.R

當然了，本教程講解的是單端測序數據的處理，如果是雙端測序，里面的很多參數是需要修改的。

不過，只要你完整的看完了我前面的流程，掌握了linux和R，以及必備的基礎生物信息學知識，我相信你肯定能hold住雙端測序數據的學習啦。

本來以為我把ChIP-seq教程寫完了： 一不小心就把ChIP-seq數據分析教程給寫完了

結果忘記了還有視頻這回事，趕緊錄制了，同樣的B站免費送給大家看，但是不提供售后，網易云課堂也有視頻，點擊下面的閱讀原文直達！

總結

以上是生活随笔為你收集整理的ChIP-seq数据处理流程(附赠长达5小时的视频指导)的全部內容，希望文章能夠幫你解決所遇到的問題。

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