学员分享-Chip-seq 实战分析流程

Chip-seq 实战分析流程

本次实践是根据Jimmy 和洲更两位大神的笔记，还有Y叔的Chipseeker包，第一次跑Chip-seq流程，主要关注了流程能否跑通，怎么解决自己实践中的Bug。

下面根据我的实践，先介绍下分析流程和跑流程时的注意事项。

---

Chip-seq 分析流程可以分为两个模块：

前期的准备工作, 包括软件的安装，数据的下载（或是自己测序的数据）；

Chip实验和数据获取：

第一步： 将蛋白交联到DNA上。 也就是保证蛋白和DNA能够结合，找到互作位点。

第二步： 通过超声波剪切DNA链。

第三步： 加上附上抗体的磁珠用于免疫沉淀靶蛋白。抗体很重要

第四步： 接触蛋白交联；纯化DNA

第五步：送去测序。

Chip-seq的分析流程：

质量控制, 用到的是FastQC

序列比对, Bowtie2或这BWA

peak calling，此处用的MACS

peak注释, 这里用的Y叔的ChIPseeker

---

软件的下载和环境的配置

所需要的软件有sratoolkit, fastqc,bowtie2,samtools,macs2,htseq-count,bedtools,由于我用的是所里的大服务器，软件都已经安装，我只需要配置自己的环境变量。具体的安装代码可以在生信技能树公众号后台回复老司机获得，或者参考转录组实战分析中 浙大植物学小白的转录组笔记。

数据的下载

下载样本数据

