生物信息學/MethylMix

MethylMix

MethylMix是一種用於識別疾病的高甲基化和低甲基化基因的算法。MethylMix基於Beta混合模型來識別甲基化狀態，並將其與正常的DNA甲基化狀態進行比較。MethylMix使用一種新的統計量，差異甲基化值或DM值，定義為甲基化狀態與正常甲基化狀態之間的差異。最後，匹配的基因表達數據用於識別除差異外，功能性甲基化狀態，重點關注影響基因表達的甲基化變化。

計算基因甲基化化水平和基因表達相關性，找出顯著相關的基因。
對樣品甲基化程度構建混合模型( mixture model )。
Wilcoxon rank test計算正常組和處理組的差異甲基化。
對計算的p值進行矯正。

MethylMix軟件官網：https://bioconductor.org/packages/release/bioc/html/MethylMix.html

使用手冊

https://bioconductor.org/packages/release/bioc/vignettes/MethylMix/inst/doc/vignettes.html

安裝MethylMix

在R語言IDE命令行中運行如下代碼進行安裝：

# 检查是否安装BiocManager软件包
if (!requireNamespace("BiocManager", quietly=TRUE))
  install.packages("BiocManager")
  
# 安装MethylMix扩展包
BiocManager::install("MethylMix")

MethylMix包的函數與參數

用於下載和整理TCGA數據的函數

Download_DNAmethylation：下載TCGA指定癌症類型的甲基化數據
- Preprocess_DNAmethylation：處理Download_DNAmethylation函數下載的甲基化數據
Download_GeneExpression：下載TCGA指定癌症類型的基因表達數據
- Preprocess_GeneExpression：處理Download_GeneExpressio函數下載的基因表達數據
GetData：下載和整理基因表達和甲基化數據

**Download_*和Getdata**參數
`CancerSite`	長度為1的字符類型，包含TCGA癌症編碼。
`TargetDirectory`	字符類型，用於創建儲存下載文件的目錄名。
`downloadData`	邏輯值，指示是否應下載數據（默認值：TRUE）。如果為false，則返回所需數據的URL。

Download_*返回值

列出 27k 和 450k 甲基化數據的下載文件路徑。

Getdata返回值

將在目標目錄中創建以下文件：

gdac：包含從 TCGA 下載的原始數據的文件夾。
MET_CancerSite_Processed.rds：在 CpG 位點級別（未聚類）處理的甲基化數據。
GE_CancerSite_Processed.rds：處理過的基因表達數據。
data_CancerSite.rds：包含基因表達和甲基化數據的列表。甲基化數據在基因水平上聚集和呈現。包含從 CpG 位點到基因映射的矩陣。

Preprocess_*參數

`CancerSite`	長度為1的字符類型，包含TCGA癌症編碼。
`MAdirectories`	帶有下載數據目錄的字符向量。它可以是 Download_DNAmethylation 函數返回的對象。
`MissingValueThresholdGene`	每個基因缺失值的閾值。去除 NA 百分比大於此閾值的基因。默認值為 0.3。
`MissingValueThresholdSample`	每個樣本缺失值的閾值。去除 NA 百分比大於此閾值的樣本。默認值為 0.1。
`METdirectories`	帶有下載數據目錄的字符向量。它可以是 Download_DNAmethylation 函數返回的對象。
`MissingValueThreshold`	刪除具有缺失值的樣本或基因的閾值。默認值為 0.2

Preprocess_*返回值

列出癌症和正常樣本的預處理數據矩陣。

`ClusterProbes：`其他來源數據整理注釋

此函數使用 Illumina 甲基化芯片的注釋將每個探針映射到一個基因。然後，對於每個基因，它使用層次聚類和 Pearson 相關作為距離和完整連鎖對其所有 CpG 位點進行聚類。如果提供了正常樣本的數據，則僅使用癌症和正常樣本之間的重疊探針。去除帶有 SNP 的探針。如果用戶註冊了一個並行結構，這個函數準備並行運行，否則它會順序運行。此函數還整理樣本名稱，將它們轉換為 12 位數字格式。

參數


`MET_Cancer`	癌症樣本的甲基化數據矩陣。
`MET_Normal`	正常樣本的甲基化數據矩陣。
`CorThreshold`	分類的相關閾值。

返回值

列出聚類數據集以及探針和基因之間的映射。

內置示例數據：METcancer、METnormal和GEcancer

> data(METcancer)
> METcancer[1:4,1:4]
       TCGA.02.0001 TCGA.02.0003 TCGA.02.0006 TCGA.02.0007
