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

生信技能树

发布于 2021-10-12 12:15:41

4.8K0

发布于 2021-10-12 12:15:41

文章被收录于专栏：生信技能树生信技能树

咱们《生信技能树》的B站有一个ChIP-

其中中国医科大的“小高”同学给大家带来的就是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-region和reference-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 文件。