FEM是一個整合RANSEQ和DNA甲基化數據的R包，由Andrew E. Teschendorff 和 Zhen Yang 開發、維護。
不多說，下面介紹如何使用FEM整合RANSEQ和DNA甲基化數據分析。
1、安裝、下載FEM
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("FEM")
library(FEM)
2、下載數據adj.m
adj.m數據存儲在網址（https://sourceforge.net/projects/signalentropy/files/?source=navbar）的“hprdAsigH-13Jun12.Rd”上，下載“hprdAsigH-13Jun12.Rd”即可。
“hprdAsigH-13Jun12.Rd”文件包含三個數據："hprdAsigH.m"、"sigHclassA.v"、"sigHclassA2.v"
"hprdAsigH.m"為我們后續分析需要的數據。
“hprdAsigH-13Jun12.Rd”數據也可以通過公眾號bio生物信息后台發生關鍵字“FEM”獲得。
3、准備DNA甲基化數據
DNA甲基化數據取得是beta值，這里我們命名為“beta”，示例圖如下所示：
行名為每個CpG位點的ID，列名為每個樣本的ID。
4、准備DNA甲基化數據對應的表型文件
DNA甲基化數據對於的表型文件，我們命名為group，其示例圖如下所示：
表示為第一個樣本sample1是control，第二個樣本sample2是control，第三個樣本sample3是case，以此類推。beta文件的sample和group是一一對應的。
5、准備RANSEQ基因表達數據
RANSEQ基因表達數據，我們命名為need，其示例如下所示：
行名是每一個基因的entrez gene IDs，列名是每一個樣本名。
6、准備RANSEQ基因表達數據對應的表型文件
RANSEQ基因表達數據對應的表型文件，我們命名為gg, 示例圖如下所示：
表示的是每一個樣本對應的是case還是control。與DNA甲基化的情況一樣，need文件的sample和group是一一對應的。
7、生成DNA差異甲基化統計量
如果是850k，則用以下命令：
statM.o=GenStatM(beta,group,"EPIC")
如果是450K，則用以下命令：
statM.o=GenStatM(beta,group,"450K")
生成的statM.o結果包含三個數據："top"、 "cont"、"avbeta"
"top"是差異甲基化的統計結果,top是一個list，差異甲基化結果一般存儲在top的第一個元素（item）中；
"cont"是差異甲基化分析時構建的case-control；
"avbeta"是DNA甲基化數據；
8、生成差異表達的統計量
使用命令：
statR.o=GenStatR(need,gg)
生成的statR.o結果包含三個數據： "top"、"cont"、"avexp"
與差異甲基化的結果類似， "top"是差異表達的統計結果；
"cont"是差異表達分析時構建的case-control；
"avbeta"是表達數據；
9、整合差異表達和差異甲基化數據
load("/data/chenwenyan/hprdAsigH-13Jun12.Rd")
re=DoIntFEM450k(statM.o,statR.o,hprdAsigH.m,1,1,"avbeta")
解釋一下，statM.o和statR.o分別是步驟7和8產生的結果文件，hprdAsigH.m是步驟2下載的“hprdAsigH-13Jun12.Rd”文件包含的數據,這里我存儲在/data/chenwenyan/路徑下，請讀者們根據各自存儲的路徑自行修改，不要完成照抄我的路徑。
兩個1分別指的是statM.o和statR.o的top數據的第一個文件，即步驟7和8生成的差異甲基化和差異表達結果。
這里需要注意的是，如果你感興趣的分組結果存儲在top數據的第二個元素，則代碼需要改成re=DoIntFEM450k(statM.o,statR.o,hprdAsigH.m,2,2,"avbeta")
10、鑒定甲基化與表達之間存在負相關的基因
DoFEMbi=DoFEMbi(re, nseeds = 100, gamma = 0.5, nMC = 1000, sizeR.v = c(1,100), minsizeOUT = 10, writeOUT = TRUE, nameSTUDY = "TEST", ew.v = NULL)
這里所有參數均可以使用默認值。
輸出的DoFEMbi結果包含以下文件：
這里我們主要關注fem和topmod這兩個元素,分別指的是模塊以及模塊對應的統計數據，如下所示：
11、可視化結果
可視化SGMS2模塊信息：
SGMS2=FemModShow(DoFEMbi$topmod$SGMS2,name="SGMS2", DoFEMbi)
畫出來的圖如下所示,可以看到，這個模塊的基因主要是高甲基化低表達：
