探索生信之美，解构每段代码背后的逻辑
大家好呀，我是风间琉璃，上一期我们讲解了m6A-seq(MeRIP-Seq)分析在五分文章中的应用（“这波方案火了！下一个国自然申请方向？用这些分析当你的前期工作，拿到直接用！”）。这一期我们就开始进行代码实操环节进一步讲解怎么分析m6A-seq(MeRIP-Seq)数据。
本期我们讲解的是MeRIP-Seq数据的分析利器——exomePeak2。exomePeak2对RNA甲基化免疫沉淀测序数据（MeRIP-Seq）提供了偏差量化和峰值检测的方法。MeRIP-Seq能够应用测序手段在给定细胞系的条件下对RNA修饰位点的丰度进行测量。
但是MeRIP-Seq存在几个弊端：
1. 和其他二代测序的技术一样，对PCR放大的偏倚格外敏感；
2. 在生物学重复样本间，测序的count数据具有显著的生物学改变。
而exomePeak2内置的算法则能够解决这一系列问题。通过exomePeak2包能够进行RNA甲基化位点peak鉴定、peak定量以及差异分析。这些步骤能够通过一站式，也可以通过多步骤进行分析。为了演示具体的分析流程，我们以多步骤解析先进行讲解~
在正式开始之前我们还需要大致讲解一下MeRIP-Seq，首先我们从样本中分离纯化出RNA，接下来将RNA裂解为50-100碱基长度的片段。接下来以分为二，其中一部分通过特异性抗体将甲基化的RNA序列片段进行沉淀后进行测序,这就是我们的IP数据；而另一部分则直接进行测序作为Input contol数据。后者的作业相当于进行背景校正，从而对甲基化的peak峰值进行定量。而我们经过上游的片段基因组比对后得到bam文件就可以开始我们的下游分析了。我们通过输入bam文件，得到我们peak的bed文件。我们就开始来看看吧~
讲完大致流程我们直接开始吧~
一、数据准备
我们进行MeRIP-seq分析时首先需要进行数据准备，这时需要将片段比对储存的bam文件，每一个样本都有bam_input和bam_ip两个文件，分别是背景测序文件和免疫沉淀测序文件。除此之外我们还需要转录本的注释文件，这里我们输入的是第7号染色体的GFF文件格式。同样我们也可以使用TxDb文件进行分析。
# 安装包 if(!requireNamespace('BiocManager', quietly = TRUE))# install.packages('BiocManager') BiocManager::install('exomePeak2') 加载包library(exomePeak2)# 基因注释文件目录GENE_ANNO_GTF = system.file("extdata", "example.gtf", package = "exomePeak2")# 获取IP样本的bam文件目录f1 = system.file("extdata", "IP1.bam", package = "exomePeak2")f2 = system.file("extdata", "IP2.bam", package = "exomePeak2")f3 = system.file("extdata", "IP3.bam", package = "exomePeak2")f4 = system.file("extdata", "IP4.bam", package = "exomePeak2")# 目录进行整合IP_BAM = c(f1, f2, f3, f4)# 获取input样本的bam文件目录f1 = system.file("extdata", "Input1.bam", package = "exomePeak2")f2 = system.file("extdata", "Input2.bam", package = "exomePeak2")f3 = system.file("extdata", "Input3.bam", package = "exomePeak2")# 目录进行整合INPUT_BAM = c(f1, f2, f3)# treat组的IP数据目录f1 = system.file("extdata", "treated_IP1.bam", package = "exomePeak2")TREATED_IP_BAM = c(f1)# treat组的Input数据目录f1 = system.file("extdata", "treated_Input1.bam", package = "exomePeak2")TREATED_INPUT_BAM = c(f1)
二、检查BAM文件
exomPeak2的scanMeripBAM()将BAM进行整合，并且将FLAG进行筛除。
MeRIP_Seq_Alignment <- scanMeripBAM(bam_ip = IP_BAM, bam_input = INPUT_BAM, bam_treated_input = TREATED_INPUT_BAM, bam_treated_ip = TREATED_IP_BAM, paired_end = FALSE #是否为双端测序)# 查看一下MeRIP_Seq_Alignment## MeripBamFileList of length 9## names(9): IP1.bam IP2.bam IP3.bam ... treated_IP1.bam treated_Input1.bam
三、外显子区域peak鉴定
通过exomePeakCalling()函数首先能够根据GFF注释文件确定外显子区间，并且对read进行计数以及广义线性模型（GLM）从而检测出具有显著改变的peak。除此之外，通过提供我们的基因组数据，比如hg19，从而校正不同样本的GC含量偏倚。
SummarizedExomePeaks <- exomePeakCalling(merip_bams = MeRIP_Seq_Alignment, gff_dir = GENE_ANNO_GTF, genome = "hg19")# 查看一下格式SummarizedExomePeaks## class: SummarizedExomePeak ## dim: 23 9 ## metadata(0):## assays(1): counts## rownames(23): peak_10 peak_11 ... control_5 control_6## rowData names(2): GC_content feature_length## colnames(9): IP1.bam IP2.bam ... treated_IP1.bam treated_Input1.bam## colData names(2): design_IP design_Treatment# 获得read count数据assays(SummarizedExomePeaks)## List of length 1## names(1): counts# 查看研究设计colData(SummarizedExomePeaks)## DataFrame with 9 rows and 2 columns## design_IP design_Treatment## <logical> <logical>## IP1.bam TRUE FALSE## IP2.bam TRUE FALSE## IP3.bam TRUE FALSE## IP4.bam TRUE FALSE## Input1.bam FALSE FALSE## Input2.bam FALSE FALSE## Input3.bam FALSE FALSE## treated_IP1.bam TRUE TRUE## treated_Input1.bam FALSE TRUE# GC含量以及feature lengthelementMetadata(SummarizedExomePeaks)## DataFrame with 23 rows and 2 columns## GC_content feature_length## <numeric> <integer>## 1 0.485000 600## 2 0.505000 200## 3 0.592000 250## 4 0.698462 325## 5 0.574815 675## ... ... ...## 19 0.687898 157## 20 0.593750 96## 21 0.608939 179## 22 0.680000 125## 23 0.536000 125# 修饰位点（Peak）的基因组起止rowRanges(SummarizedExomePeaks)## GRangesList object of length 23:## $peak_10## GRanges object with 1 range and 1 metadata column:## seqnames ranges strand | gene_id## <Rle> <IRanges> <Rle> | <character>## 10 chr7 150910205-150910654 - | ABCF2## -------## seqinfo: 1 sequence from an unspecified genome; no seqlengths## ## $peak_11## GRanges object with 1 range and 1 metadata column:## seqnames ranges strand | gene_id## <Rle> <IRanges> <Rle> | <character>## 11 chr7 150909780-150909829 - | ABCF2## -------## seqinfo: 1 sequence from an unspecified genome; no seqlengths## ## $peak_12## GRanges object with 2 ranges and 1 metadata column:## seqnames ranges strand | gene_id## <Rle> <IRanges> <Rle> | <character>## 12 chr7 150817227-150817232 + | AGAP3## 12 chr7 150819812-150819905 + | AGAP3## -------## seqinfo: 1 sequence from an unspecified genome; no seqlengths## ## ...## <20 more elements># 或者获得单碱基RNA修饰注释文件f2 = system.file("extdata", "mod_annot.rds", package = "exomePeak2")MOD_ANNO_GRANGE <- readRDS(f2)SummarizedExomePeaks <- exomePeakCalling(merip_bams = MeRIP_Seq_Alignment, gff_dir = GENE_ANNO_GTF, genome = "hg19", mod_annot = MOD_ANNO_GRANGE)
四、预测校正因素
在进行peak定量或者差异分析之前，我们还需要进行标准化从而校正不同样本因为测序深度而导致的差异。这里有estimateSeqDepth()和normalizeGC()两种功能。
# 绘制GC和log2FC之间的散点图plotLfcGC(SummarizedExomePeaks)
# 进行测序深度预测以及GC校正SummarizedExomePeaks<- estimateSeqDepth(SummarizedExomePeaks)SummarizedExomePeaks<- normalizeGC(SummarizedExomePeaks)# 再次绘制绘制GC和log2FC之间的散点图plotLfcGC(SummarizedExomePeaks)
# 有/无统计学意义的差异修饰位点用不同颜色进行标记# 接下来展示Background、ALL、Modification三种类型的测序深度plotSizeFactors(SummarizedExomePeaks)
五、输出peak的计算结果
最后分析的结果能够以CSV和bed文件格式进行输出。输出到exomePeak2_outpu文件夹中。
exportResults(SummarizedExomePeaks, format = "BED")
六、一站式分析
前面我们讲解了如果多个步骤进行分析，其实我们可以非常简单一步完成。我们可以直接通过一行代码进行差异分析，exomePeak2能够对外显子区域的peak鉴定定量后通过GLM方法进行差异分析
exomePeak2(bam_ip = IP_BAM, bam_input = INPUT_BAM, bam_treated_input = TREATED_INPUT_BAM, bam_treated_ip = TREATED_IP_BAM, gff_dir = GENE_ANNO_GTF, genome = "hg19", paired_end = FALSE, save_dir = "exomePeak_output")## class: SummarizedExomePeak ## dim: 23 9 ## metadata(0):## assays(2): counts GCsizeFactors## rownames(23): peak_10 peak_11 ... control_5 control_6## rowData names(2): GC_content feature_length## colnames(9): IP1.bam IP2.bam ... treated_IP1.bam treated_Input1.bam## colData names(3): design_IP design_Treatment sizeFactor# 保存在exomePeak_output文件夹中
我们可以看到结果输出在DiffMod文件中，有csv、rds和bed文件格式。我们可以进一步打开查看。包括每一个peak的基因组位置、不同分组的read count数、一步差异倍数和P值等等。大家可以打开详细看一下。
好啦，本次的分析到这儿就结束啦，这一期我们讲解了如何通过输入BAM文件，进行peak鉴定、peak定量以及peak差异分析，从而得到我们想要的结果。大家可以自己上手跑一下。
本期就到此结束啦，我是风间琉璃，我们下期见~
往期传送门
手把手教你甲基化生信分析——甲基化minfi包的使用（一）
手把手教你甲基化生信分析—甲基化minfi包的使用（二）
手把手教你甲基化分析——甲基化CHAMP包的使用（一）
五分钟学会甲基化芯片处理，快上车！！！
一首歌的时间，我就掌握了这个生信技能
解锁高阶甲基化分析流程？你的SCI又多了一张图
2020年了！！这样的文章也能发近5分！
甲基化分析实战，将你的数据用在刀刃上！
小白都能学会的基因组可视化，不来看看吗？
万字长文教你入门biocondutor基石，小白都能看懂～
进阶！！！biocondutor基石的熟练应用
顶刊青睐的生信套路，小白都能学会，你不掌握下吗？
原来高分文章的图是这么做的？真简单！！！！
什么？UCSC gene browser 美图R也能画？新年送美图，牛年红红火火
IF 14+ 非肿瘤的顶刊单细胞套路，亮点在这里！
2021年，加上转录因子还能发3+分！！
转录因子+LncRNA=6+分？
掌握这个技能，探针注释再也不是问题！！！
震惊！甲基化和转录因子可以这样结合？
这个R包，学会省你半年的实验时间！别再傻傻用实验验证甲基化和转录因子结合了！
太香了！10+顶刊里都爱的美图，一刻钟教你学会！瞬间变高级！（附代码）
不容错过！CNS级别的多组学Circos图是怎么做的？简单又实用，这个美图技能你必须get！
实用又美观！频频登上高分生信期刊，这个图，让你脱颖而出！（附美图代码）
甲基化最常用的分析！没有之一！
多组学分析工具一键式搞定？高分文章的秘诀，学会能给你的SCI加分！
秀到飞起！生信中高逼格的多组学分析怎么做的？半小时轻松教你搞定！（附代码）
好东西啊！审稿人直呼内行！多组学，多平台分析终极解决方案！不学你就out了！
宝，在吗？520送你个分子分型的绝招 ，包你在组会上被导师狠狠夸！
太猛了！2021年这个套路发了134篇5+生信文章，还不快学！
逆天！生信也能做出CNS湿实验的效果！真的惊艳到我了！
这波方案火了！下一个国自然申请方向？用这些分析当你的前期工作，拿到直接用！
—
END
—

撰文丨风间琉璃
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~

继续阅读
阅读原文
关键词
数据
基因组
结果
文章
步骤