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
查看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 两种分析结果
knownResults.html为在公共数据库中已有的 Motif 的计算结果汇总表
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 文件夹中存放的是该项结果的分析数据和图片。
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 实战系列差不多就结束啦！