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

Original 生信技能树生信技能树 2022-06-06

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

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

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

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

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

下游Chip-seq数据分析

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

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

示例：

library(ChIPseeker)
library(ggplot2)
library(dplyr)
library(org.Hs.eg.db)

require(ChIPseeker)
f <- getSampleFiles()[[4]]#.bed文件,这里用包里的测试数据
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene#TxDb文件
#注释主体
x <- annotatePeak(f, tssRegion=c(-1000, 1000), TxDb=txdb)#x为输出对象

这里需要注意的是，启动子区域是没有明确的定义的，所以你可能需要指定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.查看输出文件

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

   > x
   Annotated peaks generated by ChIPseeker
   1331/1331  peaks were annotated
   Genomic Annotation Summary:
                Feature  Frequency
   9           Promoter 48.1592787
   4             5' UTR  0.7513148
   3             3' UTR  4.2073629
   1           1st Exon  0.7513148
   7         Other Exon  3.9068370
   2         1st Intron  3.6814425
   8       Other Intron  7.7385424
   6 Downstream (<=3kb)  1.1269722
   5  Distal Intergenic 29.6769346

as.GRanges(x) %>% head(5)
查看具体信息，会显示出每个peak的位置，所在的区域

1570289448158

1.2.2.输出文件解读

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

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

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

1.3.annotatePeak函数的参数

1.3.1.搜索peak上下游基因

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

addFlankGeneInfo=TRUE联合使用flankDistance=5000
peak位置-5k到5k的距离内所包含的基因
解析：输出中多三列: flank_txIds, flank_geneIds和flank_gene_distances，在指定范围内所有的基因都被列出，但是ID很不友好，看不出什么有用的信息，使用下面的参数会解决这个问题。

1.3.2.注释peak上下游基因

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

注释后会多出三列信息

ENSEMBL：不解释
SYMBOL：不解释
GENENAME：基因描述信息

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

sameStrand=TRUE

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

1.3.4.nearest相对整个基因

ignoreOverlap=TRUE

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

overlap=”TSS”

最近基因是相对于TSS

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.载入必要的包

library(ggplot2)
library(dplyr)
library(org.Hs.eg.db)
library(ChIPseeker)
require(TxDb.Hsapiens.UCSC.hg19.knownGene)

#Vennerable包没有镜像，需要从github下载安装
require(Vennerable, character.only = TRUE)
library(devtools)
install_github("js229/Vennerable")


txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

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

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

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

原始示例数据3：

1570329971545

把peaks数目只留下10行

1570330011350

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

1570330123157

#这里用的是包内置数据
library(ChIPseeker)
library(ggplot2)

files <- getSampleFiles()
#绘制第四个样本peaks在染色体上的分布
covplot(files[[4]])

#分别绘制不同样本上的peaks分布图
peak<-GenomicRanges::GRangesList(CBX6=readPeakFile(files[[4]]),
                                 CBX7=readPeakFile(files[[5]]))
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)

#具体到指定染色体上的指定范围
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7))+ facet_grid(chr ~ .id)

#covplot选项：
#chrs    指定想要展示染色体，默认为全部
#xlim    指定x轴坐标范围，即染色体的范围
#如果需要ggplot刻面展示，不同样本是用id来表示

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

1570330279933

2.3.固定窗口的peaks分布

library(ChIPseeker)
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
files <- getSampleFiles()
peak <- readPeakFile(files[[4]])
#选择promoter区作为窗口，使用转录起始位点，然后指定上下游，使用getPromoters函数，通过提取数据库里的数据，找到promoters区域上下游3000bp范围，准备好窗口
promoter <- getPromoters(TxDb=txdb, 
    upstream=3000, downstream=3000)
#使用getTagMatrix函数，把peak比对到这个窗口，并生成矩阵，供分析、可视化
tagMatrix <- getTagMatrix(peak,windows=promoter)

tagHeatmap(tagMatrix, xlim=c(-3000, 3000), 
    color="red")

#peakHeatmap函数只要给文件名，就可以直接出peaks的热图，是针对promoter区上下游的窗口进行绘图
peakHeatmap(files,TxDb=txdb,
            upstream=3000, downstream=3000, 
            color=rainbow(length(files)))

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

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

1570332880801

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

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

1570333070461

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

