ChIP-Seq分析小实战（三）

来源：Kai评论1,585

ChIP-Seq分析小实战（一）

ChIP-Seq分析小实战（二）

ChIP-Seq分析小实战（三）

在ChIP-Seq分析小实战（一）中以将数据做好了准备，接着就是按照文章中所述将ChIP-seq数据用botwie比对至NCBIM37基因组上，但是这次我还是按照教程的方法，使用bowtie2比对至mm10基因组上

The millions of reads produced by ChIP-seq of RYBP, Cbx7, Ring1B, Suz12,
and H2AK119Ub were aligned with the mouse genome (version NCBIM37)using the Bowtie (Langmead et al., 2009) tool version 0.12.7; two mismatches
were allowed within the seed alignment

Mapping reads

由于用fastqc质控后，发现Ring1B、Suz12和IgG_old样本的3’端有大约5个bp长度的碱基质量不是太好（大概在Q25左右），按照原来是惯例应该用软件将其剪去，但是看了教程发现可以用bowtie2中的参数在比对时进行处理，去官网手册上查了下

-3/--trim3 <int> Trim bases from 3’ (right) end of each read before alignment (default: 0)

--local In this mode, Bowtie 2 does not require that the entire read align from one end to the other. Rather, some characters may be omitted (“soft clipped”) from the ends in order to achieve the greatest possible alignment score

其实我们在用bowtie2时默认是使用--end-to-end，但是其是要求 entire read align from one end to the other, without any trimming (or “soft clipping”) of characters from either end，所以在使用-3/--trim3参数时，需要再加个--local参数

即比对代码如下，由于机子较差，因此将样本一个个运行，写了简单的shell脚本，同时将比对结果输出到log文件中，然后用samtools将sam文件转化成bam文件，并排序

#/bin/bash

GENOME=/home/anlan/reference/index/bowtie2/mm10/mm10

bowtie2 -p 5 --trim3 5 --local -x $GENOME -U SRR620204.fastq 2> ring1B_bowtie2.log |samtools sort -O bam -@ 5 -o ring1B.bam
bowtie2 -p 5 --trim3 5 --local -x $GENOME -U SRR620205.fastq 2> cbx7_bowtie2.log |samtools sort -O bam -@ 5 -o cbx7.bam
bowtie2 -p 5 --trim3 5 --local -x $GENOME -U SRR620206.fastq 2> suz12_bowtie2.log |samtools sort -O bam -@ 5 -o suz12.bam
bowtie2 -p 5 --trim3 5 --local -x $GENOME -U SRR620207.fastq 2> RYBP_bowtie2.log |samtools sort -O bam -@ 5 -o RYBP.bam
bowtie2 -p 5 --trim3 5 --local -x $GENOME -U SRR620208.fastq 2> IgGold_bowtie2.log |samtools sort -O bam -@ 5 -o IgGold.bam
bowtie2 -p 5 --trim3 5 --local -x $GENOME -U SRR620209.fastq 2> IgG_bowtie2.log |samtools sort -O bam -@ 5 -o IgG.bam

从比对结果上来看，比对率均并不高，如ring1B_bowtie2.log

19687027 reads; of these:
      19687027 (100.00%) were unpaired; of these:
        7960096 (40.43%) aligned 0 times
        8614711 (43.76%) aligned exactly 1 time
        3112220 (15.81%) aligned >1 times
59.57% overall alignment rate

最后别忘记对每个bam文件建索引，不然后续步骤会报错

samtools index ring1B.bam

Peak calling

Peak calling可以理解为寻找peak出现的位置，而这些位置可能是一些靶蛋白结合的位点，我们做ChIP-seq也就是为了研究捕获到的序列以及其差异。至于其peak calling的原理以及我接下来使用的MACS2软件的原理可以参照Chip-seq处理流程，本来想贴这篇文章作者的博客链接的，结果发现博客打不开了！看了这文章后，大概能了解一点其统计学的含义

软件的下载及安装，一看是python的软件，如果只是尝试使用，我一般选择bioconda安装（我已配置好了conda环境了）

