一个MeDIP-seq实战-超级简单-2小时搞定!

Featured

请不要直接拷贝我的代码,需要自己理解,然后打出来,思考我为什么这样写代码。
软件请用最新版,尤其是samtools等被我存储在系统环境变量的,考虑到读者众多,一般的软件我都会自带版本信息的!
我用两个小时,不代表你是两个小时就学会,有些朋友反映学了两个星期才 学会,这很正常,没毛病,不要异想天开两个小时就达到我的水平。

MeDIP-seq 跟ChIP-seq的分析手段是一模一样的,同理hMeDIP-seq,caMeDIP-seq等等,都没有本质上的区别,只是用的抗体不一样而已,请自行搜索基础知识,我只讲数据分析。

一个ChIP-seq实战-超级简单-2小时搞定!

一个RNA-seq实战-超级简单-2小时搞定!

Continue reading

10

一个ChIP-seq实战-超级简单-2小时搞定!

请不要直接拷贝我的代码,需要自己理解,然后打出来,思考我为什么这样写代码。
软件请用最新版,尤其是samtools等被我存储在系统环境变量的,考虑到读者众多,一般的软件我都会自带版本信息的!
我用两个小时,不代表你是两个小时就学会,有些朋友反映学了两个星期才 学会,这很正常,没毛病,不要异想天开两个小时就达到我的水平。
本次讲解选取的文章是为了探索PRC1,PCR2这样的蛋白复合物,不是转录因子或者组蛋白的CHIP-seq,请注意区别!
这是一个系列帖子,你可以先看:
文章是:RYBP and Cbx7 define specific biological functions of polycomb complexes in mouse embryonic stem cells
RYBP and Cbx7都是Polycomb repressive complex 1 (PRC1)的组分:
所以用脚本在ftp里面批量下载即可:
1

Continue reading

01

ngsplot辅助CHIP-seq数据分析-可视化

最近在忙一些chip-seq的数据分析项目,它的可视化展现比较复杂一点,自己写程序将会耗费挺长时间的,就想着利用现成的工具,前面试用了deeptools,挺好 的,但是有点慢,是python程序,如下:
现在换一个R程序,这个非常快速,而且绘图个人觉得稍微美观一点,大家也可以都试试看。
首先软件的github里面有源代码,然后作者还四处宣讲这个包的神奇之处,下面的ppt非常言简意赅的描述了它的功能和强大之处。

Continue reading

十二 30

一个RNA-seq实战-超级简单-2小时搞定!

请不要直接拷贝我的代码,需要自己理解,然后打出来,思考我为什么这样写代码。
软件请用最新版,尤其是samtools等被我存储在系统环境变量的,考虑到读者众多,一般的软件我都会自带版本信息的!
我用两个小时,不代表你是两个小时就学会,有些朋友反映学了两个星期才 学会,这很正常,没毛病,不要异想天开两个小时就达到我的水平。
转录组如果只看表达量真的是超级简单,真是超级简单,而且人家作者本来就测是SE50,这种破数据,也就是看表达量用的!
首先作者分析结果是:
1

Continue reading

十二 28

生物信息数据分析文章就是看图写作文

首先是从测试原始数据里面得到汇总数据
然后把各种统计汇总数据可视化成图表
最后根据图表来写作文即可。
来源:Genome-wide Mapping of HATs and HDACs Reveals Distinct Functions in Active and Inactive Genes

Continue reading

十一 25

hisat2+stringtie+ballgown

早在去年九月,我就写个博文说 RNA-seq流程需要进化啦! http://www.bio-info-trainee.com/1022.html  ,主要就是进化成hisat2+stringtie+ballgown的流程,但是我一直没有系统性的讲这个流程,因为我觉真心木有用。我只用了里面的hisat来做比对而已!但是群里的小伙伴问得特别多,我还是勉为其难的写一个教程吧,你们之间拷贝我的代码就可以安装这些软件的!然后自己找一个测试数据,我的脚本很容易用的! Continue reading

十一 15

htseq-counts跟bedtools的区别

我以前写过bedtools和htseq-counts的教程,它们都可以用来对比对好的bam文件进行计数,正好群里有小伙伴问我它们的区别,我就简单做了一个比较,大家可以先看看我以前写的软件教程。写的有的挫:

使用Bedtools对RNA-seq进行基因计数 ,

转录组HTseq对基因表达量进行计数

言归正传,我这里没精力去探究它们的具体原理,只是看看它们数一个read是否属于某个基因的时候,区别在哪里,大家看下图: Continue reading

15

用lumi包来处理illumina的bead系列表达芯片

表达芯片大家最熟悉的当然是affymetrix系列芯片啦,而且分析套路很简单,直接用R的affy包,就可以把cel文件经过RMA或者MAS5方法得到表达矩阵。illumina出厂的芯片略微有点不一样,它的原始数据有3个层级,一般拿到的是Processed data (示例), 当仍然需要一系列的统计学方法才能提取到表达矩阵。我比较喜欢用bioconductor,所以下面讲一讲如何用lumi包来处理这个芯片数据!

