當(dāng)前位置：首頁 >

ChIP-seq数据处理流程(附赠长达5小时的视频指导)

發(fā)布時(shí)間：2023/12/14 29 豆豆

生活随笔收集整理的這篇文章主要介紹了 ChIP-seq数据处理流程(附赠长达5小时的视频指导) 小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

本次給學(xué)徒講解的文章是 : 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).

查看文章發(fā)現(xiàn)數(shù)據(jù)是： 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個(gè)樣本。

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

在SRA數(shù)據(jù)庫可以下載：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 這里需要警覺了，參考基因組應(yīng)該是鼠。

我這里： 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文件，要想利用這個(gè)數(shù)據(jù)，只能自己重新處理，這就是為什么需要學(xué)會(huì)ChIP-seq數(shù)據(jù)處理的原因。不過作者給了bw文件，所以可以勉強(qiáng)跟自己的結(jié)果相互驗(yàn)證。

這里作者使用的是 Illumina Genome Analyzer II 測序儀，有點(diǎn)過時(shí)了，測序策略是 se50。

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

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先進(jìn)行檢索

conda search trim_galore

保證所有的軟件都是安裝在 wes 這個(gè)環(huán)境下面

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(‘a(chǎn)irway’,‘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叔的包檢查會(huì)有版本的問題，包括 ChIPseeker

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(ChIPpeakAnno)
library(ChIPseeker)
下載sra并且轉(zhuǎn)換為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命令創(chuàng)建 srr.list 文件。

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

默認(rèn)下載目錄：~/ncbi/public/sra/

ls -lh ~/ncbi/public/sra/

下載耗時(shí)，自行解決，學(xué)員使用現(xiàn)成數(shù)據(jù)：/public/project/epi/Chipseq-OS25_Esc/OS25_Esc/sra

假如提前下載好了數(shù)據(jù)。

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]}?\t$ F[1]”}’ > config
得到內(nèi)容如下：

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轉(zhuǎn)fq文件：

下面需要用循環(huán)

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 $i d a r r = ($ id)
srr= ${arr[1]} sample=$ {arr[0]}

單端測序數(shù)據(jù)的sra轉(zhuǎn)fasq

nohup $dump -A $sample -O $analysis_dir --gzip --split-3 sra/$ 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測序數(shù)據(jù)如下：

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軟件進(jìn)行質(zhì)控
這個(gè)時(shí)候選擇trim_galore軟件進(jìn)行過濾，單端測序數(shù)據(jù)的代碼如下；

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應(yīng)該走一波

cd ~/project/epi/qc

相對目錄需要理解

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

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

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

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

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 $i d) s a m p l e =$ {file%%.}
echo $file $sample

比對過程3分鐘一個(gè)樣本

$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文件進(jìn)行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文件
因?yàn)橐粋€(gè)樣品分成了多個(gè)lane進(jìn)行測序，所以在進(jìn)行peaks calling的時(shí)候，需要把bam進(jìn)行合并。

如果不用循環(huán)：

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個(gè)fq測序數(shù)據(jù)只剩下8個(gè)樣本啦。（我下載的時(shí)候漏掉了2個(gè)sra文件，也就是漏掉了一個(gè)樣本。）

假如需要去除PCR重復(fù)
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重復(fù)前后比較：

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進(jìn)行找peaks
macs2包含一系列的子命令，其中最主要的就是callpeak，官方提供了使用實(shí)例

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

-t: 實(shí)驗(yàn)組的輸出結(jié)果

-c: 對照組的輸出結(jié)果

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

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

-n: 輸出文件的前綴名

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

-q: q值，也就是最小的PDR閾值，默認(rèn)是0.05。q值是根據(jù)p值利用BH計(jì)算，也就是多重試驗(yàn)矯正后的結(jié)果。

-p：這個(gè)是p值，指定p值后MACS2就不會(huì)用q值了。

-m: 和MFOLD有關(guān)，而MFOLD和MACS預(yù)構(gòu)建模型有關(guān)，默認(rèn)是5：50，MACS會(huì)先尋找100多個(gè)peak區(qū)構(gòu)建模型，一般不用改，因?yàn)槟愫艽蟾怕噬喜粫?huì)懂。

所以我這里給學(xué)徒講解的實(shí)戰(zhàn)代碼是：

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
其實(shí)上面的-B 參數(shù)意義也不大，得到的bedgraph文件沒啥用。

得到的bed格式的peaks文件的行數(shù)如下：

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
因?yàn)镸ockIP是control，所以它自己跟自己比較，肯定是沒有peaks的。

值得注意的是S5P并不是一個(gè)樣本多個(gè)lane，而是本身樣本有重復(fù)，其實(shí)是需要分開的。

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

而且很有趣的是我前幾個(gè)月處理這個(gè)數(shù)據(jù)集的時(shí)候使用的過濾低質(zhì)量reads參數(shù)是短于 35bp的全部丟棄，現(xiàn)在是短于25bp的全部拋棄，導(dǎo)致了得到的peaks從數(shù)量上千差別不小。

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

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

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進(jìn)行可視化。這里的參數(shù)僅使用了-e/–extendReads和-bs/–binSize即拓展了原來的read長度，且設(shè)置分箱的大小。其他參數(shù)還有

–filterRNAstrand {forward, reverse}: 僅統(tǒng)計(jì)指定正鏈或負(fù)鏈

–region/-r CHR:START:END: 選取某個(gè)區(qū)域統(tǒng)計(jì)