两种方式：

#方式一：直接从原始文件出相同的图
plotAvgProf2(files[[4]], TxDb=txdb, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency")

#方式二：
peak <- readPeakFile(files[[4]])
#选择promoter区作为窗口，使用转录起始位点，然后指定上下游，使用getPromoters函数，通过提取数据库里的数据，找到promoters区域上下游3000bp范围，准备好窗口
promoter <- getPromoters(TxDb=txdb, 
                         upstream=3000, downstream=3000)
#使用getTagMatrix函数，把peak比对到这个窗口，并生成矩阵，供分析、可视化
tagMatrix <- getTagMatrix(peak,windows=promoter)

plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
            xlab="Genomic Region (5'->3')", 
            ylab = "Read Count Frequency")

#添加置信区间
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), 
            conf = 0.95, resample = 1000)

#峰高证明结合的紧密

1570334187029

加上了置信区间：

1570334202334

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

2.5.注释信息的可视化

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

#整个对象
peakAnno <- annotatePeak(files[[4]], 
    tssRegion=c(-3000, 3000),
    TxDb=txdb, annoDb="org.Hs.eg.db")
#饼图
plotAnnoPie(peakAnno)
#柱状分布图
plotAnnoBar(peakAnno)
#vennpie
vennpie(peakAnno)
#UpSetR
upsetplot(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)#融合vennpie

plotAnnoPie

1570340148847

plotAnnoBar

1570340181116

vennpie

1570340197551

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

1570340214709

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

library(ChIPseeker)
library(clusterProfiler)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

bedfile=getSampleFiles()
seq=lapply(bedfile, readPeakFile)
genes=lapply(seq, function(i) 
  seq2gene(i, c(-1000, 3000), 3000, TxDb=txdb))

cc = compareCluster(geneClusters = genes, 
                    fun="enrichKEGG", organism="hsa")

dotplot(cc, showCategory=10)

4.自己结果的peaks注释

1570088544387

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

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

得到的文件有：

2.5K Oct  4 17:44 Xu_MUT_rep2_rmdup_summits.bed
2.4K Oct  4 17:20 Xu_MUT_rep2_sort_summits.bed
3.3K Oct  4 17:43 Xu_WT_rep1_rmdup_summits.bed
3.2K Oct  4 17:19 Xu_WT_rep1_sort_summits.bed

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

chr1    121485178   121485179   Xu_MUT_rep2_rmdup_peak_1    169.64561
chr10    42385559    42385560    Xu_MUT_rep2_rmdup_peak_2    28.26733
chr10   42387315    42387316    Xu_MUT_rep2_rmdup_peak_3    21.20116
chr10   42392730    42392731    Xu_MUT_rep2_rmdup_peak_4    24.90387
chr10   42599794    42599795    Xu_MUT_rep2_rmdup_peak_5    13.26007
chr11   18127886    18127887    Xu_MUT_rep2_rmdup_peak_6    5.11958
chr12   107104197    107104198   Xu_MUT_rep2_rmdup_peak_7    4.89626
chr16   46392027    46392028    Xu_MUT_rep2_rmdup_peak_8    18.40010
#~~~内容过多，无法全部展示~~~

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

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

library(ChIPseeker)
library(ggplot2)
library(dplyr)
library(org.Hs.eg.db)

rm(list = ls())
options(stringsAsFactors = F)

file=list.files("bedfile/")
filepath=file.path("bedfile",file)
filepath

require(ChIPseeker)
require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene#TxDb文件
require(clusterProfiler) 

##循环比对
for (bedPeaksFile in filepath) {
  peak=readPeakFile( bedPeaksFile )
  peakAnno=annotatePeak(peak, tssRegion=c(-3000, 3000), 
                           TxDb=txdb, annoDb="org.Hs.eg.db") 
  peakAnno_df=as.data.frame(peakAnno)
  write.table(peakAnno_df,sep = "\t",row.names = F,file = paste0(basename(bedPeaksFile),".xls"))
}

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

首先下载数据库：

perl ~/miniconda2/share/homer-4.10-0/configureHomer.pl -install hg19 &

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

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

利用homer进行寻找motif：

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

awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' Xu_MUT_rep2_rmdup_summits.bed >homer_peaks.tmp

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

