学员分享-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叔。。


推荐阅读更多精彩内容