探索生信之美，解构每一段代码的故事
大家好呀，我是风间琉璃，上一期我们详细介绍了ELMER的运行原理以及其创建MAE对象、筛选distalprob的流程（详见“震惊！甲基化和转录因子可以这样结合？”）。今天我们继续向下探索。完成对ELMER包的下游分析。
三、确定prob-gn对
这一步我们将会把远端探针（distalprob）的甲基化水平和靶基因的表达水平进行相关性分析，从而构建出prob-gn对。这个包具体的做法是分别找到差异甲基化远端探针（distalprob）的上游和下游最近的10个基因，分析其余对应探针甲基化水平是否存在负相关，从而筛选出来的潜在prob-gn对。
这一步ELMER包同样提供了两类方式：
（1）suprvisd，即通过上游分析疾病组vs.正常组差异远端探针，进一步纳入所有样本通过非参数检验的方式，找出同时在疾病组高表达的基因；
（2）unsuprvisd，根据基因的distalprob甲基化水平进行排序，提取出甲基化水平分别在前20%和后20%的样本，分别为M和U组，通过U检验方式比较两组间基因表达水平,如果在M组的表达水平显著低于U组，并通过迭代的方式计算出所有远端探针（distalprob）和靶基因的相关性Pvalu,从而筛选出具有意义的prob-gn关系对。
#加载我们之前的数据ma-gt(load("ma.rda"))#加载我们筛选得到的差异甲基化远端探针（distalprob）sig.diff-rad.csv("rsult/gtMthdiff.hypo.probs.significant.csv")#找到差异甲基化远端探针（distalprob）上下游10个基因narGns-GtNarGns(data=ma,probs=sig.diff$prob,numFlankingGns=20)#分别是上游10个和下游10基因##Sarchingforth20nargns##Idntifyinggnpositionforachprob#因为前面使用unsuprvisd，我们这而也用unsuprvisdHypo.pair-gt.pair(data=ma,group.col="dfinition",group1="PrimarysolidTumor",group2="SolidTissuNormal",narGns=narGns,mod="unsuprvisd",prmu.dir="rsult/prmu",prmu.siz=,#一般要设置迭代次数为00，处于演示目的减少次数raw.pvalu=0.05,P=0.01,filtr.probs=TRUE,filtr.prcntag=0.05,filtr.portion=0.3,dir.out="rsult",#结果同样输出到rsult目录cors=1,labl="hypo")##SlctingU(unmthylatd)andM(mthylatd)groups.Eachgroupshas20%ofsampls##-------------------##*Filtringprobs##-------------------##FormorinformationsfunctionprAssociationProbFiltring##Makingsurwhavatlast5%ofbtavaluslssrthan0.3and5%ofbtavalusgratr0.3.##Rmovingprobsoutof##CalculatingPp(prob-gn)forallnarbygnsCompltdaftr16s##Filcratd:rsult/gtPair.hypo.all.pairs.statistic.csv##CalculatingPr(randomprob-gn).Prmutatingprobsforgns##CalculatmpiricalPvalu.#展示一下prob-gn对had(Hypo.pair)##ProbGnIDSymbolDistancSids##cg.ENSGcgENSGLMO40L1##cg.ENSGcgENSGSYT1486R1##cg.ENSGcgENSGMLK4-L4##cg.ENSGcgENSGADCYR3##cg.ENSGcgENSGLMO40L1##cg.ENSGcgENSGGPRR6##Raw.pP##cg.ENSG7.-.##cg.ENSG1.-.##cg.ENSG2.-.##cg.ENSG4.-.##cg.ENSG4.-.##cg.ENSG8.-.
接下来就到了最激动人心的时刻，我们将开始我们最核心的一步——构建TF-gn调控网络。
四、筛选后探针的motif富集分析
这一步我们通过上一步分析得到的prob-gn对中的探针用以富集分析。提取探针上下游bp区域的碱基序列，从而找到富集的motif。接下来再通过转录因子（TF）结合motif数据库，从而预测出motif结合的转录因子。
nrichd.motif-gt.nrichd.motif(data=ma,probs=Hypo.pair$Prob,dir.out="rsult",labl="hypo",min.incidnc=10,lowr.OR=1.1)##Loadingobjct:Probs.motif.hg19.KCompltdaftr27s##RtrivingTFClassfamilyclassificationfromELMER.data.##RtrivingTFClasssubfamilyclassificationfromELMER.data.##------------------------------------##**Filtringmotifsbasdonquality##------------------------------------##Numbrofnrichdmotifswithquality:##-----------##=A:4##=B:1##=C:1##=D:0##=S:0##-----------##ConsidringonlymotifswithqualityfromAuptoDS:6motifsarnrichd.##---------------------------------------------------##*Addingnrichdmotifstosignificantpairsfil##---------------------------------------------------##Addingcoordinatsforprobsandgnsfromthprovidddata##Joining,by="GnID"##Joining,by="Prob"##Savingfil:rsult/gtPair.hypo.pairs.significant.withmotif.csv这一步里，我们需要输入我们的MAE对象，筛选出prob-gn对中的probs，并且还需要设置motif最低95%置信区间OR阈值以及最少几个prob进行富集motif。#展示一下名称nams(nrichd.motif)##[1]"HXA9_HUMAN.H11MO.0.B""FOXH1_HUMAN.H11MO.0.A""PO5F1_HUMAN.H11MO.1.A"##[4]"HXB13_HUMAN.H11MO.0.A""IRF9_HUMAN.H11MO.0.C""SOX10_HUMAN.H11MO.1.A"#展示一下每个motif是通过哪些探针进行富集的nrichd.motif[[1]]##[1]"cg""cg""cg""cg""cg94320"##[6]"cg""cg""cg""cg""cg"##[11]"cg""cg""cg""cg"
除此之外，ELMER包还提供了两个表：“gtMotif.hypo.nrichd.motifs.rda”“gtMotif.hypo.motif.nrichmnt.csv”，大家可以在rsult目录看一下。
五、确定调控作用转录因子
这一步，ELMER包对motif以及上游转录因子的关系对进行筛选，从而得到具有调控作用的转录因子。如果某一特定亚组中基因的nhancr发生改变，其上有调控的转录因子同样也会发生改变。基于此，ELMER包将上一步我们通过富集分析找到富集的motif以及对应motif的转录因子进行分析，筛选出对应转录水平发生改变的转录因子。
同样地，ELMER提供了两种方式：
（1）unsuprvisd，即将存在相同motif的远端探针（distalprob）根据甲基化水平分为高甲基化组（一般取前20%）和低甲基化组（后20%），比较其两组间对应TF的表达值是否存在差异。
（2）Suprvisd，即对疾病组vs.对照组相同motif对应TF表达值进行分析。接下来筛选出前5%Pvalu最小的TFs，并视为潜在的上游调控转录因子。
##找到调控motif对应的转录因子TF-gt.TFs(data=ma,group.col="dfinition",group1="PrimarysolidTumor",group2="SolidTissuNormal",mod="unsuprvisd",nrichd.motif=nrichd.motif,dir.out="rsult",cors=1,labl="hypo")##SlctingU(unmthylatd)andM(mthylatd)groups.Eachgroupshas20%ofsampls##-------------------------------------------------------------------------------------------------------------------##**DownloadingTFlistfromLambrt,SamulA.,tal.Thhumantranscriptionfactors.Cll.4():-.##-------------------------------------------------------------------------------------------------------------------##AccssingTFfamilisfromTFClassdatabastoindntifyknownpotntialTF##RtrivingTFClassfamilyclassificationfromELMER.data.##RtrivingTFClasssubfamilyclassificationfromELMER.data.##Calculatingthavragmthylationatallmotif-adjacntprobs#展示一下结果had(TF)##motiftop.potntial.TF.family##HXA9_HUMAN.H11MO.0.BHXA9_HUMAN.H11MO.0.BPOU6F2##FOXH1_HUMAN.H11MO.0.AFOXH1_HUMAN.H11MO.0.AFOXE1##PO5F1_HUMAN.H11MO.1.APO5F1_HUMAN.H11MO.1.APOU6F2##HXB13_HUMAN.H11MO.0.AHXB13_HUMAN.H11MO.0.AHOXD11##IRF9_HUMAN.H11MO.0.CIRF9_HUMAN.H11MO.0.CIRF6##SOX10_HUMAN.H11MO.1.ASOX10_HUMAN.H11MO.1.ASOX2##top.potntial.TF.subfamily##HXA9_HUMAN.H11MO.0.BHOXD12##FOXH1_HUMAN.H11MO.0.AFOXH1##PO5F1_HUMAN.H11MO.1.ANA##HXB13_HUMAN.H11MO.0.AHOXD11##IRF9_HUMAN.H11MO.0.CNA##SOX10_HUMAN.H11MO.1.ANA##potntial.TF.family##HXA9_HUMAN.H11MO.0.BPOU6F2;HOXD12##FOXH1_HUMAN.H11MO.0.AFOXE1;FOXN1;FOXH1;FOXD1##PO5F1_HUMAN.H11MO.1.APOU6F2;HOXA10##HXB13_HUMAN.H11MO.0.AHOXD11;HOXD12;HOXD13;HOXD10;HOXA10;POU6F2;HOXA11##IRF9_HUMAN.H11MO.0.CIRF6##SOX10_HUMAN.H11MO.1.ASOX2;SOX21;SOX15;SOX6##potntial.TF.subfamily##HXA9_HUMAN.H11MO.0.BHOXD12##FOXH1_HUMAN.H11MO.0.AFOXH1##PO5F1_HUMAN.H11MO.1.ANA##HXB13_HUMAN.H11MO.0.AHOXD11;HOXD12;HOXD13;HOXD10;HOXA10;HOXA11##IRF9_HUMAN.H11MO.0.CNA##SOX10_HUMAN.H11MO.1.ANA##top_5prcnt_TFs##HXA9_HUMAN.H11MO.0.BZFP64;ZNF74;ZNF;PATZ1;SOX2;OTX1;BCL11A;AHCTF1;ZNF;TFAP4;ZNF;PAX6;ZBTB33;SIX4;ZIC2;TFCP2;ZNF77;E2F6;TGIF2;ZNF;ZBTB12;TCF20;POU6F2;ZNF;MYNN;DMRT2;SOX12;DLX2;TP63;TBPL1;ZNF18;SP3;ZNF;SAFB;ZNF3;DLX6;ZNF;ZBTB39;RFX7;CREB3L4;ZNF;ZNF;DLX5;VAX2;ZNF;ZNF;HOXD12;SIX1;NEUROG2;ZNF26;SIM2;NEUROD2;ZNF;MAZ;DMRT1;FOXM1;DMRT3;DLX1;AEBP2##FOXH1_HUMAN.H11MO.0.ASOX2;OTX1;ZFP64;TP63;KLF5;NFE2L2;ZIC2;ZNF;FOXE1;SOX21;TFAP4;BCL11A;ZNF77;TBPL1;DMRT2;ZNF;PLAG1;DMRT3;DLX5;ZNF18;TCF20;ZNF;POU6F2;SOX6;ZNF;FOXN1;AHCTF1;OTX2;MYCL;SP3;ZBTB33;DLX6;ZNF3;FOXH1;ZNF;AEBP2;ZNF;PAX6;ZNF74;ZNF;GTF2IRD1;ZNF;FOXD1;MYNN;ZNF;PITX1;HOXD13;GRHL3;ZNF26;ZNF;SOX15;NRF1;DMTF1;RFX7;PRDM13;ZNF;SIX4;HOXD12;DMRT1##PO5F1_HUMAN.H11MO.1.ATFAP4;ZFP64;ZNF;OTX1;BCL11A;ZBTB33;PATZ1;TP63;DMRT2;SOX2;TCF20;ZNF;HES6;ZNF74;FOXN1;ZNF77;DLX5;POU6F2;KLF5;GTF2IRD1;PAX6;TFAP2B;NKRF;ZIC2;ZNF;SOX12;ZNF;TFCP2;AHCTF1;OTX2;HMGA1;TFAP2C;SOX15;CREB3L4;ZNF;MYNN;HOXA10;FOXH1;NFE2L2;ZBTB39;MYCL;ZNF;ZNF18;DMRT3;FOXE1;DLX6;ZNF26;ZNF;ZNF;MAZ;NRF1;PITX1;GLI4;ZBTB12;RFX7;SAFB;SOX21;SIX1;AEBP2##HXB13_HUMAN.H11MO.0.ABCL11A;ZFP64;ZNF74;OTX1;TP63;PATZ1;HOXD11;ZIC2;ZNF77;ZNF;ZNF18;PAX6;SOX2;AEBP2;TCF20;SIM2;TFAP4;HOXD12;TGIF2;KLF5;ZNF;HOXD13;HOXD10;ZNF26;ZNF;DLX5;HOXA10;POU6F2;ZBTB33;HOXA11;MYNN;FOXM1;DLX6;DMRT2;TBPL1;RFX7;AHCTF1;ZNF;ZNF;SOX21;GTF2IRD1;SIX4;ZNF;ZNF;TFCP2;ZBTB39;FOXE1;ZNF;E2F6;DMRT1;ZNF;MAZ;SOX6;TLX1;NRF1;EMX1;NFE2L2;ZNF;ZNF##IRF9_HUMAN.H11MO.0.CZFP64;OTX1;TFAP4;ZNF;ZNF74;TCF20;BCL11A;PATZ1;TP63;ZIC2;SOX2;SAFB;ZNF77;DMRT2;POU6F2;ZBTB33;AHCTF1;ZNF18;TFCP2;PAX6;SIM2;SIX1;ZNF;SREBF1;IRF6;SIX4;ZNF;ZNF;HOXD12;TBPL1;ZNF;TCF3;ZNF26;ZNF3;ZNF7;DMRT3;FOXH1;NKRF;HOXA11;HES6;AEBP2;ZNF;ZNF;DMRT1;ZBTB12;ZNF;PLAG1;ZNF;PAX9;TFAP2B;KLF5;FOXE1;HMGA1;FOXN1;SOX21;ZNF;ELK4;DLX5;SREBF2##SOX10_HUMAN.H11MO.1.ATP63;SOX2;KLF5;OTX1;DLX5;BCL11A;DMRT2;ZNF77;TFAP4;NFE2L2;ZNF;TCF20;MYNN;ZFP64;FOXE1;ZNF74;ZIC2;ZBTB33;ZNF;ZNF;SOX21;DLX6;POU6F2;ZNF18;DMRT3;DMRT1;AEBP2;ELK4;PITX1;ZNF;FOXN1;GRHL3;ZNF;RFX7;GTF2IRD1;ZNF;E2F6;GRHL1;ZNF;TBPL1;ZNFA;HOXD12;SP3;EMX1;SOX15;AHCTF1;PRDM13;ZNF3;ZBTB39;ZNF;HOXA10;ZNF;SOX6;PATZ1;PAX6;HMGA1;MYCL;MAZ;SMAD3
好啦，我们的分析就结束了，总结一下，我们通过同一批样本的甲基化数据和转录组数据，首先找到具有差异的远端探针（distalprob），并根据远端探针和对应基因的表达水平构建出prob-gn对。并将prob-gn对中的prob进一步进行motif富集分析，找出潜在富集motif，最后根据富集的motif推测出上游的调控TFs。这一条逻辑链比较复杂，大家可以再翻出上一期再品鉴一下。下一期我们讲介绍ELMER包的可视化以及高分文章上使用ELMER包的示例。
好啦，我是风间琉璃，咱们下期见~
往期传送门
手把手教你甲基化生信分析——甲基化minfi包的使用（一）
手把手教你甲基化生信分析—甲基化minfi包的使用（二）
手把手教你甲基化分析——甲基化CHAMP包的使用（一）
五分钟学会甲基化芯片处理，快上车！！！
一首歌的时间，我就掌握了这个生信技能
解锁高阶甲基化分析流程？你的SCI又多了一张图
年了！！这样的文章也能发近5分！
甲基化分析实战，将你的数据用在刀刃上！
小白都能学会的基因组可视化，不来看看吗？
万字长文教你入门biocondutor基石，小白都能看懂～
进阶！！！biocondutor基石的熟练应用
顶刊青睐的生信套路，小白都能学会，你不掌握下吗？
原来高分文章的图是这么做的？真简单！！！！
什么？UCSCgnbrowsr美图R也能画？新年送美图，牛年红红火火
IF14+非肿瘤的顶刊单细胞套路，亮点在这里！
年，加上转录因子还能发3+分！！
转录因子+LncRNA=6+分？
掌握这个技能，探针注释再也不是问题！！！
震惊！甲基化和转录因子可以这样结合？
—END—撰文丨风间琉璃排版丨四金兄值班风风主编丨小雪球
欢迎大家