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

　　探索生信之美，解构每一段代码的故事

　　更多生信干货教程

　　私信我回复“数据库”即可免费领取哦

　　大家好呀，我是风间琉璃，上一期我们详细介绍了ELMER的运行原理以及其创建MAE对象、筛选distal probe的流程（详见“震惊！甲基化和转录因子可以这样结合？”）。今天我们继续向下探索。完成对ELMER包的下游分析。

　　三、确定probe-gene对

　　这一步我们将会把远端探针（distal probe）的甲基化水平和靶基因的表达水平进行相关性分析，从而构建出probe-gene对。这个包具体的做法是分别找到差异甲基化远端探针（distal probe）的上游和下游最近的10个基因，分析其余对应探针甲基化水平是否存在负相关，从而筛选出来的潜在probe-gene对。

　　这一步ELMER包同样提供了两类方式：

　　（1）supervised，即通过上游分析疾病组vs.正常组差异远端探针，进一步纳入所有样本通过非参数检验的方式，找出同时在疾病组高表达的基因；

　　（2）unsupervised，根据基因的distal probe甲基化水平进行排序，提取出甲基化水平分别在前20%和后20%的样本，分别为M和U组，通过U检验方式比较两组间基因表达水平,如果在M组的表达水平显著低于U组，并通过迭代的方式计算出所有远端探针（distal probe）和靶基因的相关性P value,从而筛选出具有意义的probe-gene关系对。

　　# 加载我们之前的数据

　　mae

　　#加载我们筛选得到的差异甲基化远端探针（distal probe）

　　sig.diff

　　#找到差异甲基化远端探针（distal probe）上下游10个基因

　　nearGenes

　　probes=sig.diff$probe,

　　numFlankingGenes=20) # 分别是上游10个和下游10基因

　　## Searching for the 20 near genes

　　## Identifying gene position for each probe

　　#因为前面使用unsupervised，我们这而也用unsupervised

　　Hypo.pair

　　group.col="definition",

　　group1="Primary solid Tumor",

　　group2="Solid Tissue Normal",

　　nearGenes=nearGenes,

　　mode="unsupervised",

　　permu.dir="result/permu",

　　permu.size=100, # 一般要设置迭代次数为10000，处于演示目的减少次数

　　raw.pvalue=0.05,

　　Pe=0.01,

　　filterbes=TRUE,

　　filter.percentage=0.05,

　　filter.portion=0.3,

　　dir.out="result",#结果同样输出到result目录

　　cores=1,

　　label="hypo")

　　## Selecting U (unmethylated) and M (methylated) groups. Each groups has 20% of samples

　　## -------------------

　　## * Filtering probes

　　## -------------------

　　## For more information see function preAssociationProbeFiltering

　　## Making sure we have at least 5% of beta values lesser than 0.3 and 5% of beta values greater 0.3.

　　## Removing 592 probes out of 1642

　　## Calculating Pp (probe - gene) for all nearby genes

　　Completed after 16 s

　　## File created: result/getPair.hypo.all.pairs.statistic.csv

　　## Calculating Pr (random probe - gene). Permutating 100 probes for 747 genes

　　## Calculate empirical P value.

　　#展示一下probe-gene对

　　head(Hypo.pair)

　　## Probe GeneID Symbol Distance Sides

　　## cg20701183.ENSG00000143013 cg20701183 ENSG00000143013 LMO4 0 L1

　　## cg19403323.ENSG00000143469 cg19403323 ENSG00000143469 SYT14 87476 R1

　　## cg12213388.ENSG00000143674 cg12213388 ENSG00000143674 MLK4 -984182 L4

　　## cg26607897.ENSG00000143199 cg26607897 ENSG00000143199 ADCY10 267784 R3

　　## cg10574861.ENSG00000143013 cg10574861 ENSG00000143013 LMO4 0 L1

　　## cg26607897.ENSG00000143147 cg26607897 ENSG00000143147 GPR161 543156 R6

　　## Raw.p Pe

　　## cg20701183.ENSG00000143013 7.453984e-14 0.00990099

　　## cg19403323.ENSG00000143469 1.671937e-12 0.00990099

　　## cg12213388.ENSG00000143674 2.527644e-12 0.00990099

　　## cg26607897.ENSG00000143199 4.593610e-12 0.00990099