ERBB2      0.066976     0.056493      0.16083     0.059322
FAAH       0.251890     0.193620      0.29216     0.856140
FNDC3B     0.691420     0.563600      0.54231     0.306230
FOXD1      0.084173     0.091946      0.19121     0.752810
> data(METnormal)
> METnormal[1:4,1:4]
       X4321207019_B X4321207026_C X4321207042_E X4321207042_K
ERBB2        0.14834      0.105520      0.088934      0.084731
FAAH         0.21961      0.279150      0.196890      0.220890
FNDC3B       0.76382      0.784840      0.796010      0.752130
FOXD1        0.13591      0.074815      0.085960      0.088573
> data(GEcancer)
> GEcancer[1:4,1:4]
       TCGA.02.0001 TCGA.02.0003 TCGA.02.0006 TCGA.02.0007
ERBB2     -0.191840     -1.24440    -0.026743     1.524600
FAAH      -0.075087     -0.70167    -0.580130    -1.237400
FNDC3B     0.033780     -0.64668    -0.356230    -0.050206
FOXD1      0.102790      0.11901     0.155030    -1.002300

`MethylMix`：癌症中 DNA 甲基化數據的混合模型

MethylMix 通過模擬癌症與正常的 DNA 甲基化數據並尋找同質亞群來識別 DNA 甲基化驅動基因。此外，匹配的基因表達數據（例如來自芯片技術或 RNA 測序）用於通過要求特定基因的甲基化和基因表達之間存在負相關來識別功能性 DNA 甲基化事件。請參閱下面的參考資料。

Gevaert 0. MethylMix: an R package for identifying DNA methylation-driven genes. Bioinformatics (Oxford, England). 2015;31(11):1839-41. doi:10.1093/bioinformatics/btv020.

Gevaert O, Tibshirani R, Plevritis SK. Pancancer analysis of DNA methylation-driven genes using MethylMix. Genome Biology. 2015;16(1):17. doi:10.1186/s13059-014-0579-8.

Pierre-Louis Cedoz, Marcos Prunello, Kevin Brennan, Olivier Gevaert; MethylMix 2.0: an R package for identifying DNA methylation genes. Bioinformatics. doi:10.1093/bioinformatics/bty156.

參數

`METcancer`	癌組織甲基化數據的矩陣，行中有基因，列中有樣本。
`GEcancer`	癌組織的基因表達數據，行中有基因，列中有樣本。
`METnormal`	矩陣與 METcancer 中相同基因的正常甲基化數據。同樣是行中的基因和列中的樣本。樣本不必與癌症數據相匹配。如果此參數為 NULL，MethylMix 將在不與正常樣本進行比較的情況下運行。
`listOfGenes`	帶有待評估基因名稱的向量，名稱必須與 METcancer 的行名稱一致。
`filter`	邏輯指示是否應執行線性回歸以選擇在甲基化和基因表達之間具有顯着線性負相關關係的基因（默認值：TRUE）。
`NoNormalMode`	邏輯指示是否應將癌症樣本中發現的甲基化狀態與正常樣本進行比較（默認值：FALSE）。
`OutputRoot`	存儲 MethylMix 結果對象的路徑。

返回值

`MethylationDrivers`	被 MethylMix 鑑定為轉錄預測和差異甲基化的基因。
`NrComponents`	為每個驅動基因找到的甲基化數目。
`MixtureStates`	每個驅動基因的 DM 值列表。對於給定基因，差異甲基化值（DM 值）定義為癌症樣本的一種混合物成分的甲基化平均值與正常樣本的甲基化平均值之間的差異。
`MethylationStates`	具有所有驅動基因（行）和所有樣本（列）的 DM 值的矩陣。
`Classifications`	對於每個基因，帶有整數的矩陣指示每個癌症樣本被分配到哪個混合物成分。
`Models`	每個驅動基因的 Beta 混合模型參數。

對`MethylMix`的輸出結果進行操作的函數

MethylMix_ModelGeneExpression：使用簡單的線性回歸模型將基因表達建模為基因表達的函數。返回在 DNA 甲基化和基因表達之間具有顯着負線性關聯的基因。
MethylMix_PlotModel：生成圖來表示 MethylMix 的輸出。
MethylMix_Predict：給定一個包含甲基化數據的新數據集，此函數可預測每個新樣本和驅動基因的混合成分。預測基於使用 MethylMix'x 擬合混合模型計算的後驗概率。

參數