mkdir /project/home/lyang/ChIP-seq_test/motif && cd /project/home/lyang/ChIP-seq_test/motif
ln -s ../6.macs2-callpeak/bedfile ./


for id in bedfile/*.bed;
do
echo $id
file=$(basename $id )
sample=${file%%.*} 
echo $sample

awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >homer_peaks.tmp  
findMotifsGenome.pl homer_peaks.tmp hg19 ${sample}_motifDir -p 10 -len 8,10,12 2>${sample}.motifrr.log 
annotatePeaks.pl    homer_peaks.tmp hg19  -p 10 1>${sample}.peakAnn.xls 2>${sample}.ann_err.log 
done


#Usage: findMotifsGenome.pl <input_pos_file> <genome> <output_directory> [options]
#-p：线程
#-len：motif length, default=8,10,12
#

#Usage: annotatePeaks.pl <input_peak_file | tss>  <genome version>  [options...]
#-p：线程
#

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

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

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

1570679919720

下图是R包的注释结果：

1570680164949

下面是homer的注释结果：

1570680374490

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

homer输出结果：

homer找motif结果文件：

1570347503028

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

1570347642069

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

一定要继续关注哦，下期更精彩！

学徒写在最后

首先感谢Jimmy老师的教程和代码，基本上只要跟着一步步学下来，肯定能复现漂亮的图，但是其中的原理需要自己仔细研究和领会。
另外，非常感谢jimmy老师对我耐心的指导和引导，当我遇到问题时，能一下了解我遇到的代码问题在哪里，比我百度谷歌半天教程都有用。

我写在后面

学徒已经做的很优秀了，一个月的时间总是短暂的，但学习的脚步不能停下，希望他回去以后能有更多的学习成果跟大家分享！

当然，如果你感兴趣学徒培养，也可以看看招学徒的通知：

你可以先看看我是如何培养学徒的：

高三女生醉酒后被强奸致死？检方回应

常德悲剧：让谴责无差别杀戮之声更加响亮一点

2024【公共营养师】培训报名通道已开启，不限学历，23岁及以上可报！还能领2000补贴

【惊】"以为要写遗书"! 飞温哥华航班遇炸弹气旋华人崩溃大哭连空姐都吐了; 客机颠簸盘旋3小时

女人最偏爱的十种男人

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

下游Chip-seq数据分析

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

1.1.ChIPseeker包输入文件

1.2.ChIPseeker包输出文件

1.2.1.查看输出文件

1.2.2.输出文件解读

1.3.annotatePeak函数的参数

1.3.1.搜索peak上下游基因

1.3.2.注释peak上下游基因

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

1.3.4.nearest相对整个基因

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

1.4.注释文件导出

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

2.1.载入必要的包

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

2.3.固定窗口的peaks分布

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

2.5.注释信息的可视化

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

4.自己结果的peaks注释

5.寻找motif

背景

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

首先下载数据库：

利用homer进行寻找motif：

学徒写在最后

您可能也对以下帖子感兴趣

高三女生醉酒后被强奸致死？检方回应

常德悲剧：让谴责无差别杀戮之声更加响亮一点

2024【公共营养师】培训报名通道已开启，不限学历，23岁及以上可报！还能领2000补贴

【惊】"以为要写遗书"! 飞温哥华航班遇炸弹气旋 华人崩溃大哭 连空姐都吐了; 客机颠簸盘旋3小时

女人最偏爱的十种男人

生成图片，分享到微信朋友圈

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

下游Chip-seq数据分析

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

1.1.ChIPseeker包输入文件

1.2.ChIPseeker包输出文件

1.2.1.查看输出文件

1.2.2.输出文件解读

1.3.annotatePeak函数的参数

1.3.1.搜索peak上下游基因

1.3.2.注释peak上下游基因

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

1.3.4.nearest相对整个基因

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

1.4.注释文件导出

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

2.1.载入必要的包

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

2.3.固定窗口的peaks分布

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

2.5.注释信息的可视化

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

4.自己结果的peaks注释

5.寻找motif

背景

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

首先下载数据库：

利用homer进行寻找motif：

学徒写在最后

您可能也对以下帖子感兴趣

【惊】"以为要写遗书"! 飞温哥华航班遇炸弹气旋华人崩溃大哭连空姐都吐了; 客机颠簸盘旋3小时