小鬼的m6A圖文復(fù)現(xiàn)06-樣本相關(guān)性檢驗(yàn)與Peak_Calling

m6A圖文復(fù)現(xiàn)已經(jīng)更新到第五期了
群里還是有非常多的小伙伴在問哪個(gè)步驟用哪個(gè)軟件，在這里，我給出一張m6A分析的常規(guī)pipeline：

image

圖出處：https://repicmod.uchicago.edu/repic/statistics.php

這個(gè)網(wǎng)站集合了非常多的目前已經(jīng)發(fā)表的m6A測(cè)序數(shù)據(jù)，Peak Calling等方法結(jié)果以及比較。

1、樣本相關(guān)性分析

我們前面已經(jīng)完成了數(shù)據(jù)比對(duì)，在進(jìn)行Peak Calling之前我們先來(lái)看看這幾個(gè)樣本之間的相關(guān)性，使用deepTools工具包來(lái)看看生物學(xué)重復(fù)是否聚集到了一起，分組信息為：

image

輸入數(shù)據(jù)可以是bam也可以是bigwig，-bs參數(shù)默認(rèn)為10000，可以適當(dāng)調(diào)整這個(gè)區(qū)間，對(duì)相關(guān)性計(jì)算出來(lái)的結(jié)果影響還挺大，作者提供的代碼用的10，耗時(shí)比較久出來(lái)結(jié)果也不是很好，在原有代碼上我還添加了 --plotNumbers 參數(shù)在圖片上顯示相關(guān)性大小。

此外corMethod可以選擇pearson或者spearman，-o 可以輸出為heatmap_pearsonCor.pearson.pdf 的pdf格式或者h(yuǎn)eatmap_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

結(jié)果如下，可以看到IP樣本聚集在一起，Input樣本聚集在一起，WT組和KO組也能看到兩個(gè)色塊。結(jié)果還行。

image

2、Peak Calling

從上面那張流程圖里我們可以看到Peak Calling有三個(gè)軟件：MACS2，exomePeak/exomePeak2，MeTPeak，也是目前市面上m6A分析數(shù)據(jù)中最常用的三款軟件。

關(guān)于其中兩個(gè)軟件有篇文獻(xiàn)做了測(cè)評(píng)和比較：

目前對(duì)MeRIP-Seq數(shù)據(jù)進(jìn)行m6A peak calling分析的軟件有兩類，一類是早先為ChIP-Seq數(shù)據(jù)分析所研發(fā)的軟件，如MACS；另一類則是專門為轉(zhuǎn)錄組m6A位點(diǎn)檢測(cè)而開發(fā)的如exomepeak，HEPeak及在前兩者基礎(chǔ)上發(fā)展的MeTPeak。

利用已發(fā)表的兩份MeRIP-Seq數(shù)據(jù)，對(duì)MACS和MeTPeak兩種方法進(jìn)行了對(duì)比，發(fā)現(xiàn)盡管MACS可以獲得更多的Peak數(shù)，但使用MeTPeak所獲得m6A Peaks更具有m6A主要發(fā)生的位點(diǎn)RRACH（R=G/A; H=A/C/U）這一motif的富集的特點(diǎn)，在MACS 單獨(dú)發(fā)現(xiàn)的Peaks中有超過20%的Peaks分布在TSS這一m6A修飾較少發(fā)生的區(qū)域。此外，MeTPeak的分析具有鏈特異性，并且能更好的在剪接的外顯子（junction exons）區(qū)域發(fā)現(xiàn)Peaks，因而在MeRIP-Seq數(shù)據(jù)分析中，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測(cè)序數(shù)據(jù)開發(fā)的軟件來(lái)進(jìn)行分析。

首先第一個(gè)問題，應(yīng)該也是大多數(shù)人會(huì)遇到的：Peak Calling的時(shí)候有兩種做法，分為Sample和Group兩種

  • 單個(gè)樣本IP與自身Input進(jìn)行Peak Calling：這樣每個(gè)樣本就會(huì)得到一個(gè)Peak結(jié)果，對(duì)于生物學(xué)重復(fù)，后面可以合并Peak取交集部分
  • 同組的多個(gè)IP放在一起，多個(gè)Input合并在一起進(jìn)行Peak Calling：這樣一個(gè)組就一個(gè)Peak結(jié)果。

第一種要求每個(gè)樣本都有自身的Input樣本；第二種可以不要求，可以是三個(gè)IP，兩個(gè)Input之類的。

exomePeak2是exomePeak的升級(jí)版本，可以參考：孟佳課題組|Peak Calling軟件exomePeak以及exomePeak2使用指北，本教程此次使用exomePeak2作為Peak Calling分析。

使用exomePeak2進(jìn)行Peak Calling

由于exomePeak2為R包，輸入數(shù)據(jù)為bam文件，數(shù)據(jù)比較大，耗時(shí)比較久，代碼就寫成傳參腳本然后在服務(wù)器上提交后臺(tái)運(yùn)行。

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

服務(wù)器提交后臺(tái)運(yùn)行：

# 手動(dòng)造一個(gè)config文件，內(nèi)容為下：
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腳本內(nèi)容一行示例如下
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 &

# 提交后臺(tái)運(yùn)行
sh exomePeak2.sh &

下一期分享運(yùn)行結(jié)果以及結(jié)果可視化！

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容