查看原文
其他

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

生信技能树 单细胞天地 2022-06-07

学徒第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(-10001000), 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.查看输出文件
  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
  1. as.GRanges(x) %>% head(5)

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


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

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(4e75e7))+ 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(-30003000), 
    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 114050346chr11 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(-30003000),
            xlab="Genomic Region (5'->3')"
            ylab = "Read Count Frequency")

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

#峰高证明结合的紧密

1570334187029

加上了置信区间:

1570334202334

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

2.5.注释信息的可视化

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

#整个对象
peakAnno <- annotatePeak(files[[4]], 
    tssRegion=c(-30003000),
    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(-10003000), 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(-30003000), 
                           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’motifRNA聚合酶结合位点的特异性序列。而且,当时的人们就知道,不是所有的结合位点都一定完美地与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老师对我耐心的指导和引导,当我遇到问题时,能一下了解我遇到的代码问题在哪里,比我百度谷歌半天教程都有用。

我写在后面


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


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

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

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

文章有问题?点此查看未经处理的缓存