ChIP-seq数据分析课程学习笔记之peaks的可视化

咱们《生信技能树》的B站有一个ChIP-seq数据分析实战视频课程,缺乏配套笔记。恰好前些天的求贤令:曾经我给你带来了十万用户,但现在祝你倒闭,走大运结识了几位优秀小伙伴!

其中中国医科大的“小高”同学给大家带来的就是ChIP-seq数据分析实战视频课程的配套笔记,希望可以帮助大家更好的吸收消化课程内容!

首先视频免费共享在B站:【生信技能树】Chip-seq测序数据分析ChIP-SEQ实战演练的素材:链接:https://share.weiyun.com/53CwQ8B 密码:ju3rrh, 包括一些公司PPT,综述以及文献 ChIP-SEQ 实战演练的思维导图:文档链接:https://mubu.com/doc/11taEb9ZYg 密码:wk29

接下来是根据生信技能树课程、思维导图、PPT和综述文献整理的笔记。

了解bam、bedgraph、bw文件

wig、bigWig和bedgraph文件详解(http://www.bio-info-trainee.com/1815.html)

sam/bam文件:把测序reads比对到参考基因组后的文件;

bam或者bed文件:主要是为了追踪我们的reads到底比对到了参加基因组的什么区域;

UCSC规定的wig、bigWig和bedgraph文件:追踪参考基因组的各个区域的覆盖度,测序深度,还可以无缝连接到UCSC的Genome Browser工具里面进行可视化!

使用deeptools进行可视化

官方文档见http://deeptools.readthedocs.io/en/latest/index.html

以下引自deeptools辅助CHIP-seq数据分析-可视化(http://www.bio-info-trainee.com/2136.html)

第一个功能,把bam文件转换为bw格式文件:

bamCoverage -b tmp.sorted.bam -o tmp.bw

里面有一个参数非常重要,就是--extendReads 在 macs软件里面也有,macs2 pileup --extsize 200 ,就算是你的reads长度可能不一致,是否需要把它们补齐成一个统一的长度,因为我们只要是测到了reads,就代表那里是有signal的,只是因为我们的测序仪限制,测到的长度不够,或者质量不好,我们QC的时候,把前后碱基给trim掉了。还可以安装基因组的有效大小来对测序深度进行normlization。

第二个功能,画所有基因附近的信号热图,tools: computeMatrix, then plotHeatmap

http://deeptools.readthedocs.io/en/latest/content/example_step_by_step.html#heatmaps-and-summary-plots

需要自行下载合适的基因坐标记录文件,BED格式的。把上面两个命令结合起来即可,代码和图形实例如下:

第3个功能,画profile的图!

use computeMatrix for all of the signal files (bigWig format) at once

use plotProfile on the resulting file

首先把bam文件转为bw文件,详情:http://www.bio-info-trainee.com/1815.html

cd  ~/project/epi/mergeBam 
source activate chipseq
ls  *.bam  |xargs -i samtools index {} 
ls *.bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw & 
done

cd dup 
ls  *.bam  |xargs -i samtools index {} 
ls *.bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.rm.bw & 
done 

bw文件的IGV
找peaks

查看TSS附件信号强度(TSS 表示转录起始位点)

首先在UCSC的Table Browser 里面下载下面这个文件:

参考这篇教程:BED文件以及如何正确的从UCSC下载BED文件

## 首先对单一样本绘图: 
## both -R and -S can accept multiple files 
mkdir -p  ~/project/epi/tss 
cd  ~/project/epi/tss 
computeMatrix reference-point  --referencePoint TSS  -p 5  \
-b 10000 -a 10000 \
-R /home/data/vip13t16/project/epi/tss/ucsc.refseq.bed  \
-S /home/data/vip13t16/project/epi/mergeBam/H2Aub1.bw  \
--skipZeros  -o matrix1_test_TSS.gz  \
--outFileSortedRegions regions1_test_genes.bed

##  both plotHeatmap and plotProfile will use the output from  computeMatrix
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.png
plotHeatmap -m matrix1_test_TSS.gz  -out test_Heatmap.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.png
plotProfile -m matrix1_test_TSS.gz  -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720 

首先画10K附近

bed=/home/data/vip13t16/project/epi/tss/ucsc.refseq.bed
for id in  /home/data/vip13t16/project/epi/mergeBam/*bw ;
do 
echo $id
file=$(basename $id )
sample=${file%%.*} 
echo $sample

computeMatrix reference-point  --referencePoint TSS  -p 5  \
-b 10000 -a 10000 \
-R  $bed \
-S $id  \
--skipZeros  -o matrix1_${sample}_TSS_10K.gz  \
--outFileSortedRegions regions1_${sample}_TSS_10K.bed
# 输出的gz为文件用于plotHeatmap, plotProfile
##  both plotHeatmap and plotProfile will use the output from  computeMatrix
plotHeatmap -m matrix1_${sample}_TSS_10K.gz  -out ${sample}_Heatmap_10K.png
plotHeatmap -m matrix1_${sample}_TSS_10K.gz  -out ${sample}_Heatmap_10K.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_${sample}_TSS_10K.gz  -out ${sample}_Profile_10K.png
plotProfile -m matrix1_${sample}_TSS_10K.gz  -out ${sample}_Profile_10K.pdf --plotFileFormat pdf --perGroup --dpi 720

done 

使用命令批量提交:nohup bash 10k.sh 1>10k.log &

横坐标是起始位置,纵坐标是信号强度。

然后画2K的

bed=/home/data/vip13t16/project/epi/tss/ucsc.refseq.bed
for id in  /home/data/vip13t16/project/epi/mergeBam/*bw ;
do 
echo $id
file=$(basename $id )
sample=${file%%.*} 
echo $sample

computeMatrix reference-point  --referencePoint TSS  -p 15  \
-b 2000 -a 2000 \
-R  $bed \
-S $id  \
--skipZeros  -o matrix1_${sample}_TSS_2K.gz  \
--outFileSortedRegions regions1_${sample}_TSS_2K.bed

##  both plotHeatmap and plotProfile will use the output from  computeMatrix
plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.png
plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.png
plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720

done 

使用命令批量提交:nohup bash 2k.sh 1>2k.log &

Control在TSS附近没有的信号。

还可以给多个bed文件来绘图,还可以画genebody的图。

上面的批量代码其实就是为了统计全基因组范围的peak在基因特征的分布情况,也就是需要用到computeMatrix计算,用plotHeatmap热图的方式对覆盖进行可视化,用plotProfile折线图的方式展示覆盖情况。

computeMatrix具有两个模式:scale-regionreference-point。前者用来信号在一个区域内分布,后者查看信号相对于某一个点的分布情况。无论是那个模式,都有有两个参数是必须的,-S是提供bigwig文件,-R是提供基因的注释信息。还有更多个性化的可视化选项。

使用R包对找到的peaks文件进行注释

还需要安装必备R包:

options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos"=c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
install.packages("devtools")
library(devtools) 
BiocManager::install
BiocManager::install(c('airway','DESeq2','edgeR','limma'))
BiocManager::install(c('ChIPpeakAnno','ChIPseeker'))
BiocManager::install('TxDb.Hsapiens.UCSC.hg19.knownGene',
ask=F,suppressUpdates=T)
BiocManager::install('TxDb.Hsapiens.UCSC.hg38.knownGene',
ask=F,suppressUpdates=T)
BiocManager::install('TxDb.Mmusculus.UCSC.mm10.knownGene',
ask=F,suppressUpdates=T)  
library(TxDb.Hsapiens.UCSC.hg19.knownGene) 
library(TxDb.Mmusculus.UCSC.mm10.knownGene) 
library(TxDb.Hsapiens.UCSC.hg38.knownGene) 
library(ChIPpeakAnno) 
library(ChIPseeker)

bedPeaksFile   = 'H2Aub1_summits.bed'; 
bedPeaksFile
## loading packages
require(ChIPseeker)
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
require(clusterProfiler) 
peak <- readPeakFile( bedPeaksFile )  
keepChr= !grepl('_',seqlevels(peak))
seqlevels(peak, pruning.mode="coarse") <- seqlevels(peak)[keepChr]
peakAnno <- annotatePeak(peak, tssRegion=c(-3000, 3000), 
TxDb=txdb, annoDb="org.Mm.eg.db") 
peakAnno_df <- as.data.frame(peakAnno)

seqnames start end width strand V4 V5 annotation geneChr geneStart geneEnd geneLength geneStrand geneId transcriptId distanceToTSS ENSEMBL SYMBOL GENENAME
1 chr1 4493147 4493147 1 * H2Aub1_peak_1 9.12358 Promoter (<=1kb) 1 4492465 4493735 1271 2 20671 ENSMUST00000191939.1 588 ENSMUSG00000025902 Sox17 SRY (sex determining region Y)-box 17
2 chr1 4496501 4496501 1 * H2Aub1_peak_2 5.41762 Promoter (<=1kb) 1 4490931 4496413 5483 2 20671 ENSMUST00000027035.9 -88 ENSMUSG00000025902 Sox17 SRY (sex determining region Y)-box 17
3 chr1 4497427 4497427 1 * H2Aub1_peak_3 7.23122 Promoter (<=1kb) 1 4491390 4497354 5965 2 20671 ENSMUST00000192650.5 -73 ENSMUSG00000025902 Sox17 SRY (sex determining region Y)-box 17
4 chr1 5916563 5916563 1 * H2Aub1_peak_4 8.73105 Promoter (<=1kb) 1 5913707 5917398 3692 2 226304 ENSMUST00000044180.4 835 ENSMUSG00000033774 Npbwr1 neuropeptides B/W receptor 1
5 chr1 5917014 5917014 1 * H2Aub1_peak_5 10.94250 Promoter (<=1kb) 1 5913707 5917398 3692 2 226304 ENSMUST00000044180.4 384 ENSMUSG00000033774 Npbwr1 neuropeptides B/W receptor 1
6 chr1 17097438 17097438 1 * H2Aub1_peak_6 6.16313 Promoter (<=1kb) 1 16994938 17097889 102952 2 57339 ENSMUST00000038382.4 451 ENSMUSG00000042686 Jph1 junctophilin 1
7 chr1 18058179 18058179 1 * H2Aub1_peak_7 4.51036 Intron (ENSMUST00000209405.1/ENSMUST00000209405.1, intron 1 of 3) 1 18115204 18137051 21848 2 78081 ENSMUST00000115340.7 78872 ENSMUSG00000025774 Crisp4 cysteine-rich secretory protein 4
8 chr1 19210429 19210429 1 * H2Aub1_peak_8 7.67158 Promoter (1-2kb) 1 19208960 19238083 29124 1 21419 ENSMUST00000187754.6 1469 ENSMUSG00000025927 Tfap2b transcription factor AP-2 beta
9 chr1 19211915 19211915 1 * H2Aub1_peak_9 5.20813 Promoter (<=1kb) 1 19212054 19234624 22571 1 21419 ENSMUST00000064976.5 -139 ENSMUSG00000025927 Tfap2b transcription factor AP-2 beta

#一些可视化
plotAnnoPie (peakAnno)

plotAnnoBar(peakAnno)

纵坐标为样本,横坐标为各功能性区域中 Peak 的比例。

Promoter:启动子到转录起始位点的区域中的 Peak 的比例;

5’UTR:5’ 非编码区域中Peak的比例;

3’UTR:3’ 非编码区域中Peak的比例;

Exon:外显子中的 Peak 的比例;

Intron:内含子中的 Peak 的比例;

Intergenic:基因间区中的 Peak 的比例。

可以载入IGV看看效果,检测软件找到的peaks是否真的合理,还可以配合rmarkdown来出自动化报告。

也可以使用其它代码进行下游分析;https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq

homer软件来寻找motif

转录因子往往倾向于结合在特定的 DNA 序列上,这些 DNA 序列通常具有高度相似的核苷酸序列模式,即每个转录因子都有一个目标 DNA 序列的 Motif,公开发表的转录因子数据库中一般记录了转录因子对应的 Motif 信息。

安装:conda install -c bioconda homer,  找到自己安装的homer,然后使用其附带的配置脚本来下载数据库。

可以参考这个教程:Homer软件的介绍-最全面而详细的找motif教程

cd ~/biosoft
mkdir homer &&  cd homer
wget http://homer.ucsd.edu/homer/configureHomer.pl 
perl configureHomer.pl -install
perl configureHomer.pl -install mm10

homer软件找motif整合了两个方法,包括依赖于数据库的查询,和de novo的推断,都是读取ChIP-seq数据上游分析得到的bed格式的peaks文件。

运行homer软件

但是使用起来很简单:http://homer.ucsd.edu/homer/ngs/peakMotifs.html

cd  ~/project/epi/motif  
for id in /home/data/vip13t16/project/epi/peaks/*.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 mm10 ${sample}_motifDir -len 8,10,12
annotatePeaks.pl homer_peaks.tmp mm10  1>${sample}.peakAnn.xls 2>${sample}.annLog.txt 
done 

把上面的代码保存为脚本runMotif.sh,然后运行:nohup bash runMotif.sh 1>motif.log &

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

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

HOMER结果

每个样本的分析结果保存在后缀为 "_MotifOutput" 的文件夹路径下,每个文件夹下均有已知 Moitf 和未知 Motif 两种分析结果

HOMER结果

knownResults.html为在公共数据库中已有的 Motif 的计算结果汇总表

![Homer Known Motif Enrichment Results (H3K36me3_summits_motifDir)](../Library/Application Support/typora-user-images/image-20210829125310823.png)

Rank:根据 P-value 对 Motif 进行排序得到的排名; Motif:该 Motif 的序列; Name: 该 Motif 对应的转录因子,包含已知 Motif 的数据来源(HOMER / Jaspar 等); P-value:该 Motif 的 pvalue 值,p 值越小,预测的转录因子可信度越高; Log P_value:P-value 以常数 e 为底取对数; q-value: 该 Motif 的 qvalue 值,q 值越小,预测的转录因子可信度越高; Target Sequence with motif:预测的转录因子在给定 DNA 序列(HOMER 分析的时候提供的序列信息)上的结合位点数目; % of Target Sequence with motif:预测的转录因子在给定 DNA 序列(HOMER分析的时候提供的序列信息)上的结合位点与所有转录因子结合位点的占比值; Background Sequence with Motif:预测的转录因子在背景序列(基因组 DNA 序列)上的结合位点数目; % of Background Sequence with Motif:预测的转录因子在背景序列(基因组 DNA 序列)上的结合位点与所有转录因子结合位点的占比值; Motif File:Motif 的 matrix 文件;SVG:Motif 序列信息的 SVG 格式文件。

homerResults.html 为未知 Motif 的计算结果汇总表,homerResults 文件夹中存放的是该项结果的分析数据和图片。

![Homer de novo Motif Results (H3K36me3_summits_motifDir/)](../Library/Application Support/typora-user-images/image-20210829125736598.png)

Rank:根据 P-value 对 Motif 进行排序得到的排名; Motif:该 Motif 的序列; P-value:该 Motif 的 pvalue 值,p 值越小,预测的转录因子可信度越高; Log P_value:P-value 以常数 e 为底取对数; % of Targets:该 Motif 在给定 DNA 序列(HOMER 分析的时候提供的序列信息)上的结合位点与所有未知 Motif 结合位点的占比值; % of Background:该 Motif 在背景序列(基因组 DNA 序列)上的结合位点与所有未知 Motif 结合位点的占比值; Best Match/Details:与该 Motif 最相似的已知转录因子 Motif; Motif File:该 Motif 的 matrix 文件。

其它高级分析

比如可以 比较不同的peaks文件,代码见:https://github.com/jmzeng1314/NGS-pipeline/blob/master/CHIPseq/step6-ChIPpeakAnno-Venn.R

大致流程

以下是这篇文章 Chip-seq 部分的大致流程, 这部分 NGS 实战系列差不多就结束啦!

(0)

相关推荐