　　## cg10574861.ENSG00000143013 4.770162e-12 0.00990099

　　## cg26607897.ENSG00000143147 8.048248e-12 0.00990099

　　接下来就到了最激动人心的时刻，我们将开始我们最核心的一步——构建TF-gene调控网络。

　　四、筛选后探针的motif富集分析

　　这一步我们通过上一步分析得到的probe-gene对中的探针用以富集分析。提取探针上下游250bp区域的碱基序列，从而找到富集的motif。接下来再通过转录因子（TF）结合motif数据库，从而预测出motif结合的转录因子。

　　enriched.motif

　　probes=Hypo.pair$Probe,

　　dir.out="result",

　　label="hypo",

　　min.incidence=10,

　　lower.OR=1.1)

　　## Loading object: Probes.motif.hg19.450K

　　Completed after 27 s

　　## Retrieving TFClass family classification from ELMER.data.

　　## Retrieving TFClass subfamily classification from ELMER.data.

　　## ------------------------------------

　　## ** Filtering motifs based on quality

　　## ------------------------------------

　　## Number of enriched motifs with quality:

　　## -----------

　　##=> A: 4

　　##=> B: 1

　　##=> C: 1

　　##=> D: 0

　　##=> S: 0

　　## -----------

　　## Considering only motifs with quality from A up to DS: 6 motifs are enriched.

　　## ---------------------------------------------------

　　## * Adding enriched motifs to significant pairs file

　　## ---------------------------------------------------

　　## Adding coordinates for probes and genes from the provided data

　　## Joining, by="GeneID"

　　## Joining, by="Probe"

　　## Saving file: result/getPair.hypo.pairs.significant.withmotif.csv

　　这一步里，我们需要输入我们的MAE对象，筛选出prob-gene对中的probes，并且还需要设置motif最低95%置信区间OR阈值以及最少几个probe进行富集motif。

　　#展示一下名称

　　names(enriched.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是通过哪些探针进行富集的

　　enriched.motif[[1]]

　　## [1] "cg13067635" "cg24873093" "cg25210796" "cg18345456" "cg14094320"

　　## [6] "cg17901924" "cg22736624" "cg12213388" "cg13305336" "cg26607897"

　　## [11] "cg11718886" "cg24593832" "cg25729466" "cg24446417"

　　除此之外，ELMER包还提供了两个表：“

　　getMotif.hypo.enriched.motifs.rda” “

　　getMotif.hypo.motif.enrichment.csv”，大家可以在result目录看一下。

　　五、确定调控作用转录因子

　　这一步，ELMER包对motif以及上游转录因子的关系对进行筛选，从而得到具有调控作用的转录因子。如果某一特定亚组中基因的enhancer发生改变，其上有调控的转录因子同样也会发生改变。基于此，ELMER包将上一步我们通过富集分析找到富集的motif以及对应motif的转录因子进行分析，筛选出对应转录水平发生改变的转录因子。

　　同样地，ELMER提供了两种方式：

　　（1）unsupervised，即将存在相同motif的远端探针（distal probe）根据甲基化水平分为高甲基化组（一般取前20%）和低甲基化组（后20%），比较其两组间对应TF的表达值是否存在差异。

　　（2）Supervised，即对疾病组vs.对照组相同motif对应TF表达值进行分析。接下来筛选出前5% P value最小的TFs，并视为潜在的上游调控转录因子。

　　##找到调控motif对应的转录因子

　　TF

　　group.col="definition",

　　group1="Primary solid Tumor",

　　group2="Solid Tissue Normal",

　　mode="unsupervised",

　　enriched.motif=enriched.motif,

　　dir.out="result",

　　cores=1,

　　label="hypo")

　　## Selecting U (unmethylated) and M (methylated) groups. Each groups has 20% of samples

　　## -------------------------------------------------------------------------------------------------------------------

　　## ** Downloading TF list from Lambert, Samuel A., et al. The human transcription factors. Cell 172.4 (2021): 650-665.

　　## -------------------------------------------------------------------------------------------------------------------

　　## Accessing TF families from TFClass database to indentify known potential TF

　　## Retrieving TFClass family classification from ELMER.data.

　　## Retrieving TFClass subfamily classification from ELMER.data.

