一篇文章学会ChIP-seq分析(上)

写在前面:《一篇文章学会ChIP-seq分析(上)》《一篇文章学会ChIP-seq分析(下)》为生信菜鸟团博客相关文章合集,共九讲内容。带领你从相关文献解读、资料收集和公共数据下载开始,通过软件安装、数据比对、寻找并注释peak、寻找motif等ChIP-seq分析主要步骤入手学习,最后还会介绍相关可视化工具。

第一讲:文献选择与解读

文献;CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis

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

作者首先实验证明了用small haripin RNA来knockout CARM1 只能达到90%的敲除效果,有趣的是,对CARM1的功能影响非常小,说明只需要极少量的CARM1就可以发挥很好的作用,因此作者通过zinc finger nuclease这种基因组编辑技术设计了100%敲除CARM1的实验材料。(当然,现在有更好的基因编辑技术啦)

这样就能比较CARM1有无时各种蛋白被催化状态了,其中SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物的一个亚基 BAF155,非常明显的只有在CARM1这个基因完好无损的细胞系里面才能被正常的甲基化。作者证明了BAF155是CARM1这个基因非常好的一个底物, 而且通过巧妙的实验设计,证明了BAF155这个蛋白的第1064位氨基酸(R) 是 CARM1的作用位点。

因为早就有各种文献说明了SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物在癌症的重要作用, 所以作者也很自然想探究BAF155在癌症的功能详情,这里作者选择的是ChIP-seq技术。BAF155是作为SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物的一个组分,必然neng 直接或者间接的结合DNA咯。而ChIP-seq技术最适合来探究能直接或者间接结合DNA的蛋白的功能,所以作者构造了一种细胞系(MCF7),它的BAF155蛋白的第1064位氨基酸(R) 突变而无法被CARM1这个基因催化而甲基化,然后比较突变的细胞系和野生型细胞系的BAF155的两个ChIP-seq结果,这样就可以研究BAF155是否必须要被CARM1这个基因催化而甲基化后才能行使生物学功能。

作者用me-BAF155特异性抗体+western bloting 证明了正常的野生型MCF7细胞系里面有~74%的BAF155被甲基化。

有一个细胞系SKOV3,可以正常表达除了BAF155之外的其余14种SWI/SNF(BAF) chromatin remodeling complex 染色质重构复合物,而不管是把突变的细胞系和野生型细胞系的BAF155混在里面都可以促进染色质重构复合物的组装,所以甲基化与否并不影响这个染色质重构复合物的组装,重点应该研究的是甲基化会影响BAF155在基因组其它地方结合。

结果显示,突变的细胞系和野生型细胞系种BAF155在基因组结合位置(peaks)还是有较大的overlap的,重点是看它们的peaks在各种基因组区域(基因上下游,5,3端UTR,启动子,内含子,外显子,基因间区域,microRNA区域)分布情况的差别,还有它们距离转录起始位点的距离的分布区别,还有它们注释到的基因区别,已经基因富集到什么通路等等。

虽然作者在人的细胞系(MCF7)上面做ChIP-seq,但是在老鼠细胞系(MDA-MB-231)做了mRNA芯片数据分析,BAF155这个蛋白的第1064位氨基酸(R) 突变细胞系和野生型细胞系,用的是Affymetrix HG U133 Plus 2.0这个常用平台。

which was hybridized to Affymetrix HG U133 Plus 2.0 microarrays containing 54,675 probesets for >47,000 transcripts and variants, including 38,500 human genes.

To identify genes differentially expressed between MDA-MB-231-BAF155WT and MDA-MB-231-BAF155R1064K

表达矩阵下载地址:http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4004525/bin/NIHMS556863-supplement-03.xlsx

我简单摘抄作者ChIP-seq数据的生物信息学分析结果

