搜索
查看: 4038|回复: 0

[CHIP-seq] ChIP-seq基础入门学习

[复制链接]

29

主题

131

帖子

1200

积分

金牌会员

Rank: 6Rank: 6

积分
1200
发表于 2017-12-17 22:03:51 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-12-17 22:12 编辑

ChIP-seq基础入门学习
首先要先感谢生信技能树Jimmy分享ChIP-seq基础入门教程,让我有机会以及动力去好好了解ChIP-seq,教程的提纲以这篇2013年发表在Cell Reports的文章RYBP and Cbx7 define specific biological functions of polycomb complexes in mouse embryonic stem cells为例,并给出了这篇文章ChIP-seq分析的每个步骤,可供我们实践。
由于文章部分代码过多,在这就不整理了,可以看有道笔记里的,排版整齐点
ChIP-seq基础入门学习

文章研究对象
本文主要讲了PRC1(Polycomb repressive complex 1)在小鼠的胚胎干细胞中有两类亚型Cbx7-PRC1和RYBP-PRC1,但是Cbx7和RYBP这两种是不能共存的,尽管两者的全基因组定位在某些gene上交叉,也就是说在基因组上结合在了同一个位置。在分子水平上,Cbx7用于招募Ring1B结合到染色质上,然而RYBP可以增强PRC1的酶活性。RYBP结合的基因,其有着了较低水平的Ring1B和H2AK119ub,但其相比结合Cbx7有较高的表达量。在功能上,RYBP结合的基因主要涉及代谢调节和细胞周期,Cbx7结合的基因主要涉及early-lineage commitment of ESCs。
在做出上述总结前,肯定需要对PRC1进行介绍,其组成看起来非常复杂:
The protein families that constitute the core of PRC1 contain several members:
Cbx (Cbx2, Cbx4, Cbx6, Cbx7, or Cbx8); Ring1A or Ring1B;
PHC (PHC1, PHC2, or PHC3); PCGF (PCGF1, PCGF2, PCGF3,
PCGF4, PCGF5, or PCGF6); and RYBP or YAF2.
接着作者抛出自己这篇文章想要说明的问题:
  • What is the genome-wide localization of these two types of
    PCR1 complexes?  
    • Is the expression of specific sets of genes differentially regulated by both complexes?
    • Do they exert common and/or unique biological functions?
    • Do they show any interdependency for localizing to chromatin?

接着就是对各个结论进行表述及说明,这个在后续分析时再说(主要我也说不清楚。。。)
我们需要先拿到这篇文章所用的ChIP-seq数据,如下:
Sequencing data have been deposited into the NCBI Gene Expression
Omnibus database under accession number GSE42466
所以我们通过https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42466网页下方的ftp进行下载界面
ChIP-seq数据有SRR620204到SRR620209总共6个文件,每个文件的对应是样本Ring1B、Cbx7、Suz12、RYBP、IgG_old和IgG。bx7、Ring1B、RYBP是复合体的组成,外加PRC2的Suz12以及2个对照组

数据的下载和整理
以SRR620204.sra为例:
ascp -T -i /home/anlan/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp-private.ncbi.nlm.nih.gov:sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620204/SRR620204.sra
转化为fq格式,检查质控fastqc
for i in $(seq 4 9);do fastq-dump --split-3 SRR62020${i}.sra ;done fastqc SRR62020*.fastq multiqc SRR62020*
由于使用bowtie2软件将reads比对至参考基因组上,所以要下载Bowtie2的mm10的索引
wget ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
最后需要mm10 refseq注释bed文件,进行网址下载http://genome.ucsc.edu/cgi-bin/hgTables,其中track和table选择RefSeq,然后输出格式为bed即可,最后下载(可以先在自己电脑上下载再传到服务器上,也可以像下面直接在服务器上下载)
curl 'http://genome.ucsc.edu/cgi-bin/hgTables?hgsid=646311755_P0RcOBvAQnWZSzQz2fQfBiPPSBen&boolshad.hgta_printCustomTrackHeaders=0&hgta_ctName=tb_ncbiRefSeq&hgta_ctDesc=table+browser+query+on+ncbiRefSeq&hgta_ctVis=pack&hgta_ctUrl=&fbQual=whole&fbUpBases=200&fbExonBases=0&fbIntronBases=0&fbDownBases=200&hgta_doGetBed=get+BED' >mm10.refseq.bed
Mapping reads
将数据准备好,接着就是按照文章中所述将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
由于用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 =
    10 -5)
    所以对每个样本的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/s ... 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进行注释,然后寻找其靶基因以及结合位点等