–smoothLength: 通過使用分箱附近的read對分箱進(jìn)行平滑化

如果為了其他結(jié)果進(jìn)行比較，還需要進(jìn)行標(biāo)準(zhǔn)化，deeptools提供了如下參數(shù)：

–scaleFactor: 縮放系數(shù)

–normalizeUsingRPKMReads: Per Kilobase per Million mapped reads (RPKM)標(biāo)準(zhǔn)化

–normalizeTo1x: 按照1x測序深度(reads per genome coverage, RPGC)進(jìn)行標(biāo)準(zhǔn)化

–ignoreForNormalization：指定那些染色體不需要經(jīng)過標(biāo)準(zhǔn)化

如果需要以100為分箱，并且標(biāo)準(zhǔn)化到1x，且僅統(tǒng)計(jì)某一條染色體區(qū)域的正鏈，輸出格式為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類似，只不過需要提供兩個(gè)樣本，并且采用SES方法進(jìn)行標(biāo)準(zhǔn)化，于是多了–ratio參數(shù)。

首先把bam文件轉(zhuǎn)為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附件信號(hào)強(qiáng)度：

首先對單一樣本繪圖：

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

如果要批處理，需要學(xué)習(xí)好linux命令。

首先畫10K附近

bed=/public/annotation/CHIPseq/mm10/ucsc.refseq.bed
for id in /home/jmzeng/project/epi/mergeBam/bw ;
do
echo $i d f i l e =$ (basename $i d) s a m p l e =$ {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: \dotstmap -m matrix1_?$ {sample}_TSS_10K.gz -out $KaTeX parse error: Expected group after '_' at position 82: \dotsfile -m matrix1_?$ {sample}_TSS_10K.gz -out $KaTeX parse error: Expected group after '_' at position 48: \dotsfile -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 $i d f i l e =$ (basename $i d) s a m p l e =$ {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: \dotstmap -m matrix1_?$ {sample}_TSS_2K.gz -out $KaTeX parse error: Expected group after '_' at position 81: \dotsfile -m matrix1_?$ {sample}_TSS_2K.gz -out $KaTeX parse error: Expected group after '_' at position 47: \dotsfile -m matrix1_?$ {sample}_TSS_2K.gz -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720

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

還可以給多個(gè)bed文件來繪圖，還可以畫genebody的圖，因?yàn)樵硪粯?#xff0c;我就不做過多介紹啦。

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

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

使用R包對找到的peaks文件進(jìn)行注釋
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來出自動(dòng)化報(bào)告。

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

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

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

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

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文件

數(shù)據(jù)庫下載取決于網(wǎng)速咯

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

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

或者用下面代碼安裝:

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整合了兩個(gè)方法，包括依賴于數(shù)據(jù)庫的查詢，和de novo的推斷,都是讀取ChIP-seq數(shù)據(jù)上游分析得到的bed格式的peaks文件。

運(yùn)行homer軟件

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

cd ~/project/epi/motif
for id in /home/jmzeng/project/epi/peaks/.bed;
do
echo $i d f i l e =$ (basename $i d) s a m p l e =$ {file%%.}
echo $sample
awk ‘{print $4"\t"$1"\t"$2"\t"$3"\t+"}’ $id >homer_peaks.tmp
findMotifsGenome.pl homer_peaks.tmp mm10 ${sample}_motifDir -len 8,10,12 annotatePeaks.pl homer_peaks.tmp mm10 1>$ {sample}.peakAnn.xls 2>${sample}.annLog.txt
done
把上面的代碼保存為腳本runMotif.sh，然后運(yùn)行：nohup bash runMotif.sh 1>motif.log &

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

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

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

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

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

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

結(jié)果忘記了還有視頻這回事，趕緊錄制了，同樣的B站免費(fèi)送給大家看，但是不提供售后，網(wǎng)易云課堂也有視頻，點(diǎn)擊下面的閱讀原文直達(dá)！

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

上一篇： android虚拟摇杆
下一篇：《大数据之路：阿里巴巴大数据实践》-第1

ChIP-seq数据处理流程(附赠长达5小时的视频指导)

可以用search先進(jìn)行檢索

保證所有的軟件都是安裝在 wes 這個(gè)環(huán)境下面

R

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

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

vim 或者cat命令創(chuàng)建 srr.list 文件。

默認(rèn)下載目錄：~/ncbi/public/sra/

下載耗時(shí)，自行解決，學(xué)員使用現(xiàn)成數(shù)據(jù)：/public/project/epi/Chipseq-OS25_Esc/OS25_Esc/sra

假如提前下載好了數(shù)據(jù)。

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

下面需要用循環(huán)

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

單端測序數(shù)據(jù)的sra轉(zhuǎn)fasq

相對目錄需要理解

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

相對目錄需要理解

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

比對過程3分鐘一個(gè)樣本

如果不用循環(huán)：

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

通常我們用批處理。

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

首先對單一樣本繪圖：

both -R and -S can accept multiple files

both plotHeatmap and plotProfile will use the output from computeMatrix

如果要批處理 ，需要學(xué)習(xí)好linux命令。

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

both plotHeatmap and plotProfile will use the output from computeMatrix

both plotHeatmap and plotProfile will use the output from computeMatrix

loading packages

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

數(shù)據(jù)庫下載取決于網(wǎng)速咯

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

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

或者用下面代碼安裝:

總結(jié)

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

如果要批處理，需要學(xué)習(xí)好linux命令。