我的第一次ChIP-seq实践

1. 软件安装

整个过程基本是从零开始,也就是说服务器没有安装任何所需软件。因为我平时会用到Python,所以第一步安装的是Anaconda,版本是Anaconda3-4.4.0-Linux-x86_64.sh. Ananconda是一个用于科学计算的Python发行版,能够方便解决多版本Python并存(后面会看到)、切换以及各种第三方包安装问题(最大的好处)。

Anaconda和整个ChIP-seq分析没关系,提到它是因为安装Anaconda后可以用BIOCONDA,能够方便安装管理生物信息软件,无需自己解决软件之间依赖关系。

安装本次实践需要的软件,包括

fastqc(0.11.5), sra-tools(2.8.1), bowtie2(2.3.2), samtools(1.5), MACS2(2.1.1.20160309), deeptools(2.5.1)

    conda config --add channels defaults    conda config --add channels conda-forge    conda config --add channels bioconda    conda install fastqc    conda install sra-tools    conda install bowtie2    conda install samtools    conda install deeptools

MACS2比较特殊,因为它需要Python2.7的环境才能安装,因为我习惯使用Python3,所以安装的Anaconda版本是Python3.6的。为了安装MACS2,首先创建一个虚拟环境,在python2.7的虚拟环境中使用。

    conda create -n py27 python=2    source activate py27    pip install MACS2

用pip没用conda,是因为conda install MACS2需要装好多其他包,我就按照MACS2网站提示pip安装,发现也能用。

2. 数据下载