这个lumi包的使用代码和说明书都有,按部就班的学一遍就好了。
如果仅仅是分析数据,那么并不难,但是每个分析步骤后面都隐含着一系列的统计学方法,想彻底搞清楚他它们, 就很难了。

Continue reading

15

illumina的bead 系列表达芯片扫盲

表达芯片大家最熟悉的当然是affymetrix系列芯片啦,而且分析套路很简单,直接用R的affy包,就可以把cel文件经过RMA或者MAS5方法得到表达矩阵。illumina出厂的芯片略微有点不一样,它的原始数据有3个层级,一般拿到的是Processed data (示例), 当仍然需要一系列的统计学方法才能提取到表达矩阵。接下来我们首先讲一讲illumina的bead 系列表达芯片基础知识吧: Continue reading

09

RNAseq数据完整生物信息分析流程第一讲之文献数据下载

我这里拿的是bioconductor里面最常用的airway数据,因为差异表达分析在bioconductor里面是重点,它们这些包在介绍自己的算法以及做示范的时候都用的这个数据。可以在GEO数据库里面看到信息描述:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778  可以看到是Illumina HiSeq 2000 (Homo sapiens) ,75bp paired-end 这个信息很重要,决定了下载sra数据之后如何解压以及如何比对。也可以看到作者把所有的测序原始数据都上传到了SRA中心:http://www.ncbi.nlm.nih.gov/sra?term=SRP033351  ,这里可以在linux服务器上面写一个简单的脚本批量下载所有的测序数据,然后根据GEO里面描述的metadata把原始数据改名。

for ((i=508;i<=523;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033351/SRR1039$i/SRR1039$i.sra;done
ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 $id;done

需要自己看SRA里面的数据记录,上面的脚本不难写出,然后因为是Illumina的双端测序,所以我们用fastq-dump --split-3命令来把sra格式数据转换为fastq,但是因为这里有16个测序数据,所以最好是同步改名,我这里用脚本批量生成改名脚本如下:

为了节省空间,我用了--gzip压缩,该文件名,用-A参数。

nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_untreated SRR1039508.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Dex SRR1039509.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb SRR1039510.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb_Dex SRR1039511.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_untreated SRR1039512.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Dex SRR1039513.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb SRR1039514.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb_Dex SRR1039515.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_untreated SRR1039516.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Dex SRR1039517.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb SRR1039518.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb_Dex SRR1039519.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_untreated SRR1039520.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Dex SRR1039521.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb SRR1039522.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb_Dex SRR1039523.sra &

可以看到这里的16个样本来源于同样的4个人,是HASM细胞系,处理详情如下:

测序基础:
HASM细胞系-human airway smooth muscle,
The Illumina TruSeq assay was used to prepare 75bp paired-end libraries for HASM cells from four white male donors under four treatment conditions:
1) no treatment;
2) treatment with a β2-agonist (i.e. Albuterol, 1μM for 18h);
3) treatment with a glucocorticosteroid (i.e. Dexamethasone (Dex), 1μM for 18h);
4) simultaneous treatment with a β2-agonist and glucocorticoid
and the libraries were sequenced with an Illumina Hi-Seq 2000 instrument.
我们这里只是先根据fastq数据比对到参考基因组,然后计算每个样本的表达量即可,后续的分组计算差异表达,就需要个性化了。

下载的sra大小如下:

-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 04:21 SRR1039508.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 05:20 SRR1039509.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 06:14 SRR1039510.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 07:05 SRR1039511.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 08:07 SRR1039512.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 09:17 SRR1039513.sra
-rw-rw-r-- 1 jmzeng jmzeng 3.1G Aug 9 10:56 SRR1039514.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 11:56 SRR1039515.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 13:02 SRR1039516.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.6G Aug 9 14:16 SRR1039517.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 15:17 SRR1039518.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 16:05 SRR1039519.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 16:56 SRR1039520.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.4G Aug 9 17:57 SRR1039521.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 18:46 SRR1039522.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 19:28 SRR1039523.sra

解压后成双端测序的fastq数据如下:

 -rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 289M Aug 9 20:44 N052611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 951M Aug 9 20:59 N052611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 954M Aug 9 20:59 N052611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 16M Aug 9 20:45 N061011_Alb.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 759M Aug 9 20:00 N061011_untreated.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 535M Aug 9 19:59 N080611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_2.fastq.gz

接下来所有的分析就基于此数据啦

 

 

07

对CHIP-seq数据call peaks应该选取unique比对的reads吗?

对于CHIP-seq数据处理完全是自学的,所以有很多细节得慢慢学习回来,这次记录的就是当我们把测序仪的fastq数据比对到参考基因组之后,应该对比对的结果文件做什么样的处理,然后去给peaks caller软件拿来call peaks呢?我看过博客 提到只保留比对质量值大于30的,也看过博客提到只保留unique比对的reads,我这里拿一篇公共数据测试了一下它们的区别!数据描述如下: Continue reading

