07

用网页版工具GREAT来对CHIP-seq的peaks进行下游功能分析

一般做完一个CHIP-seq测序,如果实验设计没有问题,测序质量也OK的话,很容易了根据序列call到符合要求的peaks,或者可以去很多文章或者roadmap里面下载到非常多有意义的peaks文件, 一般是BED格式文件,这是就需要对这些peaks进行各种各样的注释以及可视化了,还有根据peaks相关的基因可以做各种各样的下游分析,包括各种pathway数据库的富集,MsigDB数据库注释,gene ontology的注释等等,此时不得不强烈推荐一款网页版工具,是斯坦福大学的学者开发的GREAT。
此工具的出现主要是为了解决基因组上面的非编码区域注释缺乏的问题,而我们CHIP-seq实验得到的peaks结果通常就是在非编码区域
该工具每次只能上传一个文件,就是我们call出来的peaks记录文件,支持bed格式的:

Continue reading

07

用网页版工具ChIPseek来可视化CHIP-seq的peaks结果

一般做完一个CHIP-seq测序,如果实验设计没有问题,测序质量也OK的话,很容易了根据序列call到符合要求的peaks,或者可以去很多文章或者roadmap里面下载到非常多有意义的peaks文件, 一般是BED格式文件,这是就需要对这些peaks进行各种各样的注释以及可视化了,此时不得不强烈推荐一款网页版工具,是台湾学者开发的ChIPseek:

Continue reading

06

自学CHIP-seq分析第七讲~peaks注释

经过前面的CHIP-seq测序数据处理的常规分析,我们已经成功的把测序仪下机数据变成了BED格式的peaks记录文件,我选取的这篇文章里面做了4次CHIP-seq实验,分别是两个重复的野生型MCF7细胞系的 BAF155 immunoprecipitates和两个重复的突变型MCF7细胞系的 BAF155 immunoprecipitates,这样通过比较野生型和突变型MCF7细胞系的 BAF155 immunoprecipitates的结果的不同就知道该细胞系的BAF155 突变,对它在全基因组的结合功能的影响啦。

#我这里直接从GEO里面下载了peaks结果,它们详情如下:wc -l *bed
6768 GSM1278641_Xu_MUT_rep1_BAF155_MUT.peaks.bed
3660 GSM1278643_Xu_MUT_rep2_BAF155_MUT.peaks.bed
11022 GSM1278645_Xu_WT_rep1_BAF155.peaks.bed
5260 GSM1278647_Xu_WT_rep2_BAF155.peaks.bed
49458 GSM601398_Ini1HeLa-peaks.bed
24477 GSM601398_Ini1HeLa-peaks-stringent.bed
12725 GSM601399_Brg1HeLa-peaks.bed
12316 GSM601399_Brg1HeLa-peaks-stringent.bed
46412 GSM601400_BAF155HeLa-peaks.bed
37920 GSM601400_BAF155HeLa-peaks-stringent.bed
30136 GSM601401_BAF170HeLa-peaks.bed
25432 GSM601401_BAF170HeLa-peaks-stringent.bed

每个BED的peaks记录,本质是就3列是需要我们注意的,就是染色体,以及在该染色体上面的起始和终止坐标,如下:

#PeakID chr start end strand Normalized Tag Count region size findPeaks Score Clonal Fold Change
chr20 52221388 52856380 chr20-8088 41141 +
chr20 45796362 46384917 chr20-5152 31612 +
chr17 59287502 59741943 chr17-2332 29994 +
chr17 59755459 59989069 chr17-667 19943 +
chr20 52993293 53369574 chr20-7059 12642 +
chr1 121482722 121485861 chr1-995 9070 +
chr20 55675229 55855175 chr20-6524 7592 +
chr3 64531319 64762040 chr3-4022 7213 +
chr20 49286444 49384563 chr20-4482 6165 +

我们所谓的peaks注释,就是想看看该peaks在基因组的哪一个区段,看看它们在各种基因组区域(基因上下游,5,3端UTR,启动子,内含子,外显子,基因间区域,microRNA区域)分布情况,但是一般的peaks都有近万个,所以需要批量注释,如果脚本学的好,自己下载参考基因组的GFF注释文件,完全可以自己写一个,我这里会介绍一个R的bioconductor包ChIPpeakAnno来做CHIP-seq的peaks注释,下面的包自带的示例:

library(ChIPpeakAnno)
bed <- system.file("extdata", "MACS_output.bed", package="ChIPpeakAnno")
gr1 <- toGRanges(bed, format="BED", header=FALSE)
## one can also try import from rtracklayer
library(rtracklayer)
gr1.import <- import(bed, format="BED")
identical(start(gr1), start(gr1.import))
gr1[1:2]
gr1.import[1:2] #note the name slot is different from gr1
gff <- system.file("extdata", "GFF_peaks.gff", package="ChIPpeakAnno")
gr2 <- toGRanges(gff, format="GFF", header=FALSE, skip=3)
ol <- findOverlapsOfPeaks(gr1, gr2)
makeVennDiagram(ol)

