探索生信之美，解构每段代码后的逻辑

大家好呀，我是风间琉璃。在我们前面的几期甲基化分析流程讲解中相继讲解了minfi包、CHAMP包等分析方法（参照：“手把手教你甲基化生信分析——甲基化minfi包的使用（一）”、“手把手教你甲基化生信分析——甲基化minfi包的使用（二）”、“五分钟学会甲基化芯片处理，快上车！！！”、“解锁高阶甲基化分析流程？你的SCI又多了一张图”等等）。而在年，我发现新的一个进行甲基化分析的R包，又又又上新了！！如果只有前面几个包的功能，我就不会拿出来让大家重复看一遍，但是呢，这一个R包以更加简洁灵活的姿态无缝衔接到甲基化K芯片分析当中。整合多种方法对甲基化位点（DMP）、甲基化区域（DMR）以及特定区域甲基化差异多个分析方法进行分析。并提供多种可视化方式。在CNS级别的文章中，经常可以看到这样的甲基化分子实验图表。这个R包究竟有多好用，我们今天就一探究竟~废话不多说，我们就直接开始啦。

一、数据导入

#首先下载包#if(!rquirNamspac(BiocManagr,quitly=TRUE))#install.packags(BiocManagr)#BiocManagr::install(MEAL)#如果提示缺啥就下载啥

MEAL包分析的是已经预处理好的数据，比如minfi包处理之后的GnomicRatioSt格式的数据。我们就以minfiData包中的“MstEx”为例。

MstEx是MthylationRatioSt数据格式，包含了个CpG位点和6个样本，以及多个表型变量：如年龄、性别等。因为MEAL包是只适用于GnomicRatioSt，所以我们的数据“MstEx”还需要经过数据的转换。

我们来看看吧~

#加载我们分析中包library(MEAL)library(MultiDataSt)library(minfiData)library(minfi)library(ggplot2)#加载数据data("MstEx")#映射到基因组上，完成格式转换mth-mapToGnom(ratioConvrt(MstEx))#添加注释信息为我们的rowDatarowData(mth)-gtAnnotation(mth)接下来我们还需要对探针进行预处理，包括去除SNP以及包含NA等无用探针等。#去除SNP和CH开头的探针mth-dropMthylationLoci(mth)#去除带有SNP探针mth-dropLociWithSnps(mth)#去除NA，这里是只要有NA都去除mth-mth[!apply(gtBta(mth),1,function(x)any(is.na(x))),]#因为数据量太大啊，我们就选择10万个探针进行分析#设置随机数种子，不迷路st.sd(0)mth-mth[sampl(nrow(mth),1+05),]

二、数据分析

数据准备好了之后我们就可以开始差异分析了，MEAL包的流程特别简洁，我们只需要一个命令就能完成所有操作。我们直接开始吧~这里每个

