Chip-seq上游分析流程学习(四)

本次分析步骤包括：环境部署——数据下载——查看数据(非过滤)——数据质控清洗——数据比对——去除pcr重复——peaks calling——可视化

IGV可视化

首先确定基因组，如果需要研究其他物种的则需要自行下载从File列表中load from File载入bed文件可以输入不同的基因，比如TP53进行查看~

1. 参考基因组：Human(GRCh38/hg38) 表示这些数据是基于人类基因组参考序列GRCh38/hg38版本的。
2. 染色体位置：列出了染色体17号（chr17）上的不同区域，如p13.2、p13.1、p12、p11.2、p11.1等，这些是染色体的特定区域。
3. 基因名称：TP53 是一个著名的基因，与癌症相关，这里列出了它在染色体17号上的位置。
4. 基因组坐标：chr17:7,666,421-7,689,490 是基因TP53在染色体17号上的精确位置，以碱基对（base pairs）为单位。
5. 基因组大小：7,668 kb、7,670 kb、7,672 kb 等表示的是基因组的大致大小，单位是千碱基对（kilobase）。
6. p300 intersect peak：p300 intersect peak 14438 可能表示p300蛋白与DNA结合的峰值区域，编号为14438。

deeptools可视化

多样本分析

使用multiBamSummary、plotCorrelation和plotPCA三个模块。主要查看对照组和处理组中的组间差异和组内相似性(如果上一步把BAM转换成BW, 那么multiBamSummary可以用multiBigWigSummary替代)

path="/home/lm/Z_Projects/chipseq"
mkdir deeptools