　　## Calculating the average methylation at all motif-adjacent probes

　　#展示一下结果

　　head(TF)

　　## motif top.potential.TF.family

　　## HXA9_HUMAN.H11MO.0.B HXA9_HUMAN.H11MO.0.B POU6F2

　　## FOXH1_HUMAN.H11MO.0.A FOXH1_HUMAN.H11MO.0.A FOXE1

　　## PO5F1_HUMAN.H11MO.1.A PO5F1_HUMAN.H11MO.1.A POU6F2

　　## HXB13_HUMAN.H11MO.0.A HXB13_HUMAN.H11MO.0.A HOXD11

　　## IRF9_HUMAN.H11MO.0.C IRF9_HUMAN.H11MO.0.C IRF6

　　## SOX10_HUMAN.H11MO.1.A SOX10_HUMAN.H11MO.1.A SOX2

　　## top.potential.TF.subfamily

　　## HXA9_HUMAN.H11MO.0.B HOXD12

　　## FOXH1_HUMAN.H11MO.0.A FOXH1

　　## PO5F1_HUMAN.H11MO.1.A

　　## HXB13_HUMAN.H11MO.0.A HOXD11

　　## IRF9_HUMAN.H11MO.0.C

　　## SOX10_HUMAN.H11MO.1.A

　　## potential.TF.family

　　## HXA9_HUMAN.H11MO.0.B POU6F2;HOXD12

　　## FOXH1_HUMAN.H11MO.0.A FOXE1;FOXN1;FOXH1;FOXD1

　　## PO5F1_HUMAN.H11MO.1.A POU6F2;HOXA10

　　## HXB13_HUMAN.H11MO.0.A HOXD11;HOXD12;HOXD13;HOXD10;HOXA10;POU6F2;HOXA11

　　## IRF9_HUMAN.H11MO.0.C IRF6

　　## SOX10_HUMAN.H11MO.1.A SOX2;SOX21;SOX15;SOX6

　　## potential.TF.subfamily

　　## HXA9_HUMAN.H11MO.0.B HOXD12

　　## FOXH1_HUMAN.H11MO.0.A FOXH1

　　## PO5F1_HUMAN.H11MO.1.A

　　## HXB13_HUMAN.H11MO.0.A HOXD11;HOXD12;HOXD13;HOXD10;HOXA10;HOXA11

　　## IRF9_HUMAN.H11MO.0.C

　　## SOX10_HUMAN.H11MO.1.A

　　## top_5percent_TFs

　　## HXA9_HUMAN.H11MO.0.B ZFP64;ZNF74;ZNF326;PATZ1;SOX2;OTX1;BCL11A;AHCTF1;ZNF273;TFAP4;ZNF558;PAX6;ZBTB33;SIX4;ZIC2;TFCP2;ZNF77;E2F6;TGIF2;ZNF343;ZBTB12;TCF20;POU6F2;ZNF786;MYNN;DMRT2;SOX12;DLX2;TP63;TBPL1;ZNF18;SP3;ZNF639;SAFB;ZNF3;DLX6;ZNF519;ZBTB39;RFX7;CREB3L4;ZNF292;ZNF669;DLX5;VAX2;ZNF384;ZNF687;HOXD12;SIX1;NEUROG2;ZNF26;SIM2;NEUROD2;ZNF692;MAZ;DMRT1;FOXM1;DMRT3;DLX1;AEBP2

　　## FOXH1_HUMAN.H11MO.0.A SOX2;OTX1;ZFP64;TP63;KLF5;NFE2L2;ZIC2;ZNF326;FOXE1;SOX21;TFAP4;BCL11A;ZNF77;TBPL1;DMRT2;ZNF558;PLAG1;DMRT3;DLX5;ZNF18;TCF20;ZNF273;POU6F2;SOX6;ZNF589;FOXN1;AHCTF1;OTX2;MYCL;SP3;ZBTB33;DLX6;ZNF3;FOXH1;ZNF786;AEBP2;ZNF292;PAX6;ZNF74;ZNF519;GTF2IRD1;ZNF131;FOXD1;MYNN;ZNF639;PITX1;HOXD13;GRHL3;ZNF26;ZNF556;SOX15;NRF1;DMTF1;RFX7;PRDM13;ZNF343;SIX4;HOXD12;DMRT1

