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

本次给学徒讲解的文章是 : 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 '{$h{$F[0]}++ if exists $h{$F[0]}; $h{$F[0]}=1 unless exists $h{$F[0]};print "$F[0]$h{$F[0]}\t$F[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 $id
arr=($id)
srr=${arr[1]}
sample=${arr[0]}
# 单端测序数据的sra转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测序数据如下：

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 $id
file=$(basename $id )
sample=${file%%.*} 
echo $sample

computeMatrix reference-point  --referencePoint TSS  -p 15  \
-b 10000 -a 10000    \
-R  $bed \
-S $id  \
--skipZeros  -o matrix1_${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 ${sample}_Heatmap_10K.png
plotHeatmap -m matrix1_${sample}_TSS_10K.gz  -out ${sample}_Heatmap_10K.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_${sample}_TSS_10K.gz  -out ${sample}_Profile_10K.png
plotProfile -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 $id
file=$(basename $id )
sample=${file%%.*} 
echo $sample

computeMatrix reference-point  --referencePoint TSS  -p 15  \
-b 2000 -a 2000    \
-R  $bed \
-S $id  \
--skipZeros  -o matrix1_${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 ${sample}_Heatmap_2K.png
plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.png
plotProfile -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 $id
file=$(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}_motifDir -len 8,10,12
annotatePeaks.pl    homer_peaks.tmp mm10  1>${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数据分析教程给写完了

赞 (0)

相关推荐