conda install -c bioconda macs2

如果报错：CondaHTTPError: HTTP 000 CONNECTION FAILED for url <https://nanomirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda/linux-64/repodata.json\>

说明墙太高了，那么需要对conda的源文件进行点处理（前提你已经给canda添加了清华源后还报错）

vim ~/.condarc
删除 - defaults 的行

安装成功后就可以用macs2对每个bam文件进行call peaks了，简单看一下可能会用到的几个参数

  • -t/–treatment 输入文件，支持很多格式，搭配-f/–format使用
  • -f/–format 设定输入文件的格式，这里我会选择bam格式
  • -c/–control 对照组（或者模拟的）数据，需要跟-t的输出文件在相同目录下
  • -n/–name 输出文件（有很多个文件）的前缀
  • –outdir 输出文件所在目录（不设定的话就默认当前目录了）
  • -g/–gsize 提供基因组的大小，程序有默认的几个物种可以选hs,mm,ce,dm
  • -q/–qvalue 设定FDR的阈值，默认是0.05
  • -p/–pvalue 设定pvalue的阈值，如果参数设定了-p，则其会覆盖参数-q的作用
  • -B/–bdg If this flag is on, MACS will store the fragment pileup, control lambda, -log10pvalue and -log10qvalue scores in bedGraph files.

看到文章中有对macs2的软件参数的选择的

Peaks with a FDR lower than 5% from Cbx7, Ring1B, Suz12, and H2AK119ub,
and lower than 1% from RYBP, were combined to detect chromosomal
regions for further analyses.

We observed an overlap of RYBP peaks (3,918 in
total) with 14%, 42%, and 37% of Cbx7, Ring1B, and Suz12
peaks, respectively (false discovery rate [FDR] = 1%; p = 1 3
10-5) (Figure 1A

所以对每个样本的macs2代码如下（选择的对照组的bam文件为IgGold.bam，不太了解IgG.bam在文中的作用）:

macs2 callpeak -c IgGold.bam -t suz12.bam -p 1e-5 -f BAM -g mm -n suz12 2> suz12.macs2.log
macs2 callpeak -c IgGold.bam -t ring1B.bam -p 1e-5 -f BAM -g mm -n ring1B 2> ring1B.macs2.log
macs2 callpeak -c IgGold.bam -t cbx7.bam -p 1e-5 -f BAM -g mm -n cbx7 2> cbx7.macs2.log
macs2 callpeak -c IgGold.bam -t RYBP.bam -p 1e-5 -f BAM -g mm -n RYBP 2>RYBP.macs2.log

然后粗略看一下peak位置信息的文件，看看每个样本找到多少个peak（bed的格式，储存了每个peaks的位置信息）

wc -l *.bed

2198 cbx7_summits.bed
7072 ring1B_summits.bed
6872 RYBP_summits_2.bed
2 RYBP_summits.bed
6633 suz12_summits.bed

回到文章中一看，cbx7、ring1B，RYBP和suz12对应的peaks number分别对应2754、6982、3918和8054。。除了RYBP相差的也太严重了，其他的样本也有点差别。RYBP样本肯定无法进行后续分析了，所以按照教程中所说的，从https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42466把RYBP样本的peaks文件下载下来用用

wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42466/suppl/GSE42466_RYBP_peaks_5.txt.gz
gunzip GSE42466_RYBP_peaks_5.txt.gz
mv GSE42466_RYBP_peaks_5.txt RYBP_summits.bed

macs2 call peak出来各个文件的作用可以看官网介绍https://github.com/taoliu/MACS/

  • NAME_peaks.xls 包含called peaks的信息
  • NAME_peaks.narrowPeak is BED6+4 format file which contains the peak locations together with peak summit, pvalue and qvalue
  • NAME_summits.bed BED format, which contains the peak summits locations for every peaks.If you want to find the motifs at the binding sites, this file is recommended

接下来应该就是对peaks进行注释，然后寻找其靶基因以及结合位点等

发表评论

匿名网友