MethylMix_ModelGeneExpression
`METcancer`	癌症樣本甲基化數據的矩陣（行中的基因，列中的樣本）。
`GEcancer`	癌症樣本基因表達數據的矩陣（行中的基因，列中的樣本）。
`CovariateData`	矢量（數字或字符）指示要包含在模型中以對其進行調整的協變量。不用於 MethylMix 的標準運行。例如，如果樣本可以來自不同的組織類型，則可以使用它。
MethylMix_PlotModel
`GeneName`	要為其創建 MethylMix 圖的基因的名稱。
`MixtureModelResults`	MethylMix 函數返回的列表。
`METcancer`	癌組織甲基化數據的矩陣，行中有基因，列中有樣本。
`GEcancer`	癌組織的基因表達數據，行中有基因，列中有樣本（可選）。
`METnormal`	矩陣與 METcancer 中相同基因的正常甲基化數據（可選）。同樣是行中的基因和列中的樣本。
MethylMix_Predict
`newBetaValuesMatrix`	具有用於預測的新觀察結果的矩陣，行中的基因/cpg 位點，列中的樣本。雖然這個新矩陣的基因/cpg 位點數量可能與運行 MethylMix 時作為 METcancer 提供的矩陣不同，但基因/cpg 位點的命名應該是相同的。
`MethylMixResult`	MethylMix 的輸出對象

返回值


MethylMix_ModelGeneExpression	帶有基因名稱的向量，在這些基因中，甲基化和基因表達之間存在顯着的線性和負相關關係。
MethylMix_PlotModel	包含 MethylMix 圖的列表、甲基化數據的直方圖 (MixtureModelPlot) 和 DNA 甲基化與基因表達之間的散點圖（CorrelationPlot，如果未提供基因表達數據，則為 NULL）。兩個圖都顯示了已識別的不同混合物成分。
MethylMix_Predict	具有預測（混合成分的索引）、行中的驅動基因、列中的新樣本的矩陣

MethylMix運行示例

一次性下載TCGA基因表達和甲基化數據

library(MethylMix)
cancerSite <- "OV"
targetDirectory <- paste0(getwd(), "/")
GetData(cancerSite, targetDirectory)
### 并行处理
cancerSite <- "OV"
targetDirectory <- paste0(getwd(), "/")

library(doParallel)
cl <- makeCluster(5)
registerDoParallel(cl)
GetData(cancerSite, targetDirectory)
stopCluster(cl)

分別下載TCGA表達和甲基化數據

library(MethylMix)
cancerSite <- "OV"
targetDirectory <- paste0(getwd(), "/")

cl <- makeCluster(5)
registerDoParallel(cl)

# Downloading methylation data
METdirectories <- Download_DNAmethylation(cancerSite, targetDirectory)
# Processing methylation data
METProcessedData <- Preprocess_DNAmethylation(cancerSite, METdirectories)
# Saving methylation processed data
saveRDS(METProcessedData, file = paste0(targetDirectory, "MET_", cancerSite, "_Processed.rds"))

# Downloading gene expression data
GEdirectories <- Download_GeneExpression(cancerSite, targetDirectory)
# Processing gene expression data
GEProcessedData <- Preprocess_GeneExpression(cancerSite, GEdirectories)
# Saving gene expression processed data
saveRDS(GEProcessedData, file = paste0(targetDirectory, "GE_", cancerSite, "_Processed.rds"))

# Clustering probes to genes methylation data
METProcessedData <- readRDS(paste0(targetDirectory, "MET_", cancerSite, "_Processed.rds"))
res <- ClusterProbes(METProcessedData[[1]], METProcessedData[[2]])

# Putting everything together in one file
toSave <- list(METcancer = res[[1]], METnormal = res[[2]], GEcancer = GEProcessedData[[1]],
    GEnormal = GEProcessedData[[2]], ProbeMapping = res$ProbeMapping)
saveRDS(toSave, file = paste0(targetDirectory, "data_", cancerSite, ".rds"))

stopCluster(cl)

其他來源數據整理

library(MethylMix)
METcancer = matrix(data = methylation_data, nrow = nb_of_genes, ncol = nb_of_samples)
METnormal = matrix(data = methylation_data, nrow = nb_of_genes, ncol = nb_of_samples)
GEcancer = matrix(data = expression_data, nrow = nb_of_genes, ncol = nb_of_samples)
ClusterProbes(MET_Cancer, MET_Normal, CorThreshold = 0.4)

`MethylMix`構建模型

> library(MethylMix)
> library(doParallel)
Loading required package: foreach
Loading required package: iterators
Loading required package: parallel
> data(METcancer)
> data(METnormal)
> data(GEcancer)
### 不开启并行处理
> MethylMixResults <- MethylMix(METcancer, GEcancer, METnormal)
Found 251 samples with both methylation and expression data.
Correlating methylation data with gene expression...

Found 9 transcriptionally predictive genes.