用到的软件有

  • *QuEST 2.4 *(Valouev et al., 2008) was run using the recommend settings for transcription factor (TF) like binding with the following exceptions: kdebandwith=30, regionsize=600, ChIP threshold=35, enrichment fold=3, rescue fold=3.
  • *HOMER *(Heinz et al., 2010) analysis was run using the default settings for peak finding. False Discovery Rate (FDR) cut off was *0.001 (0.1%) for all peaks. * The tag density for each factor was normalized to 1x107 tags and displayed using the UCSC genome browser.
  • Motif analysis (de novo and known), was performed using the* HOMER software and Genomatix. *
  • *Peak overlaps *were processed with HOMER and Galaxy (Giardine et al., 2005).
  • *Peak comparisons *between replicates were processed with EdgeR statistical package in R

以上就是我们接下来需要学习的流程化分析步骤,下面我给一个主要流程的截图,但主要是实验是如何设计

这里有一个文章发表了关于CHIP-seq的流程的:http://biow.sb-roscoff.fr/ecolebioinfo/protected/jacques.van-helden/ThomasChollierNatProtoc2012peak-motifs.pdf

同时我还推荐大家看几篇相关文献

第二讲:资料收集

CHIP-seq的确是非常完善的NGS流程,各种资料层出不穷。

大家首先可以看下面几个完整流程的PPT来对CHIP-seq流程有个大致的印象,我对前面提到的文献数据处理的几个要点,就跟下面这个图片类似。

然后下面的各种资料,是针对CHIP-seq流程的各个环境的,还有一些是针对于表观遗传学知识

可视化工具

bioconductor系列工具和教程 :

公司教程

第三讲:公共数据下载

这一步跟自学其它高通量测序数据处理一样,就是仔细研读paper,在里面找到作者把原始测序数据放在了哪个公共数据库里面,一般是NCBI的GEO,SRA,本文也不例外,然后解析样本数,找到下载链接规律。

## step1 : download raw data> cd ~> mkdir CHIPseq_test && cd CHIPseq_test> mkdir rawData && cd rawData> ## batch download the raw data by shell script :> for ((i=593;i<601;i++)) ;do wget [ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492/SRR1042](ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492/SRR1042)$i/SRR1042$i.sra;done

很容易就下载了8个测序文件,每个样本的数据大小,测序量如下

> 621M Jun 27 14:03 SRR1042593.sra (16.9M reads)> 2.2G Jun 27 15:58 SRR1042594.sra (60.6M reads)> 541M Jun 27 16:26 SRR1042595.sra (14.6M reads)> 2.4G Jun 27 18:24 SRR1042596.sra (65.9M reads)> 814M Jun 27 18:59 SRR1042597.sra (22.2M reads)> 2.1G Jun 27 20:30 SRR1042598.sra (58.1M reads)> 883M Jun 27 21:08 SRR1042599.sra (24.0M reads)> 2.8G Jun 28 11:53 SRR1042600.sra (76.4M reads)

虽然下载的SRA格式数据也是一个很流行的标准,但它只是数据压缩的标准,几乎没有软件能直接跟SRA的格式的测序数据来进行分析,我们需要转成fastq格式,代码如下:

> ## step2 :  change sra data to fastq files.> ## cell line: MCF7 //  Illumina HiSeq 2000 //  50bp // Single ends // phred+33> ## [http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52964](http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52964)> ## [ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492](ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033492)> ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump $id;done> rm *sra

解压的详情如下,可以看到SRA格式有6~9倍的压缩了,比zip格式压缩的2~3倍高多了

##  621M --> 3.9G##  2.2G --> 14G##  541M --> 3.3G##  2.4G --> 15G

第四讲:必要软件安装及结果下载

博文的顺序有点乱,因为怕读到前面的公共测序数据下载这篇文章的朋友搞不清楚,我如何调用各种软件的,所以我这里强势插入一篇博客来描述这件事,当然也只是略过,我所有的软件理论上都是安装在我的home目录下的biosoft文件夹,所以你看到我一般安装程序都是:

cd ~/biosoftmkdir macs2 && cd macs2 ##指定的软件安装在指定文件夹里面

这只是我个人的安装习惯,因为我不是root,所以不能在linux系统下做太多事,我这里贴出我所有的软件安装代码:

## pre-step: download sratoolkit /fastx_toolkit_0.0.13/fastqc/bowtie2/bwa/MACS2/HOMER/QuEST/mm9/hg19/bedtools## http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software## http://www.ncbi.nlm.nih.gov/books/NBK158900/## Download and install sratoolkitcd ~/biosoftmkdir sratoolkit && cd sratoolkitwget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.6.3/sratoolkit.2.6.3-centos_linux64.tar.gz#### Length: 63453761 (61M) [application/x-gzip]## Saving to: "sratoolkit.2.6.3-centos_linux64.tar.gz"tar zxvf sratoolkit.2.6.3-centos_linux64.tar.gz## Download and install bedtoolscd ~/biosoftmkdir bedtools && cd bedtoolswget https://github.com/arq5x/bedtools2/releases/download/v2.25.0/bedtools-2.25.0.tar.gz## Length: 19581105 (19M) [application/octet-stream]tar -zxvf bedtools-2.25.0.tar.gzcd bedtools2make## Download and install PeakRangercd ~/biosoftmkdir PeakRanger && cd PeakRangerwget https://sourceforge.net/projects/ranger/files/PeakRanger-1.18-Linux-x86_64.zip/## Length: 1517587 (1.4M) [application/octet-stream]unzip PeakRanger-1.18-Linux-x86_64.zip~/biosoft/PeakRanger/bin/peakranger -h## Download and install bowtiecd ~/biosoftmkdir bowtie && cd bowtiewget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip/download#Length: 27073243 (26M) [application/octet-stream]#Saving to: "download" ## I made a mistake here for downloading the bowtie2mv download bowtie2-2.2.9-linux-x86_64.zipunzip bowtie2-2.2.9-linux-x86_64.zipmkdir -p ~/biosoft/bowtie/hg19_indexcd ~/biosoft/bowtie/hg19_index# download hg19 chromosome fasta fileswget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz# unzip and concatenate chromosome and contig fasta filestar zvfx chromFa.tar.gzcat *.fa > hg19.farm chr*.fa## ~/biosoft/bowtie/bowtie2-2.2.9/bowtie2-build ~/biosoft/bowtie/hg19_index/hg19.fa ~/biosoft/bowtie/hg19_index/hg19## Download and install BWAcd ~/biosoftmkdir bwa && cd bwahttp://sourceforge.net/projects/bio-bwa/files/tar xvfj bwa-0.7.12.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 filescd bwa-0.7.12makeexport PATH=$PATH:/path/to/bwa-0.7.12 # Add bwa to your PATH by editing ~/.bashrc file (or .bash_profile or .profile file)# /path/to/ is an placeholder. Replace with real path to BWA on your machinesource ~/.bashrc# bwa index [-a bwtsw|is] index_prefix reference.fastabwa index -p hg19bwaidx -a bwtsw ~/biosoft/bowtie/hg19_index/hg19.fa# -p index name (change this to whatever you want)# -a index algorithm (bwtsw for long genomes and is for short genomes)## Download and install macs2## // https://pypi.python.org/pypi/MACS2/cd ~/biosoftmkdir macs2 && cd macs2wget ~~~~~~~~~~~~~~~~~~~~~~MACS2-2.1.1.20160309.tar.gztar zxvf MACS2-2.1.1.20160309.tar.gzcd MACS2-2.1.1.20160309python setup.py install --user#################### The log for installing MACS2:Creating ~/.local/lib/python2.7/site-packages/site.pyProcessing MACS2-2.1.1.20160309-py2.7-linux-x86_64.eggCopying MACS2-2.1.1.20160309-py2.7-linux-x86_64.egg to ~/.local/lib/python2.7/site-packagesAdding MACS2 2.1.1.20160309 to easy-install.pth fileInstalling macs2 script to ~/.local/binFinished processing dependencies for MACS2==2.1.1.20160309############################################################~/.local/bin/macs2 --helpExample for regular peak calling:macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01Example for broad peak calling:macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1## Download and install homer (Hypergeometric Optimization of Motif EnRichment)## // http://homer.salk.edu/homer/## // http://blog.qiubio.com:8080/archives/3024## pre-install: Ghostscript,seqlogo,blatcd ~/biosoftmkdir homer && cd homerwget http://homer.salk.edu/homer/configureHomer.plperl configureHomer.pl -installperl configureHomer.pl -install hg19