1. for ((i=204;i<=209;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620$i/SRR620$i.sra;done

用sratookit 将.sra 文件转化为fastq文件

1. ls *sra |while read id; do /software/biosoft/software/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump --split-3 $id;done

NOTE:

这里要注意是双端测序还是单端测序。

• fastq-dump 转换sra文件成fastq/fasta 文件

pair-end: fastq-dump --split-3 *.sra

single-end: fastq-dump *.sra 或者 fastq-dump --fasta *sra

下载小鼠参考基因组的索引和注释文件

1. wget ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip

2. unzip mm10.zip

也可以自己下载基因组mm10,然后建索引：

● 下载参考基因组

1. wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz

2. tar -zxvf chromFa.tar.gz

3. cat *.fa mm10.fa

rm chr*.fa

● build index

1. mkdir -p ~/reference/index/bowtie

2. cd ~/reference/index/bowtie

3. nohup time bowtie2-build  ~/reference/mm10.fa  ~/reference/index/bowtie/mm10 1mm10.bowtie_index.log 2&1 &

下载注释文件

https://genome.ucsc.edu/cgi-bin/hgTables

---

序列比对

将得到的fastq文件用bowtie2比对小鼠参考基因组上

1. bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620204.fastq | samtools sort -O bam -o analysis/alignment/ring1B.bam

2. bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620205.fastq | samtools sort -O bam -o analysis/alignment/cbx7.bam

3. bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620206.fastq | samtools sort -O bam -o analysis/alignment/suz12.bam

4. bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620207.fastq | samtools sort -O bam -o analysis/alignment/RYBP.bam

5. bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620208.fastq | samtools sort -O bam -o analysis/alignment/IgGold.bam

6. bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620209.fastq | samtools sort -O bam -o analysis/alignment/IgG.bam

PS:本文的作者不会写shell脚本做循环，大家不要学习这个。

结果：

计算比对率 （samtools flagstat）

这里我是通过samtools flagstat 计算的比对率，结果如下：

1. ring1B : 60.49%

2. cdx7: 87.52%

3. suz12: 67.01%

4. RYBP: 84.17%

5. IgGold: 57.80%

6. IgG: 82.80%

用IGV查看

从官网下载IGV, 解压即可使用，linux下 igv.sh打开IGV界面，Windows 下点击igv.bat.

● 首先载入参考基因组，可以载入自己下载好的参考基因组，也可选择IGV中含有的参考基因组，ref_Genome 必须是fasta格式。

● 载入比对的文件，比对的文件必须先经过sort 和 index, 才可加载。

1. samtools sort a.bam a.sort

2. samtools index a.sort.bam

1.  igv.sh

● 比对可视化结果： IgG是对照组，其他组有的峰在对照组没有，即为peak.

---

用MACS call peak

Macs call peak

1. macs2 callpeak -c IgGold.bam -t suz12.bam -q 0.05 -f BAM -g mm -n suz12 2suz12.macs2.log

2. macs2 callpeak -c IgGold.bam -t ring1B.bam -q 0.05 -f BAM -g mm -n ring1B 2/ring1B.macs2.log

3. macs2 callpeak -c IgGold.bam -t cbx7.bam -q 0.05 -f BAM -g mm -n cbx7 2cbx7.macs2.log

4. macs2 callpeak -c IgGold.bam -t RYBP.bam -q 0.01 -f BAM -g mm -n RYBP 2RYBP.macs2.log

结果文件不只是上述提到的一类，还有如下格式：

● NAMEpeaks.xls: 以表格形式存放peak信息，虽然后缀是xls，但其实能用文本编辑器打开，和bed格式类似，但是以1为基，而bed文件是以0为基.也就是说xls的坐标都要减一才是bed文件的坐标 ● NAMEpeaks.narrowPeak NAMEpeaks.broadPeak 类似。后面4列表示为， integer score for display， fold-change，-log10pvalue，-log10qvalue，relative summit position to peak start。内容和NAMEpeaks.xls基本一致，适合用于导入R进行分析。 ● NAMEsummits.bed：记录每个peak的peak summits，话句话说就是记录极值点的位置。MACS建议用该文件寻找结合位点的motif。 ● NAMEmodel.r，能通过$ Rscript NAME_model.r作图，得到是基于你提供数据的peak模型

计算peak数

RYBP无数据，估计是数据上传错误，可以下载作者上传的peak数据。 ● 下载RYBP的peak数据

1. wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42466/suppl/GSE42466_RYBP_peaks_5.txt.gz

2. gzip -d GSE42466_RYBP_peaks_5.txt.gz

3. mv GSE42466_RYBP_peaks_5.txt RYBP2_summits.bed

---

结果注释与可视化

结果的注释用的是Y叔 的 Chipseeker包。

ChIPseeker的功能分为三类: ● 注释：提取peak附近最近的基因， 注释peak所在区域 ● 比较：估计ChIP peak数据集中重叠部分的显著性；整合GEO数据集，以便于将当前结果和已知结果比较 ● 可视化： peak的覆盖情况；TSS区域结合的peak的平均表达谱和热图；基因组注释；TSS距离；peak和基因的重叠。

下载chipseeker 有关包

• download packages

1. source ("https://bioconductor.org/biocLite.R")

2. biocLite("ChIPseeker")

3. biocLite("org.Mm.eg.db")

4. biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")

5. biocLite("clusterProfiler")

6. biocLite("ReactomePA")

7. biocLite("DOSE")

● loading packages

1. library("ChIPseeker")

2. library("org.Mm.eg.db")

3. library("TxDb.Mmusculus.UCSC.mm10.knownGene")

4. txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene

5. library("clusterProfiler")

读入bed文件

1. ring1B <- readPeakFile("F:/Chip-seq_exercise/ring1B_peaks.narrowPeak")

Chip peaks coverage plot

查看peak在全基因组的位置

1. covplot(ring1B,chrs=c("chr17", "chr18"))   #specific chr

ring1B

suz12

cbx7

RYBP

ring1B(chr17&18)

● Heatmap of ChIP binding to TSS regions

1. promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)

2. tagMatrix <- getTagMatrix(ring1B, windows=promoter)

3. tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

Average Profile of ChIP peaks binding to TSS region

● Confidence interval estimated by bootstrap method

1. plotAvgProf(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.

1. peakAnno <- annotatePeak(ring1B, tssRegion=c(-3000, 3000),

2. TxDb=txdb, annoDb="org.Mm.eg.db")

可视化 Pie and Bar plot

1. plotAnnoBar(peakAnno)

2. vennpie(peakAnno)

3. upsetplot(peakAnno)

饼图:

条形图：

upsetplot: upset技术适用于多于5个集合的表示情况。

可视化TSS区域的TF binding loci

1. plotDistToTSS(peakAnno,

2. title="Distribution of transcription factor-binding loci\nrelative to TSS")

多个peak的比较

多个peak set注释时，先构建list,然后用lapply.list(name1=bed_file1,name2=bed_file2) RYBP的数据有问题，这里加上去，会一直报错。

1. peaks <- list(cbx7=cbx7,ring1B=ring1B,suz12=suz12)

2. promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)

3. tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)

4. plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))

5. plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf=0.95,resample=500, facet="row")

6. tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=NULL)

ChIP peak annotation comparision

1. peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb,

2. tssRegion=c(-3000, 3000), verbose=FALSE)

3. plotAnnoBar(peakAnnoList)

4. plotDistToTSS(peakAnnoList)

Overlap of peaks and annotated genes

1. genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)

2. vennplot(genes)

---

• 这篇主要介绍了步骤

• 了解基础知识

peak calling的统计学原理 http://www.plob.org/2014/05/08/7227.html

根据比对的bam文件来对peaks区域可视化 http://www.bio-info-trainee.com/1843.html

wig、bigWig和bedgraph文件详解 http://www.bio-info-trainee.com/1815.html

• 文章综述

ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. https://www.ncbi.nlm.nih.gov/pubmed/22955991

---

初次完成ChIP-seq的实践过程，虽然学习资料和代码都很全，但是整个实践过程并不顺利，遇到了很多共性和个性的问题，在问题的解决过程中学到了很多知识。感谢整个实践过程中提供过解惑的健明师兄，徐州更同学，还有Y叔。。