rsAdj-runPiplin(st=mth,variabl_nams="status",covariabl_nams="ag",#纳入模型的协变量名称modl=NULL,#模型矩阵或者公式wights=NULL,#是否设置权重，通常不需要num_vars,#变量数，一般都不用关sva=FALSE,#是否进行btas=TRUE,#是选择bta值，否选择M值analyss=c("DiffMan","DiffVar"),#比较甲基化位点差异vrbos=TRUE,#展示运算过程warnings=TRUE,#展示warningDiffMan_params=NULL,#运行DiffMEAN的参数设置DiffVar_params=list(cofficint=1:2),#运行DiffVar参数设置rda_params=NULL,#RDA参数设置mthod="ls")#展示一下rsAdj##ObjctofclassRsultSt##.cratdwith:runPiplin##.sva:no##.#rsults:2(rror:0)##.faturData:probsx38variabls

这里需要给大家说一下，MEAL包具有两种不同的方法来分析差异甲基化位点（DMP），分别是比较平均值的limma包分析以及比较变异度的DiffVar方法，后面我们会以图形的方式呈现。而对于差异甲基化区域（DMR）的分析，MEAL包提供了三种算法，1.bumphuntr；2.blockFindr；3.DMRcat。还有就是针对特定区域的甲基化差异分析方法（RDA）。而我们分析出的结果最后得到的是RsultSt对象，并且不同的分析结果在不同的functionnams中。

#使用功能gtAssociation获得数据框格式的结果查看平均值具有差异的甲基化位点had(gtAssociation(objct=rsAdj,rid="DiffMan"))##logFCCI.LCI.RAvExprtP.Valu##cg-0.-0.-0..-42..-07##cg-0.-0.-0..-38..-06##cg-0.-0.-0..-35..-06##cg-0.-0.-0..-35..-06##cg-0.-0.-0..-34..-06##cg......-06##adj.P.ValBSE##cg0...##cg0...##cg0...##cg0...##cg0...##cg...#查看变异度不同的差异甲基化位点had(gtAssociation(objct=rsAdj,rid="DiffVar"))##logFCCI.LCI.RAvExprtP.Valu##cg-2.-3.-1..-6..##cg-1.-2.-1..-6..##cg-1.-2.-0..-6..##cg-2.-3.-1..-6..##cg-1.-2.-0..-6..##cg-1.-2.-1..-6..##adj.P.ValBSE##cg0...##cg0...##cg0...##cg0...##cg0...##cg0...

另外一个功能则是gtProbRsults，该功能能够帮助使用者分析CpG位点甲基化水平与不同变量线性回归系数，同样包含objct和rid两个参数，cof参数则是计算模型参数需要的cofficint数量，如果要同时分析不同数量的cofficint则可输入向量。fNams则是我们想从Faturdata选择的变量，最后整合到结果中。

had(gtProbRsults(objct=rsAdj,rid="DiffMan",cof=2:3,fNams=c("chromosom","start","nd")))#纳入染色体起止点##statusnormalagAvExprFP.Valu##cg-0.-0..929..928-06##cg-0.-0..730..-06##cg-0.-0..635..-06##cg-0.0..623..-06##cg-0.0..620..-06##cg.-0....-06##adj.P.ValSE.statusnormalSE.agchromosomstartnd##cg0...chr##cg0...chr##cg0...chr##cg0..390.chr##cg0..420.chr##cg...chr

同样我们可以使用gtGnVals功能，选择特定基因对应的CpG位点结果进行展示。

gtGnVals(rsAdj,"ARMS2",gncol="UCSC_RfGn_Nam",fNams=c("chromosom","start","nd"))##logFCCI.LCI.RAvExprtP.Valu##cg......##cg.-0.....##adj.P.ValBSEchromosomstartnd##cg.-4..chr##cg.-4..chr##UCSC_RfGn_Nam##cgARMS2##cgARMS2

三、结果的可视化

MEAL提供了多个可视化的结果，包括我们常用的曼哈顿图、火山图、基因组可视化图等等，我们一一来看看吧~

01

#简单全染色体可视化的曼哈顿图plot(rsAdj,rid="DiffMan",typ="manhattan")

#将靶点区域高亮展示targtRang-GRangs("23:-")plot(rsAdj,rid="DiffMan",typ="manhattan",highlight=targtRang)

#只绘制特定区间的曼哈顿图plot(rsAdj,rid="DiffMan",typ="manhattan",subst=targtRang)

#添加标准线以及其他内容plot(rsAdj,rid="DiffMan",typ="manhattan",suggstivlin=3,gnomwidlin=6,main="MycustomManhattan")ablin(h=13,col="yllow")

02

#设置阈值tPV以及tFC，分别是-log10ofPvalu以及tFCplot(rsAdj,rid="DiffMan",#diffrntMan的结果进行火山图呈现typ="volcano",#火山图tPV=14,tFC=0.4,show.labls=FALSE)+ggtitl("MycustomVolcano")

03

##展示CpG位点：cg在两组间的平均值差异plotFatur(st=mth,fat="cg",variabls="status")+ggtitl("DiffMans")

##展示CpG位点：cg在两组间的变异度差异plotFatur(st=mth,fat="cg",variabls="status")+ggtitl("DiffVars")

04

#包含需要绘制的染色体区间以及结果targtRang-GRangs("chrX:-")plotRgion(rsAdj,targtRang)

我们可以看上面面展示的是DiffMan和DiffVar两种方法得到的甲基化位点的logFC和-log10p-valu值。

#绘制特定的结果plotRgion(rsAdj,targtRang,rsults=c("DiffMan"),tPV=10)

好啦，本次的分析就到这里啦，今天我们介绍了MEAL如何简单快捷进行分析，并且提供多个可视化美图，大家赶紧用起来喔。我是风间琉璃，我们下期见~

往期传送门

手把手教你甲基化生信分析——甲基化minfi包的使用（一）

手把手教你甲基化生信分析—甲基化minfi包的使用（二）

手把手教你甲基化分析——甲基化CHAMP包的使用（一）

五分钟学会甲基化芯片处理，快上车！！！

一首歌的时间，我就掌握了这个生信技能

解锁高阶甲基化分析流程？你的SCI又多了一张图

年了！！这样的文章也能发近5分！

甲基化分析实战，将你的数据用在刀刃上！

小白都能学会的基因组可视化，不来看看吗？

万字长文教你入门biocondutor基石，小白都能看懂～

进阶！！！biocondutor基石的熟练应用

顶刊青睐的生信套路，小白都能学会，你不掌握下吗？

原来高分文章的图是这么做的？真简单！！！！

什么？UCSCgnbrowsr美图R也能画？新年送美图，牛年红红火火

IF14+非肿瘤的顶刊单细胞套路，亮点在这里！

年，加上转录因子还能发3+分！！

转录因子+LncRNA=6+分？

掌握这个技能，探针注释再也不是问题！！！

震惊！甲基化和转录因子可以这样结合？

这个R包，学会省你半年的实验时间！别再傻傻用实验验证甲基化和转录因子结合了！

太香了！10+顶刊里都爱的美图，一刻钟教你学会！瞬间变高级！（附代码）

不容错过！CNS级别的多组学Circos图是怎么做的？简单又实用，这个美图技能你必须gt！

实用又美观！频频登上高分生信期刊，这个图，让你脱颖而出！（附美图代码）

甲基化最常用的分析！没有之一！

多组学分析工具一键式搞定？高分文章的秘诀，学会能给你的SCI加分！

秀到飞起！生信中高逼格的多组学分析怎么做的？半小时轻松教你搞定！（附代码）

好东西啊！审稿人直呼内行！多组学，多平台分析终极解决方案！不学你就out了！

宝，在吗？送你个分子分型的绝招，包你在组会上被导师狠狠夸！

太猛了！年这个套路发了篇5+生信文章，还不快学！

欢迎大家