peak 注释
peak注释个人简单的理解为寻找chip-seq所对应的靶基因,先看下文章如何定义靶基因的:
Genes with a peak within the gene body, or within 2.5 kb from the TSS, were considered to be target genes
也就是说在基因区间内或距离TSS 2.5kb的peak的区域被认为是靶基因
我们先看文章中的Figure 1A:
其表示Cbx7,Ring1B,RYBP和Suz12之间peaks的overlapping的情况,那么我们就按照文章中所说的,用bedtools看看BYBP对其他3个样本peaks的重叠情况
ChIP-seq profiles around the
TSSs were generated for each IP by calculating the average coverage in
each position normalized for the total number of mapped reads with the
BED tools package
由于自己之前跑出来的peaks文件的在软件参数以及参考基因选择上跟文献会有不小的差别,因此下面统一用GSE网站上下载下来的peaks文件做后续分析
wc -l *peaks.bed 2754 Cbx7_peaks.bed 6982 Ring1b_peaks.bed 6872 RYBP_peaks.bed 8054 Suz12_peaks.bed
然后用bedtools对RYBP和Cbx7样本的peaks取overlap(才发现RYBP_peaks.bed的peaks数目跟文献的竟然不一样。。。),所以取Cbx7_peaks.bed和Ring1b_peaks.bed为例:
bedtools intersect -a Ring1b_peaks.bed -b Cbx7_peaks.bed >Ring1b_cbx7_coverage.bed
结果为2647个overlapping,跟文献的2463还是有点差距的,按照我现在的知识只能理解为工具上的使用的差异了。接着按照文献所示的图Figure 1B,找两个peaks文件的靶基因的交集(即venn图)
这里准备使用ChIPseeker包来根据peaks对靶基因进行注释。PS.这里我将使用自己MACS2跑出来的xxx_peaks.narrowPeak;不使用作者上传的是peaks文件是因为ChIPseeker包不识别。。。并且还需要将xxx_peaks.narrowPeak文件改为bed后缀,不然ChIPseeker包读进去会将第一行作为行名(我暂时还没在包说明文档中找到为啥会出现这种情况),所以我3个文件为下述:
cbx7_peaks.narrowPeak.bed
ring1B_peaks.narrowPeak.bed
suz12_peaks.narrowPeak.bed
首先读入必须的包和上述3个文件
library(ChIPseeker) library(org.Mm.eg.db) library(TxDb.Mmusculus.UCSC.mm10.knownGene) library(VennDiagram) cbx7 <- readPeakFile("cbx7_peaks.narrowPeak.bed") ring1B <- readPeakFile("ring1B_peaks.narrowPeak.bed") suz12 <- readPeakFile("suz12_peaks.narrowPeak.bed")   
接着我们需要用peakAnno函数对其进行peak注释,这里的注释主要可以分为三种,genomic注释、nearest gene注释以及peak上下游注释,具体解释可参照包作者写的文档peak注释。按照文章中对靶基因的要求,我们选择前两者注释即可
cbx7_peakAnno <- annotatePeak(cbx7, tssRegion = c(-2500, 2500), TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb="org.Mm.eg.db") ring1B_peakAnno <- annotatePeak(ring1B, tssRegion = c(-2500, 2500), TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb="org.Mm.eg.db") suz12_peakAnno <- annotatePeak(suz12, tssRegion = c(-2500, 2500), TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb="org.Mm.eg.db")   
然后提取出cbx7和ring1B的注释到的gene,绘制venn图
venn_list <- list(cbx7=as.data.frame(cbx7_peakAnno)$geneId,               ring1B=as.data.frame(ring1B_peakAnno)$geneId) venn.diagram(venn_list, imagetype = "png", fill = c("blue", "green"),           alpha = c(0.5, 0.5)  ,filename = "venn.png")
虽然结果和文章的有较大的差距,但是考虑到软件参数以及参考基因组的不同,也就算了。。