# 统计reads在全基因组范围的情况
# 双端测序
#multiBamSummary bins -bs 1000 --bamfiles ${path}/intersect/*_sort.bam --extendReads 130 -out treat_results.npz
# 单端测序
nohup multiBamSummary bins -bs 1000 --bamfiles ${path}/intersect/*_sort.bam -out ./deeptools/treat_results.npz > multibam.log &

# 相关性散点图
plotCorrelation -in ./deeptools/treat_results.npz -o ./deeptools/treat_results.pdf --corMethod spearman -p scatterplot

# 热图
plotCorrelation -in ./deeptools/treat_results.npz -o ./deeptools/treat_results_heatmap.pdf --corMethod spearman -p heatmap

# 主成分分析
plotPCA -in ./deeptools/treat_results.npz  -o ./deeptools/pca.pdf

1. multiBamSummary: 这是 deepTools 提供的一个工具，用于计算多个 BAM 文件的 reads 覆盖情况。它可以帮助比较不同样本在基因组上的相似性，通常用于质量控制和样本间的一致性分析。bins: 这个参数指定 multiBamSummary 使用的是“bins”模式。在这种模式下，基因组会被分成等大小的区块（称为 bins），然后统计每个 bin 中 reads 的覆盖情况。这种方式适合评估不同 BAM 文件之间的全基因组相似性。-bs 1000: -bs 是 --binSize 的缩写，表示 bin 的大小。在这里设置为 1000，表示将基因组划分为每个大小为 1000 bp（碱基对）的 bins，即每个 bin 代表 1 kb 的基因组区域。
2. --bamfiles {path}/intersect/*_sort.bam 通过变量 ${path} 来指定 BAM 文件所在的路径，并使用通配符 *_sort.bam 来匹配所有以 _sort.bam 结尾的 BAM 文件。这样可以同时分析多个样本的数据。这些 BAM 文件是经过排序的 BAM 文件，它们包含了 reads 在基因组上的比对信息。
3. -out treat_results.npz: -out 用于指定输出文件的名称。treat_results.npz 是输出文件的名称，.npz 是 numpy 的压缩文件格式。这种格式的文件可以直接用于 deepTools 的其他工具，例如 plotCorrelation 和 plotPCA，用于进一步的相关性分析和主成分分析。

peak分布可视化

打开UCSC点击Table browser下载不同物种所有基因的转录起始位点(TSS)区域的bed文件选择默认参数，红色方框处要与上图一致

处理单一样本，10K上下游区域

# 把上面得到的结果放进TSS文件夹中
mkdir TSS

# 设定路径
path="/home/lm/Z_Projects/chipseq"

## 处理单一样本
## both -R and -S can accept multiple files

computeMatrix reference-point --referencePoint TSS -p 15 \
-b 10000 -a 10000 \
-R ${path}/TSS/hg38_TSS.bed \
-S ${path}/bamCoverage/p300.bw \
--skipZeros -o ./TSS/p300_TSS.gz \
--outFileSortedRegions ./TSS/p300_genes.bed

## both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m ${path}/TSS/p300_TSS.gz -out ${path}/TSS/test.png
plotHeatmap -m ${path}/TSS/p300_TSS.gz -out ${path}/TSS/test.pdf --plotFileFormat pdf --dpi 720
plotProfile -m ${path}/TSS/p300_TSS.gz -out ${path}/TSS/test.png
plotProfile -m ${path}/TSS/p300_TSS.gz -out ${path}/TSS/test.pdf --plotFileFormat pdf --perGroup --dpi 720

computeMatrix reference-point 命令

这个命令用于生成一个矩阵文件，该文件包含了以转录起始位点（TSS）为中心的 p300 蛋白结合信号的分布情况。参数解释如下：

1. --referencePoint TSS: 指定参考点为转录起始位点。
2. -p 15: 使用 15 个处理器并行计算。
3. -b 10000: 指定参考点上游的区域长度为 10000 个碱基对。
4. -a 10000: 指定参考点下游的区域长度为 10000 个碱基对。
5. -R ${path}/TSS/hg38_TSS.bed: 指定包含转录起始位点的 BED 文件。
6. -S ${path}/p300.bw: 指定包含 p300 蛋白结合信号的 bigWig 文件。
7. --skipZeros: 跳过信号值为零的区域。
8. -o p300_TSS.gz: 指定输出矩阵文件的名称。
9. --outFileSortedRegions p300_genes.bed: 指定输出排序后的区域列表文件的名称。

plotHeatmap 命令

这个命令用于生成热图，展示基因组区域的信号强度。参数解释如下：

1. -m p300_TSS.gz: 指定输入的矩阵文件。
2. -out test.png: 指定输出的热图文件名，格式为 PNG。
3. --plotFileFormat pdf: 指定输出文件的格式为 PDF。
4. --dpi 720: 指定输出图像的分辨率为 720 DPI。

plotProfile 命令

这个命令用于生成信号的平均分布图，展示基因组区域的信号强度随位置的变化。参数解释如下：

1. -m p300_TSS.gz: 指定输入的矩阵文件。
2. -out test.png: 指定输出的轮廓图文件名，格式为 PNG。
3. --plotFileFormat pdf: 指定输出文件的格式为 PDF。
4. --perGroup: 为每个组（如每个基因或每个样本）生成单独的轮廓图。
5. --dpi 720: 指定输出图像的分辨率为 720 DPI。

绘制10k区域

# 把上面得到的结果放进TSS文件夹中
mkdir TSS

### 如果要批量处理
path="/home/lm/Z_Projects/chipseq"
bed="/home/lm/Z_Projects/chipseq/TSS/hg38_TSS.bed"

for id in ${path}/bamCoverage/*.bw ;
do
    echo $id
    file=$(basename $id)
    sample=${file%.*}
    echo $sample

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

    ## both plotHeatmap and plotProfile will use the output from computeMatrix
    plotHeatmap -m ${path}/TSS/matrix1_${sample}_TSS_10K.gz -out ${path}/TSS/${sample}_Heatmap_10K.png
    plotHeatmap -m ${path}/TSS/matrix1_${sample}_TSS_10K.gz -out ${path}/TSS/${sample}_Heatmap_10K.pdf --plotFileFormat pdf --dpi 720
    plotProfile -m ${path}/TSS/matrix1_${sample}_TSS_10K.gz -out ${path}/TSS/${sample}_Profile_10K.png
    plotProfile -m ${path}/TSS/matrix1_${sample}_TSS_10K.gz -out ${path}/TSS/${sample}_Profile_10K.pdf --plotFileFormat pdf --perGroup --dpi 720
done > 10k.log 2>&1 &

1. for id in ${path}/bamCoverage/*.bw: 遍历 bamCoverage 文件夹下的所有 .bw 文件。每个文件路径赋值给变量 id。
2. file=id): 提取文件的名称（不包含路径），如 /path/to/BRD4.bw 会变成 BRD4.bw。
3. sample=${file%.*}: 提取文件名的前缀（去掉扩展名），如 BRD4.bw 会变成 BRD4。这个 sample 会用作后续生成文件的名称前缀。

BRD4_Profile_10K的结果

1. 横轴：表示距离 TSS 的距离（单位为 bp），与折线图一致，从 -10Kb 到 +10Kb。
2. 纵轴：表示每个目标基因（或基因组区域）的排序。每一行代表一个基因。
3. 颜色：表示信号强度，颜色越亮信号越强（右侧的颜色条表示强度从 0 到 10 的变化）。
4. 热图中心（TSS 附近）信号强度最强（颜色最亮），进一步支持 BRD4 在 TSS 附近高度富集。
5. 各行之间的信号强度分布显示了一定的差异性，说明 BRD4 的结合可能因基因而异。
6. 热图底部或顶部的一些行信号较弱，可能是这些基因的 BRD4 结合活性较低。截取了一部分，图片实在太长了

绘制2k区域

### ### 如果要批量处理
bed=/public/annotation/CHIPseq/mm10/ucsc.refseq.bed
for id in /home/jmzeng/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

R语言注释peak

# 设置 CRAN 镜像为清华大学
options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# 设置 BioC_mirror 镜像为西湖大学
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")

# 安装R包
BiocManager::install(c("airway", "DESeq2", "edgeR", "limma", "ChIPpeakAnno", "ChIPseeker"))
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene")

1.导入

rm(list = ls())
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(ChIPpeakAnno)
library(ChIPseeker)

2.数据预处理

#以人类基因组 hg38 版本为例,通过GenomicFeatures包构建构建TxDb对象
# library(GenomicFeatures)
# genome <- makeTxDbFromGFF('Homo_sapiens.GRCh38.113.chr.gff3')
#除 makeTxDbFromGFF 外，还可以使用以下方法构建 TxDb 对象，
#makeTxDbFromUCSC，通过 UCSC 构建 TxDb
#akeTxDbFromBiomart，通过 ensembl 构建 TxDb
#makeTxDbFromGRanges，通过 GRanges 构建 TxDb

# 指定参考基因组
require(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene 
txdb

#读取 bed 文件
peak <- readPeakFile('BRD4_intersect_summits.bed')
keepChr <- !grepl("_",seqlevels(peak)) #去除带有“_”的染色体
seqlevels(peak,pruning.mode="coarse") <- seqlevels(peak)[keepChr]
#peaks 的基因组注释
peakAnno <- annotatePeak(peak, tssRegion = c(-3000, 3000), 
                         TxDb = txdb,annoDb = "org.Hs.eg.db")
#输出结果
write.table(peakAnno, 'peak.txt', sep = '\t', row.names = FALSE, quote = FALSE)
peakAnno_df <- as.data.frame(peakAnno)

3.可视化

# 查看所有peaks在染色体的分布情况
covplot(peak,weightCol = "V5")
# 查看peaks在所有基因启动子附近的分布情况
promoter <- getPromoters(TxDb = txdb,upstream = 3000, downstream = 3000)
tagMatrix <- getTagMatrix(peak,windows = promoter)
tagHeatmap(tagMatrix)

# 注释可视化
plotAnnoBar(peakAnno)
vennpie(peakAnno)
plotAnnoPie(peakAnno)
plotDistToTSS(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)

# 查看peaks长度分布，只统计1000bp一下的peaks
peaksLength=abs(peakAnno_df$end-peakAnno_df$start)
peaksLength=peaksLength[peaksLength<1000]
hist(peaksLength, breaks = 5, col = "lightblue", xlim=c(0,10), xlab="peak length", main="histogram of peak length")

展示其中一个饼图，上述的可视化需要自行修改参数4.peaks相关基因注释

# 得到peaks相关基因之后再用去做富集分析即可
genes=unique(peakAnno_df$ENSEMBL)

此外还可以homer注释peaks，后面也会尝试学习一下~

既往推文

1. 转录组上游分析流程(一）：转录组上游分析流程(一）
2. 转录组上游分析流程(二）：转录组上游分析流程(二）
3. 转录组上游分析流程(三）：转录组上游分析流程(三）
4. 转录组上游分析流程(四）：转录组上游分析流程(四）
5. Chip-seq上游分析流程学习(一)：Chip-seq上游分析流程学习(一)
6. Chip-seq上游分析流程学习(二)：Chip-seq上游分析流程学习(二)
7. Chip-seq上游分析流程学习(三)：Chip-seq上游分析流程学习(三)

参考资料

1. 生信技能树：https://mp.weixin.qq.com/s/KOFd60USQdKY2C4Sc4OIvA  https://mp.weixin.qq.com/s/8zwJKzdNhnmqMdSjZjOYTQ https://mp.weixin.qq.com/s/8kOLNNrwNdfrUE_Vwcpaag https://mp.weixin.qq.com/s/i4gdRGcDAnPwhdzse_OsDw
2. 徐洲更hoptop：https://blog.csdn.net/u012110870/article/details/102804213 https://www.jianshu.com/p/a7b6ce208f98 https://www.jianshu.com/p/a7b6ce208f98
3. 生信星球：https://mp.weixin.qq.com/s/e4smCjDmLXyVCeWuvNWKWw
4. 生信学习日记：https://mp.weixin.qq.com/s/z-qbVZqo-cKT3Tfv4LT6OQ
5. Tobin笔记：https://mp.weixin.qq.com/s/Uf8GgK4kZAWbRKKoh6IiyA
6. 老俊俊的生信笔记：https://mp.weixin.qq.com/s/pfRiLhP8ONMNGcUI06IPog https://mp.weixin.qq.com/s/JLRf50wXLKcUVRt1kZbp1Q
7. 生信媛：https://mp.weixin.qq.com/s/iHnogzM6QSu8tlImE26_Wg
8. IGV可视化工具-搬运视频：https://www.bilibili.com/video/BV19y4y147uq/?spm_id_from=333.337.search-card.all.click&vd_source=3a13860df939bc922ad1fd6099e42c1d
9. Z_bioinfo: https://www.jianshu.com/p/d6cb795af22a
10. 小白鱼的生统笔记：https://mp.weixin.qq.com/s/VEczKns3Dlmva71j59Ne_g

致谢：感谢曾老师以及生信技能树团队全体成员。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