处理过程中需要下载必要的数据,参考生物信息学常见的数据下载(http://www.biotrainee.com/thread-857-1-1.html)

mm10基因组

    mkdir -p  ~/reference/genome/mm10    cd ~/reference/genome/mm10    wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz    tar zvxf 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/genome/mm10/mm10.fa  ~/reference/index/bowtie/mm10 1>mm10.bowtie_index.log 2>&1 &

mm10 refseq注释文件,来自http://genome.ucsc.edu/cgi-bin/hgTables。其中track选择RefSeq Genes,output format选择BED,其他参数默认。 头几行得到结果如下。

3. 处理数据

处理的步骤和Jimmy一模一样。 首先需要下载数据,在Jimmy原文中写的清楚,这里也就不赘述了。但是原文截取的原始数据的图好像不太对。我得到的数据如图:

PS:原文截图的确出错了,原文用的是RNA-seq实战的图片,编辑错误。

然后用fastqc看数据质量

    ls *fastq|xargs fastqc -t 10

用Bowtie2将fastqc结果比对到mm10基因组上去

    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620205.fastq| samtools sort -O bam -o cbx7.bam    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620206.fastq| samtools sort -O bam -o suz12.bam    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620207.fastq| samtools sort -O bam -o RYBP.bam    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620208.fastq| samtools sort -O bam -o IgGold.bam    bowtie2 -p 6 -3 5 --local -x ~/reference/index/bowtie/mm10 -U ../data/SRR620209.fastq| samtools sort -O bam -o IgG.bam

得到比对结果如图。大小和原文已经不一样,可能是index,bowtie2版本不同有关。

然后call peaks,首先要进入py27的虚拟环境中

    source activate py27    macs2 callpeak -c IgGold.bam -t suz12.bam -m 10 30 -p 1e-5 -f BAM -g mm -n suz12 2>suz12.masc2.log    macs2 callpeak -c IgGold.bam -t ring1B.bam -m 10 30 -p 1e-5 -f BAM -g mm -n ring1B 2>ring1B.masc2.log    macs2 callpeak -c IgG.bam -t cbx7.bam -m 10 30 -p 1e-5 -f BAM -g mm -n cbx7 2>cbx7.masc2.log    macs2 callpeak -c IgG.bam -t RYBP.bam -m 10 30 -p 1e-5 -f BAM -g mm -n RYBP 2>RYBP.masc2.log

结果如下。和原文结果很接近。但是原文RYBP2不知道是怎么得到的。

PS:我的RYBP2得到的方式,其实原文有说!

然后是bam转换成bw文件,并画图

    ls *.bam |while read id    do    samtools index $id $id.bai    done    ls *bam |while read id    do    file=$(basename $id )    sample=${file%%.*}    echo $sample    bamCoverage -b $id -o $sample.bw -p 10 --normalizeUsingRPKM    computeMatrix reference-point --referencePoint TSS -b 10000 -a 10000 -R ~/annotation/CHIPseq/mm10/ucsc.refseq.bed -S $sample.bw --skipZeros -o matrix1_${sample}_TSS.gz --outFileSortedRegions regions1_${sample}_genes.bed    plotHeatmap -m matrix1_${sample}_TSS.gz -out ${sample}.png    done

结果如下。其中跑bamCoverage会报错ring1B.bam does not appear to have an index. You MUST index the file first!

先要跑samtools index,生成.bam.bai文件。这个不知道为什么跑原代码没有生成这个文件。

最后画出profile和heatmap图。

    computeMatrix reference-point -p 10 --referencePoint TSS -b 2000 -a 2000 -S *bw -R ~/annotation/CHIPseq/mm10/ucsc.refseq.bed --skipZeros -o tmp4.mat.gz    plotHeatmap -m tmp4.mat.gz -out tmp4.merge.png    plotProfile --dpi 720 -m tmp4.mat.gz -out tmp4.profile.pdf --plotFileFormat pdf --perGroup    plotHeatmap --dpi 720 -m tmp4.mat.gz -out tmp4.merge.pdf --plotFileFormat pdf

得到结果如下

和原文差别好多。后来发现是bamCoverage不需要 --normalizeUsingRPKM这个参数,调整后结果是

现在和原文结果基本一致了。

4.总结

本次实践过程,只是单纯把软件安装,数据下载,代码跑了一遍,其中参数的设置及含义也没有仔细查看,这个要在之后的练习中不断学习。 多谢Jimmy分享的教程,能够让读者快速进入实战分析!

编辑:思考问题的熊&&jimmy

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

原文发表时间:2017-07-14

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏FreeBuf

Android 5.x漏洞:黑客可以绕过屏幕密码进入系统

很多Android用户会选择使用锁屏密码保护设备,但最新爆出的漏洞却令人震惊:任何人无需复杂的操作即可绕过锁屏直接进入你的系统! 攻击者可以通过漏洞导获取上锁...

19710
来自专栏FreeBuf

打开文件夹就运行?COM劫持利用新姿势

*本文原创作者:菠菜,本文属FreeBuf原创奖励计划,未经许可禁止转载 打开文件夹就能运行指定的程序?这不是天方夜谭,而是在现实世界中确实存在的。利用本文探讨...

2068
来自专栏Samego开发资源

让子弹飞~利用OPcache扩展提升PHP7性能 | laravel篇

What is OPcache OPcache是PHP中的Zend扩展,可以大大提升PHP的性能。 OPcache 通过将 PHP 脚本预编译的字节码存储到...

2082
来自专栏杨建荣的学习笔记

Oracle 12C打补丁的简单尝试(r10笔记第55天)

最近在服务器盘点的时候,发现测试环境还是值得整合一下,因为服务器资源老旧,整体配置不高,服务器资源使用率不高,业务要求不高,多个实例分散在多台服务器上,要考虑灾...

3478
来自专栏杨建荣的学习笔记

gc服务器慢的原因分析 (r6笔记第14天)

在工作环境中有一台gc的服务器,已经好几年没有动过了,上面安装着gc的服务和数据库,也就说gc里面的HttpServer,数据库,webcache都在这台服务器...

2743
来自专栏Danny的专栏

重装Win7时提示“缺少所需的CD/DVD驱动器设备驱动程序”

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/huyuyang6688/article/...

5.3K2
来自专栏FreeBuf

甲方安全中心建设:代码审计系统

纵观甲方的安全体系建设,最开始和最重要的那一部分就是代码安全。甲方公司内部有很多项目,每个项目都由不同的开发人员进行开发,所以项目开发水平也是参差不齐,也就是说...

1062
来自专栏北京马哥教育

保护你的Linux系统的九个老生常谈

在现在这个世道中,保障基于Linux的系统的安全是十分重要的。但是,你得知道怎么干。一个简单反恶意程序软件是远远不够的,你需要采取其它措施来协同工作。那么试试...

2786
来自专栏程序你好

微服务和传统中间件平台

微服务与部署在中间件平台(esb、应用服务器)上的传统服务有何不同?什么是微服务体系结构模式,它解决了什么问题?本文将讨论所有这些重要的主题,并描述如何管理、管...

962
来自专栏我的安全视界观

[一起玩蛇】Python代码审计中的器II

2357

扫码关注云+社区