CHIP-seq

FROM (https://www.cnblogs.com/fengchuiguo1994/articles/6966089.html)

1:ChIP-Seq数据是基因组特异性富集的序列的测序结果,包括组蛋白修饰ChIP-Seq(H3K4me3/启动子相关/narrowpeak、H3K4me1/增强子相关/narrowpeak、H3K27ac/增强子相关/broadpeak)、转录因子ChIP-Seq(CTCF/绝缘子相关/narrowpeak、pol II/转录起始/narrowpeak)、DNA富集序列(DNase-Seq/弱DNA酶消化/活性区域、MNase-Seq/强DNA酶消化/核小体不活跃区域、ATAC-Seq//前两者的结果的集合)。

通过互补染色质分析实验分析的基因组位点揭示了染色质结构的不同方面:ChIP-seq显示特异性转录因子(TF)的结合位点; DNase-seq,ATAC-seq和FAIRE-seq显示开放染色质的区域;和MNase-seq鉴定良好定位的核小体。在ChIP-seq中,特异性抗体用于直接或通过包含靶因子的复合物中的其他蛋白质提取结合至靶蛋白的DNA片段。在DNase-seq中,染色质被DNA酶I内切核酸酶轻微消化。大小选择用于富集在DNA对DNA酶I攻击高度敏感的染色质区域产生的片段(在初期会生成包含各种长度的DNA小片段,但是一般来书保留100~300bp长度的小片段建库测序)。 ATAC-seq是DNase-seq的替代方法,其使用工程改造的Tn5转座酶来切割DNA并将引物DNA序列整合到切割的基因组DNA中(即,标记)。微球菌核酸酶(MNase)是内切核酸外切酶,其连续地消化DNA直到达到阻塞(和DNA酶相比(DNase-seq),属于强切,开放的区域全部都被消化),例如核小体。在FAIRE-seq中,甲醛用于交联染色质,并且苯酚 - 氯仿用于分离剪切的DNA。

详细介绍可参考(http://www.nature.com/nrg/journal/v15/n11/fig_tab/nrg3788_F1.html)

image.png

2:ChIP-Seq数据的作用:a:构建物种的epigenome,利用chromHMM将基因组分成一个一个的区域;b:与交互数据(HiC/chia-pet)联合分析;c:和RNA-Seq联合分析(chirp-seq)。

3:ChIP-Seq数据的分析流程:

1:预处理,步骤与RNA-Seq的一致,详情见RNA-Seq分析的1、2两步。

2:比对:DNA数据比对软件用的比较多的是bwa和bowtie,bwa比对结果相较于bowtie的比对结果更加准确,但是跑得慢,bowtie正好相反。我一般是用bwa比对,很少使用bowtie比对,慢不了多少。而且方便下游分析(比如使用GATK call snv/indel,bwa的结果明显比bowtie的要好)。

bwa比对提供了两种主流的比对aln和mem,一般序列长度小于70bp选用aln比对,大于70bp的时候选用mem。

首先第一步都是根据参考基因组建索引:bwa index genome.fa

当长度小于70bp时:

针对于单端数据:bwa aln -t 4 -f file.sai genome.fa file.fastq

bwa samse -n 1 -f file.sam genome.fa file.sai file.fastq

针对于双端数据:bwa aln -t 4 -f file1.sai genome.fa file1.fastq

bwa aln -t 4 -f file2.sai genome.fa file2.fastq

bwa sampe -n 1 -f file.sam genome.fa file1.sai file1.fastq file2.sai file2.fastq

-t 线程数目  -f 输出文件名  -n 最大输出条目(说白了就是多位点比对,为了省事可以考虑设置为1,就只输出一个比对的位点,一般不超过3)

可加参数 -r ,头文件是为样本加入分组信息。为了下游分析的正常进行,我一般不在此处加,在后续需要的时候用picard软件加上头文件信息来分组

当长度大于70bp时:

针对于单端数据:bwa mem -t 4 genome.fa file.fastq > file.sam

针对于双端数据:bwa mem -t 4 genome.fa file1.fastq file2.fastq > file.sam

针对于mem比对一般不考虑加入其他的参数。

可加参数 -r ,头文件是为样本加入分组信息。为了下游分析的正常进行,我一般不在此处加,在后续需要的时候用picard软件加上头文件信息来分组

3:sam2bam:相对来说bam格式作为sam文件二进制的形式,存在占用内存更小、软件读取效率更高等优点,所以在此处将sam文件转成bam文件

sam文件转bam文件并排序:samtools view -Su file.sam | samtools sort -o file.bam -(samtools view -q 20 -Su file.sam | samtools sort -o file.bam - 或者可以直接在该步骤完成质量过滤但不建议)。

去除duplication:java -jar ~/Tools/picard-tools-1.119/MarkDuplicates.jar I=file.bam O=file.rmd.bam CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=file.metrics

可以加入的参数:REMOVE_DUPLICATES=true。删除duplication的reads而不是修改该reads的flag(第二列的值)。但是为了保证数据的完整性不建议这么做。

提取高质量的比对结果:samtools view -F 1024 -q 20 -o file.flt.bam file.bam

此时的file.flt.bam文件就可以导入到可视化软件(IGV)中查看了,而不应该导入file.bam,因为file.bam中包含了很多的duplication和低比对的reads,导致和最终的结果有差异,但是查看file.flt.bam就不会。此外,file.bam主要是为了后续出什么问题方便来检查,而且可以作为其他分析的输入文件,最重要的是bam文件的大小是小于原始数据(fastq)的,可以用来保存原始数据,在需要的时候可以通过bedtools将ban文件转成fastq文件。

额外的分析:在该步骤的过程中可以加入样本重复性的分析,deeptools中的multiBamSummary和plotCorrelation可以用来检查样本的重复性

4:peak calling:该步骤是ChIP-Seq中最重要的一个步骤,做的好坏直接影响最终的结果,在peak calling时候重要有两个软件:macs2和F-seq。一般来说macs2用的比较多,针对于各种不同数据之间的参数选择:a:转录因子和组蛋白修饰中的H3K4me3、H3K4me3直接用默认的参数运行,生成narrowpeak文件;b:除了H3K4me3和H3K4me1以外的组蛋白修饰的ChIP-Seq 加上 -broad 参数即可;c:针对DNase-Seq数据加上参数 --nomodel --shift 100 --extsize 200;针对MNase-Seq数据加上参数 --nomodel --shift 37 --extsize 73。结果文件说明见官网

5:peak annotation:获得peak文件后能做的分析就比较多了,其中的一个就是对peak文件的注释。注释软件一般用homers和chipseeker,不过后来仔细比较了两个软件注释的结果,发现homer注释的结果不是很正确,相反chipseeker要好一些。chipseeker是一个R包,需要R3.3.0以上的版本,需要安装clusterProfiler和GenomicFeatures包。

注释的代码见下(chipseeker还可以做很多事,可以去官方网站查看):

library(ChIPseeker)
library(clusterProfiler)
library(GenomicFeatures)
files <- list(randA=c("/home/xyhuang/program/randomA20.peak"),randB=c("/home/xyhuang/program/randomB20.peak"))
txdb <- makeTxDbFromGFF("hg19.gtf",format="gtf")
peakAnno <- annotatePeak(files[[1]], tssRegion=c(-3000, 3000), TxDb=txdb)
pdf("randA.pdf")
plotAnnoPie(peakAnno)
dev.off()
data <- as.data.frame(peakAnno)
aa <- c(7,8,9,10,11,12)
data <- data[-aa]
write.table(data,"randA.txt",row.names=F,quote = F,sep="\t")

peakAnno <- annotatePeak(files[[2]], tssRegion=c(-3000, 3000), TxDb=txdb)
pdf("randB.pdf")
plotAnnoPie(peakAnno)
dev.off()
data <- as.data.frame(peakAnno)
aa <- c(7,8,9,10,11,12)
data <- data[-aa]
write.table(data,"randB.txt",row.names=F,quote = F,sep="\t")

6:检查peak是否富集:使用bedtools random 随机生成一组peak,然后重复步骤五,比较两个结果是否有差异,差异越大表示富集的程度越高越好。

samtools fadix genome.fa

awk '{print 1"\t"2;}' genome.fa.fai > genome.size

bedtools random -l 100 -n 100000 -g genome.size

-l peak的宽度  -n 生成的peak数目

7:motif分析:针对与转录因子和组蛋白修饰,目前认为是蛋白结构特异性的识别基因组上的一段特异的序列,而这段序列是具有保守性的。可以用来鉴定基因组上的调控元件

motif分析主要也是两个软件homer和meme,根据注释结果,我认为meme寻找的结果更加可靠,而且meme专门为ChIP-Seq数据开发了一个meme-chip软件包。但是考虑到需要将寻找到的motif与motif库进行比较,所以一般推荐网页上做。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 206,723评论 6 481
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,485评论 2 382
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 152,998评论 0 344
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,323评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,355评论 5 374
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,079评论 1 285
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,389评论 3 400
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,019评论 0 259
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,519评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,971评论 2 325
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,100评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,738评论 4 324
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,293评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,289评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,517评论 1 262
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,547评论 2 354
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,834评论 2 345

推荐阅读更多精彩内容