九月学徒ChIP-seq学习成果展（6万字总结）（下篇）

学徒第4周是ChIP-seq数据分析实战训练，讲义大纲文末的阅读原文，配套视频在B站：

九月学徒已经结业，表现还不错，学了几个NGS组学数据处理加上部分单细胞，随机安排的文献数据处理图表复现也完成的还不赖，昨天在生信技能树的WGCNA代码就是他写的；重复一篇WGCNA分析的文章（代码版）

因为公众号排版真的是力气活，同样我也是逼着学徒自己排版的，反正这辈子就苦这么一回！“作词作曲”都是他自己！

因为学徒提交作业实在是太长，超过了微信公众号发完限制，所以分成上下文！

上游数据分析在主页，现在我们一起来看看次页的下游分析！

下游Chip-seq数据分析

1.Y叔的R包学习——注释

参考y叔的ChIP-seq数据分析大礼包

示例：

这里需要注意的是，启动子区域是没有明确的定义的，所以你可能需要指定tssRegion，把基因起始转录位点的上下游区域来做为启动子区域。
有了这两个输入（BED文件和TxDb对象），你就可以跑注释了，然后就可以出结果了。

1.1.ChIPseeker包输入文件

• .bed文件

• 由上游测序数据处理而来

• TxDb文件

• TxDb.Hsapiens.UCSC.hg19.knownGene

  同名的包中含有相应的文件，直接引用即可，同样hg38也有同名的包，Bioconductor提供了30个TxDb包，可以供我们使用。要与测序数据比对一致

• 如果实在没有TxDb呢？

  使用GenomicFeatures包函数来制作TxDb对象（这也是为什么说ChIPseeker支持所有物种）

  示例：

  require(GenomicFeatures)
hg19.refseq.db <- makeTxDbFromUCSC(genome="hg19", table="refGene")


• makeTxDbFromUCSC：通过UCSC在线制作TxDb

• makeTxDbFromBiomart: 通过ensembl在线制作TxDb

• makeTxDbFromGRanges：通过GRanges对象制作TxDb

• makeTxDbFromGFF：通过解析GFF文件制作TxDb

1.2.ChIPseeker包输出文件

1.2.1.查看输出文件

1. X

   在R里打输出的对象，它会告诉我们ChIPseq的位点落在基因组上什么样的区域，分布情况如何。

1. as.GRanges(x) %>% head(5)

   查看具体信息，会显示出每个peak的位置，所在的区域

1.2.2.输出文件解读

1. genomic annotation注释：annotation列，注释的是peak的位置，它落在什么地方了，可以是UTR，可以是内含子或者外显子

2. nearest gene annotation（最近基因）注释：多出来的其它列，注释的是peak相对于转录起始位点（TSS）的距离，不管这个peak是落在内含子或者别的什么位置上，即使它落在基因间区上，都能够注释到一个离它最近的基因（即使它可能非常远）

• 针对于不同的转录本，一个基因可能有多个转录起始位点，所以注释是在转录本的水平上进行的，可以看到输出有一列是transcriptId

1. 两种注释策略面对不同的研究对象，第一种策略，peak所在的位置可能就是调控本身（比如你要做可变剪切的，最近基因的注释显然不是你关注的点）；而做基因表达调控的，当然promoter区域是重点，离结合位点最近的基因更有可能被调控。

1.3.annotatePeak函数的参数

1.3.1.搜索peak上下游基因

想看peak上下游某个范围内（比如说-5k到5k的距离）都有什么基因：

1. addFlankGeneInfo=TRUE联合使用flankDistance=5000

   peak位置-5k到5k的距离内所包含的基因

2. 解析：输出中多三列: flank_txIds, flank_geneIds和flank_gene_distances，在指定范围内所有的基因都被列出，但是ID很不友好，看不出什么有用的信息，使用下面的参数会解决这个问题。

1.3.2.注释peak上下游基因

annoDb='org.Hs.eg.db'参数

注释后会多出三列信息

1. ENSEMBL：不解释

2. SYMBOL：不解释

3. GENENAME：基因描述信息

1.3.3.正负链分开注释（少用）

sameStrand=TRUE

默认FALSE，ChIPseq数据通常情况下是没有正负链信息，但不排除特殊情况，可以给你的peak分别赋正负链，然后传入sameStrand=TRUE，分开做两次，你就可以分开拿到正链和负链的最近基因。总之，较少见

1.3.4.nearest相对整个基因

1. ignoreOverlap=TRUE

忽略当peak和TSS有overlap时，两种注释都指向同一个基因的冗杂信息，最近基因会去找不overlap的

1. overlap=”TSS”

最近基因是相对于TSS

1. overlap=”all”

最近基因是相对于整个基因

1.3.5.只想注释上游或者是下游的基因

只想注释上游：ignoreDownstream = T

只想注释下游：ignoreUpstream = T

注：默认都是FALSE

1.4.注释文件导出

as.data.frame就可以转成data.frame，然后就可以用write.table

2.Y叔的R包学习——可视化

参考CS6: ChIPseeker的可视化方法

2.1.载入必要的包

2.2.观察peaks在染色体的分布

covplot函数可以直接接受macs2-callpeak的输出bed文件进行画图。

函数返回的图中，每个竖线表示一个peaks在染色体的位置，线的宽度表示peaks在染色体的分布范围。

原始示例数据3：

把peaks数目只留下10行

修改其中9号染色体上peaks的范围，将其改大为chr9    2000    24381608。

不同样本的分布图（先用GenomicRanges::GRangesList构建对象）

2.3.固定窗口的peaks分布

