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)

畫出來的圖如下所示,可以看到，這個模組的基因主要是高甲基化低表達：