peak 注释可视化
Peak注释可视化即对注释到的peaks进行可视化展示,除了上面的靶基因的venn图,文章还对ChIP-seq profiles进行了展示
An area of 5 kb
surrounding each TSS that was associated to one or more ChIP-seq (RYBP,
Ring1B, Cbx7, Suz12, or H2AK119ub) was used to calculate the ChIP-seq
profile and the whole ChIP-seq coverage.
  • 使用ChIPseeker包可视化
      可以用covplot看一下ChIP-seq peaks在染色体上的分布图,比如cbx7的
      covplot(cbx7, weightCol="V5")
      做cbx7、ring1B和suz12样本在共有的靶基因TSS附近的profile图(代码续上述的靶基因注释)
      #先获得共有靶基因的  common_gene <- Reduce(intersect, list(cbx7=as.data.frame(cbx7_peakAnno)$geneId,                                ring1B=as.data.frame(ring1B_peakAnno)$geneId,                                suz12=as.data.frame(suz12_peakAnno)$geneId))  #然后提取出每个样本的在共有靶基因上的peaks,从而得到tagMatrix并绘图  cbx7_df <- as.data.frame(cbx7_peakAnno)  cbx7_peak <- GRanges(cbx7_df[cbx7_df$geneId %in% common_gene, 1:12])  ring1B_df <- as.data.frame(ring1B_peakAnno)  ring1B_peak <- GRanges(ring1B_df[ring1B_df$geneId %in% common_gene, 1:12])  suz12_df <- as.data.frame(suz12_peakAnno)  suz12_peak <- GRanges(suz12_df[suz12_df$geneId %in% common_gene, 1:12])  list_peak <- list(cbx7=cbx7_peak, ring1B=ring1B_peak, suz12=suz12_peak)  promoter <- promoters(TxDb.Mmusculus.UCSC.mm10.knownGene, upstream = 2500, downstream = 2500)  tagMatrix <- lapply(list_peak, getTagMatrix, window=promoter)  plotAvgProf(tagMatrix, xlim=c(-2500, 2499),              xlab="Genomic Region (5'->3')",               ylab = "Read Count Frequency")
      上述的代码有点冗余,但是凑合用着了;图如下所示,也跟文献的更加不像了。。。至于为啥用2499,是因为用2500报错,应该是跟我输入文件有关。。
      
      除了profile图,还可以画画热图,图就不展示了
      tagHeatmap(tagMatrix, xlim=c(-2500, 2499), color=NULL)
  • 使用deeptools可视化
      deeptools也可以一个可以对ChIP-seq结果进行可视化展示的软件,展现形式也差不多,也有profile图和热图等,下面粗略用一下   
      以cbx7样本为例,先用bamCoverage将bam文件转化为bw文件
      bamCoverage -p 4 -b cbx7.bam -o cbx7.bw
      然后用computeMatrix将peaks定位到mm10.bed文件上
      computeMatrix reference-point --referencePoint TSS -b 2500 -a 2500 -S cbx7.bw -R ~/annotation/mm10/mm10.bed --skipZeros -o cbx7_TSS.gz --outFileSortedRegions cbx7_region_genes.bed
      最后绘制全部peaks在TSS附件的热图和profile图
      plotHeatmap -m cbx7_TSS.gz -out cbx7.png  plotProfile --dpi 720 -m cbx7_TSS.gz -out cbx7_profile.pdf --plotFileFormat pdf
      如果是多个样本,可下述用computeMatrix一起读入.bw文件,再绘制图片
      computeMatrix reference-point --referencePoint TSS -b 2500 -a 2500 -S ./*.bw -R ~/annotation/mm10/mm10.bed --skipZeros -o merged_TSS.gz --outFileSortedRegions merged_region_genes.bed  plotHeatmap -m merged_TSS.gz -out merged.png  plotProfile --dpi 720 -m merged_TSS.gz -out merged_profile.pdf --plotFileFormat pdf

Summary
这次实战只是对ChIP-seq的一般流程进行初步的了解,还了解了ChIP-seq的一些基本概念。虽然从文献角度来说,这篇文献进行了些其他的研究,但是这里就不在对其重现了(一方面是作者上传数据的限制,另一方面则的能力有限~)。再者,这文献研究的是PRC1这个蛋白复合物,而ChIP-Seq似乎在组蛋白修饰和转录因子方面的研究应用更多点,可能还需找文献对这两方面进行尝试,加深下对ChIP-seq方法的理解
参考资料[size=0em]​




上一篇:科学家提倡乳腺癌单细胞诊断
下一篇:回顾基因表达芯片分析
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-3-25 10:24 , Processed in 0.037430 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.