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

Posted on by

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

数据的下载

下载样本数据

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

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

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

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

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

● 下载参考基因组

wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz 
tar -zxvf chromFa.tar.gz 
cat *.fa mm10.fa

rm chr*.fa

● build index 

mkdir -p ~/reference/index/bowtie 
cd ~/reference/index/bowtie 
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比对小鼠参考基因组上

bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620204.fastq | samtools sort -O bam -o analysis/alignment/ring1B.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620205.fastq | samtools sort -O bam -o analysis/alignment/cbx7.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620206.fastq | samtools sort -O bam -o analysis/alignment/suz12.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620207.fastq | samtools sort -O bam -o analysis/alignment/RYBP.bam
bowtie2 -p 6 -3 5 --local -x reference/mm10 -U SRR620208.fastq | samtools sort -O bam -o analysis/alignment/IgGold.bam
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 计算的比对率，结果如下：

ring1B : 60.49%
cdx7: 87.52%
suz12: 67.01%
RYBP: 84.17%
IgGold: 57.80%
IgG: 82.80%

用IGV查看

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

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

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

samtools sort a.bam a.sort 
samtools index a.sort.bam
  igv.sh

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

用MACS call peak

Macs call peak

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

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

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

计算peak数

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

wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42466/suppl/GSE42466_RYBP_peaks_5.txt.gz
gzip -d GSE42466_RYBP_peaks_5.txt.gz
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
source ("https://bioconductor.org/biocLite.R")
biocLite("ChIPseeker")
biocLite("org.Mm.eg.db")
biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")
biocLite("clusterProfiler")
biocLite("ReactomePA")
biocLite("DOSE")

● loading packages

library("ChIPseeker")
library("org.Mm.eg.db")
library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
library("clusterProfiler")

读入bed文件

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

Chip peaks coverage plot

查看peak在全基因组的位置

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

suz12

cbx7

RYBP

ring1B(chr17&18)

● Heatmap of ChIP binding to TSS regions

promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(ring1B, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

Average Profile of ChIP peaks binding to TSS region

● Confidence interval estimated by bootstrap method

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.

peakAnno <- annotatePeak(ring1B, tssRegion=c(-3000, 3000),
TxDb=txdb, annoDb="org.Mm.eg.db")

可视化 Pie and Bar plot 

plotAnnoBar(peakAnno)
vennpie(peakAnno)
upsetplot(peakAnno)

饼图:

条形图：

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

可视化TSS区域的TF binding loci

plotDistToTSS(peakAnno,
title="Distribution of transcription factor-binding loci\nrelative to TSS")

多个peak的比较

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

peaks <- list(cbx7=cbx7,ring1B=ring1B,suz12=suz12)
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrixList <- lapply(peaks, getTagMatrix, windows=promoter)
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf=0.95,resample=500, facet="row")
tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=NULL)

ChIP peak annotation comparision

peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb,
tssRegion=c(-3000, 3000), verbose=FALSE)
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)

Overlap of peaks and annotated genes
genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)

一篇文章学会ChIP-seq分析（上）

一篇文章学会ChIP-seq分析（下）

  • 这篇主要介绍了步骤

ChIP-seq实战分析

  • 了解基础知识

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叔。。

Comments are closed.