##还可以用binOverFeature来根据特定的GRanges对象(通常是TSS)来画分布图
## Distribution of aggregated peak scores or peak numbers around transcript start sites.

可以看到这个包使用起来非常简单,只需要把我们做好的peaks文件(GSM1278641_Xu_MUT_rep1_BAF155_MUT.peaks.bed等等)用toGRanges或者import读进去,成一个GRanges对象即可,上面的代码是比较两个peaks文件的overlap。然后还可以根据R很多包都自带的数据来注释基因组特征:

data(TSS.human.GRCh37) ## 主要是借助于这个GRanges对象来做注释,也可以用getAnnotation来获取其它GRanges对象来做注释
## featureType : TSS, miRNA, Exon, 5'UTR, 3'UTR, transcript or Exon plus UTR
peaks=MUT_rep1_peaks
macs.anno <- annotatePeakInBatch(peaks, AnnotationData=TSS.human.GRCh37,
output="overlapping", maxgap=5000L)

## 得到的macs.anno对象就是已经注释好了的,每个peaks是否在基因上,或者距离基因多远,都是写的清清楚楚
if(require(TxDb.Hsapiens.UCSC.hg19.knownGene)){
aCR<-assignChromosomeRegion(peaks, nucleotideLevel=FALSE,
precedence=c("Promoters", "immediateDownstream",
"fiveUTRs", "threeUTRs",
"Exons", "Introns"),
TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)
barplot(aCR$percentage)
}

得到的条形图如下,虽然很丑,但这就是peaks注释的精髓,搞清楚每个peaks在基因组的位置特征:

ChIPpeakAnno-genomic-region-distribution

 

同理,对每个peaks文件,都可以做类似的分析!

但是对多个peaks文件,比如本文中的,想比较野生型和突变型MCF7细胞系的 BAF155 immunoprecipitates的结果的不同,就需要做peaks之间的差异分析,已经后续的差异基因注释啦

当然,值得一提的是peaks注释我更喜欢网页版工具,反正peaks文件非常小,直接上传到别人做好的web tools,就可立即出一大堆可视化图表分析结果啦,大家可以去试试看:

虽然我花费了大部分篇幅来描述ChIPpeakAnno这个包的用法,但是真正的重点是你得明白peaks记录了什么,要注释什么,已经把这3个网页工具的可视化图表分析结果全部看懂,这网页版工具才是重点!!!
05

用R包BayesPeak来对CHIP-seq数据call peaks

BayesPeak也是peaks caller家族一员,用的人也不少,我这次也试了一下,因为是R的bioconductor系列包,所以直接在R里面安装就好,但是有几个点需要注意,我比对的基因组不只是Chr1~22,X,Y,M,还有一些contig和scaffold,需要在bam文件里面去除的,而且BayesPeak比较支持读取BED文件,可以直接转为GRanges对象,虽然它号称可以使用多核,但是计算速度还是非常慢。 Continue reading

05

用PeakRanger软件来对CHIP-seq数据call peaks

此文专门讲这个软件如何用,但是跟我以前写的软件说明书又不大一样,主要是因为我用MACS2这个软件call peaks并没有达到预期的结果,所以就多使用了几个软件,其中PeakRanger尤其值得一提,安装特别简单,而且处理数据的速度特别快,结果也非常容易理解,更重要的是它给出一个网页版的报告,里面有所有找到的符合要求的peaks的可视化图片!!!! Continue reading

05

自学CHIP-seq分析第二讲~学习资料的搜集

我只能说,CHIP-seq的确是非常完善的NGS流程了,各种资料层出不穷,大家首先可以看下面几个完整流程的PPT来对CHIP-seq流程有个大致的印象,我对前面提到的文献数据处理的几个要点,就跟下面这个图片类似:
chip-seq-workflow-all-5-steps
QuEST is a statistical software for analysis of ChIP-Seq data with data and analysis results visualization through UCSC Genome Browser.  http://www-hsc.usc.edu/~valouev/QuEST/QuEST.html
MeDIP-seq and histone modification ChIP-seq analysis  http://crazyhottommy.blogspot.com/2014/01/medip-seq-and-histone-modification-chip.html
一个实际的CHIP-seq数据分析例子: http://www.biologie.ens.fr/~mthomas/other/chip-seq-training/

Continue reading

05

自学CHIP-seq分析第一讲~文献选择与解读

文章:CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis 

我很早以前想自学CHIP-seq的时候就关注过这篇文章,那时候懂得还不多,甚至都没有仔细看这篇文章就随便下载了数据进行分析,也只是跑一些软件而已,这次仔细阅读这篇文章才发现里面的门道很多,尤其是CHIP-seq的实验基础,以及表观遗传学的生物学基础知识,我有时间一定要把这篇文章翻译一下。学习这篇文章前一定要温习一些生物学知识,见我上一篇博客

Continue reading