02

根据比对的bam文件来对peaks区域可视化

之前分析了好几个公共项目,拿到的peaks都很诡异,搞得我一直怀疑是不是自己分析错了。终于,功夫不负有心人,我分析了一个数据,它的peaks非常完美!!!可以证明,我的分析流程以及peaks绘图代码并没有错!数据来自于http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74311,是关于H3K27ac_ChIP-Seq_LOUCY,组蛋白修饰的CHIP-seq数据,很容易就下载了作者上传的测序数据,然后跑了我的流程!https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq Continue reading

28

4种方式下载roadmap计划的所有数据

精选的129个细胞系,细胞系的介绍如下:http://www.broadinstitute.org/~anshul/projects/roadmap/metadata/EID_metadata.tab
对每个细胞系,都至少处理了5个核心组蛋白修饰数据,还有其它若干转录因子数据。
官网介绍的很详细,我就不翻译了:

Continue reading

28

6种方式下载ENCODE计划的所有数据

DNA元件百科全书(Encyclopedia of DNA Elements, ENCODE)ENCODE计划的重要性我就不多说了,如果大家还不是很了解,可以直接跳到本文末尾去下载一下ENCODE教程,好好学习。该计划采用以下几种高通量测序技术来刻画了超过100种不同的细胞系或者组织内的全基因组范围内的基因调控元件信息。本来只是针对人类的,后来对mouse以及fly等模式生物也开始测这些数据并进行分析了, 叫做 modENCODE
  • chromatin structure (5C)
  • open chromatin (DNase-seq and FAIRE-seq)
  • histone modifications and DNA-binding of over 100 transcription factors (ChIP-seq)
  • RNA transcription (RNAseq and CAGE)

Continue reading

26

用UCSC提供的Genome Browser工具来可视化customTrack

customTrack,我这里翻译为自定义的测序片段示踪文件,可以追踪我们的reads到底比对到了参加基因组的什么区域,或者追踪参考基因组的各个区域的覆盖度,测序深度!翻译自:http://genome.ucsc.edu/goldenPath/help/customTrack.html  这个非常有用!!!
UCSC提供的Genome Browser工具非常好用,可以很方便的浏览我们的测序数据在参考基因组的比对情况,由于定义好了一系列track的文件格式,用户可以非常方便的上传自己的track文件,但是如果用户超过48小时没有浏览自己的数据,UCSC会默认删除掉这些数据,除非用户已经保存在session里面。或者用户可以分享这些自定义的reads示踪文件customTrack。

Continue reading

26

wig、bigWig和bedgraph文件详解

我们一般会熟悉sam/bam格式文件,就是把测序reads比对到参考基因组后的文件!bam或者bed格式的文件主要是为了追踪我们的reads到底比对到了参加基因组的什么区域,而UCSC规定的这几个文件格式(wig、bigWig和bedgraph)用处不一样,仅仅是为了追踪参考基因组的各个区域的覆盖度,测序深度!而且这些定义好的文件,可以无缝连接到UCSC的Genome Browser工具里面进行可视化!
这个网站提供了这几种数据格式的构造及转换脚本:http://barcwiki.wi.mit.edu/wiki/SOPs/coordinates
对SE数据,可以用macs2 pileup --extsize 200 -i $sample.bam -o $sample.bdg 把bam文件转换为bedgraph文件,不需要call peaks这一步骤。
而UCSC的ftp里面可以下载bedGraphToBigWig $sample.bdg ~/reference/genome/mm10/mm10.chrom.sizes $sample.bw 把bedgraph文件转换为bw文件,其余的转换工具都可以下载。

Continue reading

14

ChIP-Seq文献数据重新分析解读第二例

这篇文章是朋友推荐的, 我觉得作为CHIP-seq学习材料再好不过了,所以推荐给大家。是全基因组范围的BRCA1和PALB2的转录共激活机制的探究。请务必先看我的CHIP-seq自学系列教程,跟着好好学习!数据如下:
GSM997540    BRCA1    SRR553473.sra    Read 18878514 spots
GSM997541    PALB2    SRR553474.sra    Read 17615498 spots
GSM997542    P_Ser2    SRR553475.sra    Read 35396009 spots

Continue reading

13

ChIP-Seq文献数据重新分析解读第一例

文章是:Genome-wide maps of H3K4me2/3 in prostate cancer cell line LNCaP,数据在GEO可以下载。GSE20042,下面的所有分析,需要26G的空间。
作者想看看用 dihydrotestosterone (雄激素)处理了 cancer cell line LNCaP 这个细胞系之后,看看组蛋白甲基化修饰变化,主要是看H3K4me2和H3K4me3这两种组蛋白甲基化区别,分成三组,分别是处理前,处理后4H和16H,共有5个条件的数据,但是有7个fastq文件。