一般来说,对我这样水平的人来说,软件安装就跟家常便饭一样,没有什么问题了,但如果你是初学者呢,肯定没那么轻松,所以请加强学习,我无法在这里讲解太具体的知识了。

所有软件安装完毕后就可以下载文章对这些ChIP-seq的处理结果了,这个很重要,检验我们是否重复了人家的数据分析过程。

## step3 : download the results from paper## http://www.bio-info-trainee.com/1571.htmlmkdir paper_results && cd paper_resultswget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE52nnn/GSE52964/suppl/GSE52964_RAW.tartar xvf GSE52964_RAW.tarls *gz |xargs gunzip## step4 : run FastQC to check the sequencing quality.##这里可以看到我们下载的原始数据已经被作者处理好了,去了接头,去了低质量序列ls *.fastq | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;done## Sequence length 51## %GC 39## Adapter Content passedThe quality of the reads is pretty good, we don't need to do any filter or trimmkdir QC_resultsmv *zip *html QC_results/

编辑校对:思考问题的熊

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-05-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏CDA数据分析师

Python告诉你:单词软件火了,但真的有那么多人在背单词吗?

0x00 前言 你想知道背单词软件有大概多少人注册第一天都没有背完嘛? 你想知道背单词软件这么火,这么多人在使用,真的有多少人真的在背诵嘛? 别急,Python...

1817
来自专栏PPV课数据科学社区

【完整案例】如何用R实现空间数据可视化

image.png 流行病学的数据讲究“三间分布”,即人群分布、时间分布和空间分布。其中的“空间分布”最好是在地图上展示,才比较清楚。R软件集统计分析与高级...

4787
来自专栏FreeBuf

利用威胁情报数据平台拓展APT攻击线索一例

当我们说起APT攻击线索的发现,似乎是一个挺神秘的事,安全厂商往往说得云山雾罩,如果现在你问如何知道某件事情的时侯,得到的回答往往是:”嗯,我们用了机器学习”,...

3076
来自专栏人工智能头条

无需一行代码就能搞定机器学习的开源神器

1142
来自专栏AI科技大本营的专栏

无需一行代码就能搞定机器学习的开源神器

作者 | Shantanu Kumar 责编 | 魏伟 对于机器学习和数据科学的初学者来说,最大的挑战之一是需要同时学习太多知识,特别是如果你不知道如何编码。你...

4487
来自专栏大数据挖掘DT机器学习

R语言绘制中国地图,并展示流行病学数据

本文作者:姜晓东,博士毕业于上海交通大学,目前任教于湖南师范大学医学院,专业神经毒理学。 流行病学的数据讲究“三间分布”,即人群分布、时间分布和空间分布。其中...

6586
来自专栏Data Analysis & Viz

手把手教你完成一个数据科学小项目(8):Emoji提取与分布图谱

请先阅读“中国年轻人正带领国家走向危机”,这锅背是不背? 一文,以对“手把手教你完成一个数据科学小项目”系列有个全局性的了解。

962
来自专栏安富莱嵌入式技术分享

【安富莱专题教程第1期】基于STM32的硬件RGB888接口实现emWin的快速刷新方案,32位色或24

说明: 1. 首先感谢ST终于推出了ARGB格式的emWin库,可谓千呼万唤始出来,使用STM32的硬件RGB888接口刷新图片慢的问题终于得到解决。 2. 这...

701
来自专栏生信技能树

【直播】我的基因组52:X和Y染色体的同源区域探索

很久以前,我其实就遇到过通过NGS测序数据来判定性别的难题(搜索我博客即可查看详情),本次探究自己的基因组得到的统计结果与常识不符,所以我可以肯定是我们的常识太...

28811
来自专栏生信技能树

【直播】我的基因组58:用R包SNPRelate来对我的基因型跟hapmap计划数据比较

hapmap计划的人群分布结果和千人基因组计划的分布结果来分析是一样的!【直播】我的基因组55:简单的PCA分析千人基因组的人群分布 这两个计划里面收集的样本的...

4116

扫码关注云+社区