　　## PO5F1_HUMAN.H11MO.1.A TFAP4;ZFP64;ZNF326;OTX1;BCL11A;ZBTB33;PATZ1;TP63;DMRT2;SOX2;TCF20;ZNF273;HES6;ZNF74;FOXN1;ZNF77;DLX5;POU6F2;KLF5;GTF2IRD1;PAX6;TFAP2B;NKRF;ZIC2;ZNF692;SOX12;ZNF343;TFCP2;AHCTF1;OTX2;HMGA1;TFAP2C;SOX15;CREB3L4;ZNF589;MYNN;HOXA10;FOXH1;NFE2L2;ZBTB39;MYCL;ZNF519;ZNF18;DMRT3;FOXE1;DLX6;ZNF26;ZNF687;ZNF558;MAZ;NRF1;PITX1;GLI4;ZBTB12;RFX7;SAFB;SOX21;SIX1;AEBP2

　　## HXB13_HUMAN.H11MO.0.A BCL11A;ZFP64;ZNF74;OTX1;TP63;PATZ1;HOXD11;ZIC2;ZNF77;ZNF326;ZNF18;PAX6;SOX2;AEBP2;TCF20;SIM2;TFAP4;HOXD12;TGIF2;KLF5;ZNF273;HOXD13;HOXD10;ZNF26;ZNF639;DLX5;HOXA10;POU6F2;ZBTB33;HOXA11;MYNN;FOXM1;DLX6;DMRT2;TBPL1;RFX7;AHCTF1;ZNF384;ZNF786;SOX21;GTF2IRD1;SIX4;ZNF343;ZNF519;TFCP2;ZBTB39;FOXE1;ZNF131;E2F6;DMRT1;ZNF558;MAZ;SOX6;TLX1;NRF1;EMX1;NFE2L2;ZNF556;ZNF692

　　## IRF9_HUMAN.H11MO.0.C ZFP64;OTX1;TFAP4;ZNF326;ZNF74;TCF20;BCL11A;PATZ1;TP63;ZIC2;SOX2;SAFB;ZNF77;DMRT2;POU6F2;ZBTB33;AHCTF1;ZNF18;TFCP2;PAX6;SIM2;SIX1;ZNF292;SREBF1;IRF6;SIX4;ZNF131;ZNF273;HOXD12;TBPL1;ZNF556;TCF3;ZNF26;ZNF3;ZNF7;DMRT3;FOXH1;NKRF;HOXA11;HES6;AEBP2;ZNF519;ZNF558;DMRT1;ZBTB12;ZNF786;PLAG1;ZNF343;PAX9;TFAP2B;KLF5;FOXE1;HMGA1;FOXN1;SOX21;ZNF687;ELK4;DLX5;SREBF2

　　## SOX10_HUMAN.H11MO.1.A TP63;SOX2;KLF5;OTX1;DLX5;BCL11A;DMRT2;ZNF77;TFAP4;NFE2L2;ZNF273;TCF20;MYNN;ZFP64;FOXE1;ZNF74;ZIC2;ZBTB33;ZNF326;ZNF639;SOX21;DLX6;POU6F2;ZNF18;DMRT3;DMRT1;AEBP2;ELK4;PITX1;ZNF398;FOXN1;GRHL3;ZNF519;RFX7;GTF2IRD1;ZNF556;E2F6;GRHL1;ZNF558;TBPL1;ZNF385A;HOXD12;SP3;EMX1;SOX15;AHCTF1;PRDM13;ZNF3;ZBTB39;ZNF343;HOXA10;ZNF786;SOX6;PATZ1;PAX6;HMGA1;MYCL;MAZ;SMAD3

　　好啦，我们的分析就结束了，总结一下，我们通过同一批样本的甲基化数据和转录组数据，首先找到具有差异的远端探针（distal probe），并根据远端探针和对应基因的表达水平构建出probe-gene对。并将probe-gene对中的probe进一步进行motif富集分析，找出潜在富集motif，最后根据富集的motif推测出上游的调控TFs。这一条逻辑链比较复杂，大家可以再翻出上一期再品鉴一下。下一期我们讲介绍ELMER包的可视化以及高分文章上使用ELMER包的示例。

　　好啦，我是风间琉璃，咱们下期见

　　—END—