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
结果文件不只是上述提到的一类，还有如下格式：
● 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数据
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)
原文中此处为链接，暂不支持采集
原文中此处为链接，暂不支持采集
-
这篇主要介绍了步骤
原文中此处为链接，暂不支持采集
-
了解基础知识
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叔。。
网友评论