经过修改示例文件3，只保留前30行，发现tagMatrix结果没有数据。说明了上述步骤是需要找到在promoter区域上下游3000bp范围内的peaks进行展示。

修改这个文件第一行chr11 114049934 114050346为chr11 1140 114050346后在进行同样的操作：

可以认为由于修改了这个peak值的范围，所以找到了这个peak，但是因为范围比较大，就把在Promoter上下游3000bp范围内全部给覆盖了。如果有空，可以考虑先从TxDb.Hsapiens.UCSC.hg19.knownGene找到promoter的坐标，然后修改一个peak的范围再次进行展示实验。

使用正常的示例文件4进行绘制热图展示：

2.4.谱图的形式来展示结合的强度

两种方式：

加上了置信区间：

ChIPseeker还提供了getBioRegion函数，你可以指定’exon’或’intron’，它会像getPromoters一样为你准备好数据，即使这不是你感兴趣的，你也依然可以用它来准备窗口，然后可视化看一下。

2.5.注释信息的可视化

主要就是注释一下，peaks峰落在的区域的分布，看看都落在了哪里

plotAnnoPie

plotAnnoBar

vennpie

upsetplot，还可以使用vennpie=T这个参数将上面的vennpie图放到upset图的左上角（但是我画的时候一直报错，说空间太小，把画布拉大了也还是不行，没找到原因）

3.Y叔的R包学习——peak注释到的基因KEGG

4.自己结果的peaks注释

经过前面的ChIP-seq测序数据处理的常规分析，得到了BED格式的peaks记录文件，我选取的这篇文章里面做了2次ChIP-seq实验，分别是：

一个野生型MCF7细胞系的 BAF155 immuno precipitates和一个突变型MCF7细胞系的 BAF155 immuno precipitates，这样通过比较野生型和突变型MCF7细胞系的 BAF155 immuno precipitates的不同结果就知道该细胞系的BAF155 突变，对它在全基因组的结合功能的影响。

得到的文件有：

每个BED的peaks记录，本质是就3列是需要我们注意的，就是染色体，以及在该染色体上面的起始和终止坐标，如下，而peak calling的MACS输出bed文件默认是前5列，其中第五列指的是片段堆积的峰高：

我们所谓的peaks注释，就是想看看该peaks在基因组的哪一个区段，看看它们在各种基因组区域(基因上下游，5,3端UTR，启动子，内含子，外显子，基因间区域，microRNA区域)分布情况，但是一般的peaks都有近万个，所以需要批量注释。

这里我们使用Y叔的Chipseeker包：

5.寻找motif

背景

motif是比较有特征的短序列，会多次出现的，一般认为它的生物学意义重大，做完CHIP-seq分析之后，一般都会寻找motif 。

motif的定义如下：

motif: recurring pattern. eg, sequence motif, structure motif or network motif

DNA sequence motif: short, recurring patterns in DNA that are presumed to have a biological function.

从上边的定义可以看出，其实motif这个单词就是形容一种反复出现的模式，而序列motif往往是DNA上的反复出现的模式，并被假设拥有生物学功能。而且，经常是一些具有序列特异性的蛋白的结合位点（如，转录因子）或者是涉及到重要生物过程的（如，RNA 起始，RNA 终止， RNA 剪切等等）。

motif最先是通过实验的方法发现的，换句话说，不是说有了ChIP-seq才有了motif分析，起始很早人们就开始研究motif了！例如，'TATAAT’ box在1975年就被pribnow发现了，它与上游的‘TTGACA’motif是RNA聚合酶结合位点的特异性序列。而且，当时的人们就知道，不是所有的结合位点都一定完美地与motif匹配，大部分都只匹配了12个碱基中的7-9个。结合位点与motif的匹配程度往往也与蛋白质与DNA的结合强弱有关。

目前被人们识别出来的motif也越来越多，如TRANSFAC和JASPAR数据库都有着大量转录因子的motif。而随着ChIP-seq数据的大量产出，motif的研究会进一步深入，有一些课题组会将现有的ChIP-seq数据进行整合，提供更全面，更准确的motif数据库。

一篇文献列出了2014年以前的近乎所有知名的A survey of motif finding Web tools for detecting binding site motifs in ChIP-Seq data 链接见：https://biologydirect.biomedcentral.com/articles/10.1186/1745-6150-9-4

利用homer进行peaks的注释和寻找motif

首先下载数据库：

下载成功后会多出 ~/miniconda2/share/homer-4.10-0/data/genomes/hg19/ 文件夹, 共 5.9G

这个文件夹取决于你把homer这个软件安装到了什么地方。

利用homer进行寻找motif：

因为homer对输入文件有格式要求，所以需要对call peaks的输出文件进行格式转换：

这里可以格式修改一起写到脚本中，利用循环在批量处理。

把上面的代码保存为脚本runMotif.sh，然后运行：

nohup bash runMotif.sh 1>motif.log 2>runmotif_err.log &

不仅仅找了motif，还顺便把peaks注释了一下。得到的后缀为peakAnn.xls 的文件就可以看到和使用R包注释的结果是差不多的。

下图是R包的注释结果：

下面是homer的注释结果：

还可以使用meme来找motif，需要通过bed格式的peaks的坐标来获取fasta序列。MEME，链接：http://meme-suite.org/

homer输出结果：

homer找motif结果文件：

homer注释结果文件（文件过长，找了几个比较短的可以看全，大致了解下输出文件的格式即可）：

对了，昨天的九月学徒转录组学习成果展（3万字总结）（上篇） 也是我写的！