m6A图文复现06-样本相关性检验与Peak Calling
m6A图文复现06-样本相关性检验与Peak Calling
前面我们分享了 跟着Nature Medicine学MeDIP-seq数据分析，数据和代码都是公开，这个2G的压缩包文件，足以学习3个月，写60篇教程。同时也分享了 全套MeRIP-seq文章图表复现代码，其实MeRIP-seq其实就是RNA水平的，又叫做m6a测序。
下面是MeRIP-seq 图表复现笔记
m6A图文复现已经更新到第五期了：（额因为排版问题没有第四期）
群里还是有非常多的小伙伴在问哪个步骤用哪个软件，在这里，我给出一张m6A分析的常规pipeline：
图出处：https://repicmod.uchicago.edu/repic/statistics.php
这个网站集合了非常多的目前已经发表的m6A测序数据，Peak Calling等方法结果以及比较。
1、样本相关性分析
我们前面已经完成了数据比对，在进行Peak Calling之前我们先来看看这几个样本之间的相关性，使用deepTools工具包来看看生物学重复是否聚集到了一起，分组信息为：
输入数据可以是bam也可以是bigwig，-bs参数默认为10000，可以适当调整这个区间，对相关性计算出来的结果影响还挺大，作者提供的代码用的10，耗时比较久出来结果也不是很好，在原有代码上我还添加了 --plotNumbers 参数在图片上显示相关性大小。
此外corMethod可以选择pearson或者spearman，-o 可以输出为heatmap_pearsonCor.pearson.pdf 的pdf格式或者heatmap_pearsonCor.pearson.png的png格式
# pearson correlation between samples
multiBamSummary bins -bs 1000 -p 10 --bamfiles *bam -o results.npz
plotCorrelation \
-in results.npz \
--corMethod pearson \
--skipZeros \
--plotNumbers \
--labels KO1_MeRIP KO1_NoIP KO2_MeRIP KO2_NoIP KO3_MeRIP KO3_NoIP WT1_MeRIP WT1_NoIP WT2_MeRIP WT2_NoIP WT3_MeRIP WT3_NoIP \
--plotTitle "m6A Correlation across samples" \
--whatToPlot heatmap --colorMap Blues \
-o heatmap_pearsonCor.pearson.pdf --removeOutliers \
--outFileCorMatrix pearsonCorr.rmOut.tab
结果如下，可以看到IP样本聚集在一起，Input样本聚集在一起，WT组和KO组也能看到两个色块。结果还行。
2、Peak Calling
从上面那张流程图里我们可以看到Peak Calling有三个软件：MACS2，exomePeak/exomePeak2，MeTPeak，也是目前市面上m6A分析数据中最常用的三款软件。
关于其中两个软件有篇文献做了测评和比较：
目前对MeRIP-Seq数据进行m6A peak calling分析的软件有两类，一类是早先为ChIP-Seq数据分析所研发的软件，如MACS；另一类则是专门为转录组m6A位点检测而开发的如exomepeak，HEPeak及在前两者基础上发展的MeTPeak。 利用已发表的两份MeRIP-Seq数据，对MACS和MeTPeak两种方法进行了对比，发现尽管MACS可以获得更多的Peak数，但使用MeTPeak所获得m6A Peaks更具有m6A主要发生的位点RRACH（R=G/A; H=A/C/U）这一motif的富集的特点，在MACS 单独发现的Peaks中有超过20%的Peaks分布在TSS这一m6A修饰较少发生的区域。此外，MeTPeak的分析具有链特异性，并且能更好的在剪接的外显子（junction exons）区域发现Peaks，因而在MeRIP-Seq数据分析中，MeTPeak更为适用。 ref：[1] Zeng, Y. , et al. "Refined RIP-seq protocol for epitranscriptome analysis with low input materials." PLoS Biology 16.9(2018):e2006092.
本教程也不太推荐使用MACS2，更加倾向于选择exomePeak/exomePeak2，MeTPeak这些专门为m6A测序数据开发的软件来进行分析。
首先第一个问题，应该也是大多数人会遇到的：Peak Calling的时候有两种做法，分为Sample和Group两种
- 单个样本IP与自身Input进行Peak Calling：这样每个样本就会得到一个Peak结果，对于生物学重复，后面可以合并Peak取交集部分
- 同组的多个IP放在一起，多个Input合并在一起进行Peak Calling：这样一个组就一个Peak结果。
第一种要求每个样本都有自身的Input样本；第二种可以不要求，可以是三个IP，两个Input之类的。
exomePeak2是exomePeak的升级版本，可以参考：孟佳课题组|Peak Calling软件exomePeak以及exomePeak2使用指北，本教程此次使用exomePeak2作为Peak Calling分析。
使用exomePeak2进行Peak Calling
由于exomePeak2为R包，输入数据为bam文件，数据比较大，耗时比较久，代码就写成传参脚本然后在服务器上提交后台运行。
R传参脚本：
rm(list=ls())
options(stringsAsFactors = F)
library(getopt)
spec <- matrix(c(
'help', 'h', 0, 'logical',
'gtf', 'g', 1, 'character',
'bamdir', 'b', 1, 'character',
'ip' , 'i', 1, 'character',
'input' , 'I', 1, 'character',
'paired', 'p', 2, 'character',
'lib', 'l', 2, 'character',
'mode', 'm', 2, 'character',
'pvalue', 'P', 2, 'numeric',
'fc', 'f', 2, 'numeric',
'od', 'o', 1, 'character'),
byrow = TRUE, ncol = 4)
opt <- getopt(spec)
print_usage <- function(spec=NULL){
cat(' Rscript exomePeak2.R --gtf Mus_musculus.gtf --bamdir Hisat2 --ip ip1,ip2 --input input1,input2 --od ./
Options:
--help -h NULL
--gtf character gtf file, genome annotation [forced]
--bamdir character mapping file, bam dir [forced]
--ip character IP Sample, used for peak calling only [optional]
--input character INPUT Sample, used for peak calling only [optional]
--paired character logical of whether the data comes from the Paired-End Library
default: TRUE
--lib character the protocal type of the RNA-seq library, can be
one in c("unstranded", "1st_strand", "2nd_strand");
default = "1st_strand"
--mode character a character specifies the scope of peak calling on genome, can be
one of c("exon", "full_transcript", "whole_genome");
Default = "whole_genome"
--pvalue numeric the cutoff on p values in peak calling, default = 1e-05
--fc numeric the cutoff on IP over input fold changes in peak calling, default = 2
--od character outdir [forced]
\n')
q(status=1)
}
if ( !is.null(opt$help) | is.null(opt$gtf) | is.null(opt$bamdir) ) { print_usage(spec) }
if ( is.null(opt$ip) |is.null(opt$input)) { print_usage(spec) }
if ( is.null(opt$paired) ) { opt$paired <- TRUE }
if ( is.null(opt$lib) ) { opt$lib <- "1st_strand" }
if ( is.null(opt$mode) ) { opt$mode <- "whole_genome" }
if ( is.null(opt$pvalue) ) { opt$pvalue <- 1e-05 }
if ( is.null(opt$fc) ) { opt$fc <- 2 }
if (!file.exists(opt$od)) { dir.create(opt$od) }
library(exomePeak2)
## peak calling for every sample or groups
ip <- unlist(strsplit(opt$ip,split=','))
input <- unlist(strsplit(opt$input,split=','))
ip_bam <- paste0(opt$bamdir,'/',ip,'/',ip,'.Hisat_aln.sorted.bam')
input_bam <- paste0(opt$bamdir,'/',input,'/',input,'.Hisat_aln.sorted.bam')
print(paste('ip_bam:',ip_bam))
print(paste('input_bam',input_bam))
txdb <- GenomicFeatures::makeTxDbFromGFF(file = opt$gtf, format="gtf", dataSource="Ensembl")
result <- exomePeak2(bam_ip = ip_bam,
bam_input = input_bam,
# data type
paired_end = opt$paired,
library_type = opt$lib,
txdb = txdb,
# cut off
p_cutoff = opt$pvalue,
log2FC_cutoff = log2(as.numeric(opt$fc)),
peak_calling_mode = opt$mode,
# save
save_dir = opt$od,
)
## result write out default
服务器提交后台运行：
# 手动造一个config文件，内容为下：
cat group.txt
#MeRIP No IP NAME
SRR866991 SRR866992 KO1
SRR866993 SRR866994 KO2
SRR866995 SRR866996 KO3
SRR866997 SRR866998 WT1
SRR866999 SRR867000 WT2
SRR867001 SRR867002 WT3
# 生成sh脚本
cat group.txt |awk 'NR>1{print"nohup Rscript exomePeak2.R --gtf Mus_musculus.GRCm38.104.gtf.gz --bamdir mapping/ --ip "$1" --input " $2 " --paired FALSE --lib unstranded --od Peak_Calling/"$3 " >Peak_Calling/$3.log &"}' >exomePeak2.sh
# sh脚本内容一行示例如下
nohup Rscript exomePeak2.R --gtf Mus_musculus.GRCm38.104.gtf.gz --bamdir mapping/ --ip SRR866991 --input SRR866992 --paired FALSE --lib unstranded --od Peak_Calling/KO1 >Peak_Calling/$3.log &
# 提交后台运行
sh exomePeak2.sh &
下一期分享运行结果以及结果可视化！