m6ASeq RNA甲基化测序

m6A图文复现06-样本相关性检验与Peak_Calling

2022-01-06  本文已影响0人  信你个鬼

m6A图文复现已经更新到第五期了
群里还是有非常多的小伙伴在问哪个步骤用哪个软件,在这里,我给出一张m6A分析的常规pipeline:

image

图出处:https://repicmod.uchicago.edu/repic/statistics.php

这个网站集合了非常多的目前已经发表的m6A测序数据,Peak Calling等方法结果以及比较。

1、样本相关性分析

我们前面已经完成了数据比对,在进行Peak Calling之前我们先来看看这几个样本之间的相关性,使用deepTools工具包来看看生物学重复是否聚集到了一起,分组信息为:

image

输入数据可以是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组也能看到两个色块。结果还行。

image

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.

image

本教程也不太推荐使用MACS2,更加倾向于选择exomePeak/exomePeak2,MeTPeak这些专门为m6A测序数据开发的软件来进行分析。

首先第一个问题,应该也是大多数人会遇到的:Peak Calling的时候有两种做法,分为Sample和Group两种

第一种要求每个样本都有自身的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, 'logical',
  'lib',    'l',  2, 'character',
  'mode',   'm',  2, 'character',
  'cores',  'c',  2, 'numeric',
  'pvalue', 'P',  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
    --input  character    INPUT Sample, used for peak calling only

    --paired logical      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_tx", "whole_genome"); 
                          Default = "exon"

    --cores  numeric      a numeric value specifying the number of cores used for parallel computing; 
                          default = 6.

    --pvalue numeric      the cutoff on p values in peak calling
    --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 <- "exon" }
if ( is.null(opt$cores) )   { opt$cores <-  6 }
if ( is.null(opt$pvalue) )   { opt$pvalue <-  1e-5 }
if (!file.exists(opt$od))   { dir.create(opt$od) }

library(exomePeak2)
package.version("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")

if( !is.null(opt$pvalue) ) {
  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,
                     peak_calling_mode = opt$mode, 

                     # parallel running
                     parallel = opt$cores, save_dir = opt$od
                     )
}


## result write out default

服务器提交后台运行:

# 手动造一个config文件,内容为下:
cat group.txt
#MeRIP  NoIP    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 --pvalue 1e-05 --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 --pvalue 1e-05 --od Peak_Calling/KO1 >Peak_Calling/$3.log &

# 提交后台运行
sh exomePeak2.sh &

下一期分享运行结果以及结果可视化!

上一篇下一篇

猜你喜欢

热点阅读