Starting Beta mixture modeling.
Running Beta mixture model on 9 genes and on 251 samples.
ERBB2 :  2  components are best.
FAAH :  2  components are best.
FOXD1 :  2  components are best.
ME1 :  2  components are best.
MGMT :  2  components are best.
OAS1 :  2  components are best.
SOX10 :  2  components are best.
TRAF6 :  2  components are best.
ZNF217 :  2  components are best.
Warning message:
executing %dopar% sequentially: no parallel backend registered 
### 并行处理
> library(doParallel)
> cl <- makeCluster(5)
> registerDoParallel(cl)
> MethylMixResults <- MethylMix(METcancer, GEcancer, METnormal)
Found 251 samples with both methylation and expression data.
Correlating methylation data with gene expression...

Found 9 transcriptionally predictive genes.

Starting Beta mixture modeling.
Running Beta mixture model on 9 genes and on 251 samples.
> stopCluster(cl)
### 查看结果
> MethylMixResults$MethylationDrivers
[1] "ERBB2"  "FAAH"   "FOXD1"  "ME1"    "MGMT"   "OAS1"   "SOX10" 
[8] "TRAF6"  "ZNF217"
> MethylMixResults$NrComponents
 ERBB2   FAAH  FOXD1    ME1   MGMT   OAS1  SOX10  TRAF6 ZNF217 
     2      2      2      2      2      2      2      2      2 
> MethylMixResults$MixtureStates
$ERBB2
         [,1]
[1,] 0.000000
[2,] 0.318555

$FAAH
          [,1]
[1,] 0.0000000
[2,] 0.3973068

$FOXD1
          [,1]
[1,] 0.0000000
[2,] 0.2217989

$ME1
          [,1]
[1,] 0.1418296
[2,] 0.5834757

$MGMT
          [,1]
[1,] 0.0000000
[2,] 0.3141419

$OAS1
           [,1]
[1,] -0.3391054
[2,]  0.0000000

$SOX10
          [,1]
[1,] 0.0000000
[2,] 0.2076643

$TRAF6
          [,1]
[1,] 0.0000000
[2,] 0.2452776

$ZNF217
           [,1]
[1,] -0.1724214
[2,]  0.1875210

> MethylMixResults$MethylationStates[, 1:5]
       TCGA.02.0001 TCGA.02.0003 TCGA.02.0006 TCGA.02.0007
ERBB2     0.0000000    0.0000000    0.0000000    0.0000000
FAAH      0.0000000    0.0000000    0.0000000    0.3973068
FOXD1     0.0000000    0.0000000    0.2217989    0.2217989
ME1       0.5834757    0.5834757    0.5834757    0.5834757
MGMT      0.0000000    0.0000000    0.3141419    0.0000000
OAS1     -0.3391054    0.0000000   -0.3391054   -0.3391054
SOX10     0.2076643    0.2076643    0.2076643    0.2076643
TRAF6     0.2452776    0.0000000    0.0000000    0.0000000
ZNF217   -0.1724214   -0.1724214   -0.1724214   -0.1724214
       TCGA.02.0009
ERBB2     0.0000000
FAAH      0.3973068
FOXD1     0.2217989
ME1       0.1418296
MGMT      0.0000000
OAS1      0.0000000
SOX10     0.2076643
TRAF6     0.0000000
ZNF217   -0.1724214
> MethylMixResults$Classifications[, 1:5]
       TCGA.02.0001 TCGA.02.0003 TCGA.02.0006 TCGA.02.0007
ERBB2             1            1            1            1
FAAH              1            1            1            2
FOXD1             1            1            2            2
ME1               2            2            2            2
MGMT              1            1            2            1
OAS1              1            2            1            1
SOX10             2            2            2            2
TRAF6             2            1            1            1
ZNF217            1            1            1            1
       TCGA.02.0009
ERBB2             1
FAAH              2
FOXD1             2
ME1               1
MGMT              1
OAS1              2
SOX10             2
TRAF6             1
ZNF217            1

結果可視化

plots <- MethylMix_PlotModel("MGMT", MethylMixResults, METcancer)
plots$MixtureModelPlot
plots <- MethylMix_PlotModel("MGMT", MethylMixResults, METcancer, METnormal = METnormal)
plots$MixtureModelPlot
plots <- MethylMix_PlotModel("ZNF217", MethylMixResults, METcancer, METnormal = METnormal)
plots$MixtureModelPlot
plots <- MethylMix_PlotModel("MGMT", MethylMixResults, METcancer, GEcancer, METnormal)
plots$MixtureModelPlot
plots$CorrelationPlot
for (gene in MethylMixResults$MethylationDrivers) {
    MethylMix_PlotModel(gene, MethylMixResults, METcancer, METnormal = METnormal)
}