学徒第4周是ChIP-seq数据分析实战训练，讲义大纲文末的阅读原文，配套视频在B站：
因为公众号排版真的是力气活，同样我也是逼着学徒自己排版的，反正这辈子就苦这么一回！“作词作曲”都是他自己！
因为学徒提交作业实在是太长，超过了微信公众号发完限制，所以分成上下文！
上游数据分析在主页，现在我们一起来看看次页的下游分析！
示例：
这里需要注意的是，启动子区域是没有明确的定义的，所以你可能需要指定tssRegion，把基因起始转录位点的上下游区域来做为启动子区域。有了这两个输入（BED文件和TxDb对象），你就可以跑注释了，然后就可以出结果了。
.bed文件
由上游测序数据处理而来
TxDb文件
同名的包中含有相应的文件，直接引用即可，同样hg38也有同名的包，Bioconductor提供了30个TxDb包，可以供我们使用。要与测序数据比对一致
如果实在没有TxDb呢？
使用GenomicFeatures包函数来制作TxDb对象（这也是为什么说ChIPseeker支持所有物种）
示例：
makeTxDbFromUCSC：通过UCSC在线制作TxDb
makeTxDbFromBiomart: 通过ensembl在线制作TxDb
makeTxDbFromGRanges：通过GRanges对象制作TxDb
makeTxDbFromGFF：通过解析GFF文件制作TxDb
在R里打输出的对象，它会告诉我们ChIPseq的位点落在基因组上什么样的区域，分布情况如何。
查看具体信息，会显示出每个peak的位置，所在的区域
genomic annotation注释：annotation列，注释的是peak的位置，它落在什么地方了，可以是UTR，可以是内含子或者外显子
nearest gene annotation（最近基因）注释：多出来的其它列，注释的是peak相对于转录起始位点（TSS）的距离，不管这个peak是落在内含子或者别的什么位置上，即使它落在基因间区上，都能够注释到一个离它最近的基因（即使它可能非常远）
针对于不同的转录本，一个基因可能有多个转录起始位点，所以注释是在转录本的水平上进行的，可以看到输出有一列是transcriptId
两种注释策略面对不同的研究对象，第一种策略，peak所在的位置可能就是调控本身（比如你要做可变剪切的，最近基因的注释显然不是你关注的点）；而做基因表达调控的，当然promoter区域是重点，离结合位点最近的基因更有可能被调控。
想看peak上下游某个范围内（比如说-5k到5k的距离）都有什么基因：
addFlankGeneInfo=TRUE联合使用flankDistance=5000
peak位置-5k到5k的距离内所包含的基因
解析：输出中多三列: flank_txIds, flank_geneIds和flank_gene_distances，在指定范围内所有的基因都被列出，但是ID很不友好，看不出什么有用的信息，使用下面的参数会解决这个问题。
注释后会多出三列信息
ENSEMBL：不解释
SYMBOL：不解释
GENENAME：基因描述信息
sameStrand=TRUE
默认FALSE，ChIPseq数据通常情况下是没有正负链信息，但不排除特殊情况，可以给你的peak分别赋正负链，然后传入sameStrand=TRUE，分开做两次，你就可以分开拿到正链和负链的最近基因。总之，较少见
ignoreOverlap=TRUE
忽略当peak和TSS有overlap时，两种注释都指向同一个基因的冗杂信息，最近基因会去找不overlap的
overlap=”TSS”
最近基因是相对于TSS
overlap=”all”
最近基因是相对于整个基因
只想注释上游：ignoreDownstream = T
只想注释下游：ignoreUpstream = T
注：默认都是FALSE
covplot函数可以直接接受macs2-callpeak的输出bed文件进行画图。
函数返回的图中，每个竖线表示一个peaks在染色体的位置，线的宽度表示peaks在染色体的分布范围。
原始示例数据3：
把peaks数目只留下10行
修改其中9号染色体上peaks的范围，将其改大为chr9 2000 24381608。
不同样本的分布图（先用GenomicRanges::GRangesList构建对象）
经过修改示例文件3，只保留前30行，发现tagMatrix结果没有数据。说明了上述步骤是需要找到在promoter区域上下游3000bp范围内的peaks进行展示。
修改这个文件第一行chr11 114049934 114050346为chr11 1140 114050346后在进行同样的操作：
使用正常的示例文件4进行绘制热图展示：
两种方式：
加上了置信区间：
ChIPseeker还提供了getBioRegion函数，你可以指定’exon’或’intron’，它会像getPromoters一样为你准备好数据，即使这不是你感兴趣的，你也依然可以用它来准备窗口，然后可视化看一下。
主要就是注释一下，peaks峰落在的区域的分布，看看都落在了哪里
plotAnnoPie
plotAnnoBar
vennpie
upsetplot，还可以使用vennpie=T这个参数将上面的vennpie图放到upset图的左上角（但是我画的时候一直报错，说空间太小，把画布拉大了也还是不行，没找到原因）
经过前面的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包：
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数据库。
下载成功后会多出 ~/miniconda2/share/homer-4.10-0/data/genomes/hg19/ 文件夹, 共 5.9G
这个文件夹取决于你把homer这个软件安装到了什么地方。
因为homer对输入文件有格式要求，所以需要对call peaks的输出文件进行格式转换：
这里可以格式修改一起写到脚本中，利用循环在批量处理。
下图是R包的注释结果：
下面是homer的注释结果：
homer输出结果：
homer找motif结果文件：
homer注释结果文件（文件过长，找了几个比较短的可以看全，大致了解下输出文件的格式即可）：
一定要继续关注哦，下期更精彩！
首先感谢Jimmy老师的教程和代码，基本上只要跟着一步步学下来，肯定能复现漂亮的图，但是其中的原理需要自己仔细研究和领会。
另外，非常感谢jimmy老师对我耐心的指导和引导，当我遇到问题时，能一下了解我遇到的代码问题在哪里，比我百度谷歌半天教程都有用。
学徒已经做的很优秀了，一个月的时间总是短暂的，但学习的脚步不能停下，希望他回去以后能有更多的学习成果跟大家分享！