MethylMix: an R package for identifying DNA methylation-driven genes.

本文参考：识别甲基化驱动的癌症基因

GetData函数，怎么爬取数据：https://www.bioconductor.org/developers/how-to/web-query/

参考download.file函数：https://github.com/MSQ-123/MethylMix/blob/master/R/Download_Preprocess.R

这个包是用纯R写的，可以参考一下其结构：怎么写函数爬取TCGA的data

background

这个包的功能是鉴定潜在的甲基化驱动的癌基因。也就是说，如果一个基因在癌症样本中相比正常样本低表达，且这个基因在癌症样本中被高度甲基化，那么我们认为这个基因就是一个driver gene

driver gene 的差异甲基化值(DM-value)就是该基因的差异甲基化程度，对于混合的多个相同癌症类型的样本，还可以根据它们的DM-value，对该癌症类型进行分亚型的操作。

Run

第一步获取数据：包提供了获取数据的函数，无须自己下载，但速度感人：

Step 1: Automated Downloading from TCGA: DNA methylation datasets and Gene expression Datasets are downloaded automatically from TCGA by supplying the TCGA cancer code. We have provided the functionality to study any of the 33 TCGA cancer sites that are currently available.

以卵巢癌为例：

library(doParallel)
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)

函数会自动捕获27k or 450k的甲基化芯片数据，并进行适当的合并；Download函数下载的数据不能直接拿来分析，需要预处理：预处理包括校正批次效应，缺省值模拟，CpG探针检测位点clustering，缺省值过多的数据直接过滤等操作，出来以后直接是matrix对象，可以直接用来下游分析，很方便。

这里我们用包的内置数据运行一下：

#examples
data(METcancer)
data(METnormal)
data(GEcancer)
head(METcancer[, 1:4])
head(GEcancer)

基本需要的input就是METcancer（癌症样本的甲基化数据）、METnormal（正常样本的甲基化数据）以及GEcancer（癌症样本的基因表达量数据）。核心函数是MethylMix：

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.

解释一下Results里面每个参数的含义：

• MethylationDrivers: Driver genes，利用beta分布模拟出来的驱动基因。
• NrComponents: The number of methylation states found for each driver gene. 正常来说，methylation states有两个：hypo和hyper
• MixtureStates: 每个驱动基因的DM-values，应该是取平均的值
• MethylationStates: 一个矩阵，每个驱动基因在每个样本中的DM-values
• Classifications: 对于每个样本每个基因的methylation states，值为1对应hypo，值为2对应hyper
• Models: Beta分布模型的参数

获得的Results可以进行可视化：选择需要可视化的基因名，

# Plot the most famous methylated gene for glioblastoma
plots <- MethylMix_PlotModel("MGMT", MethylMixResults, METcancer)
plots$MixtureModelPlot

# Plot MGMT also with its normal methylation variation
plots <- MethylMix_PlotModel("MGMT", MethylMixResults, METcancer, METnormal = METnormal)
plots$MixtureModelPlot

# Plot a MethylMix model for another gene
plots <- MethylMix_PlotModel("ZNF217", MethylMixResults, METcancer, METnormal = METnormal)
plots$MixtureModelPlot

Rplot.png

同样地，我们可以结合表达量数据，探究基因的表达量和甲基化水平的负相关关系：

# Also plot the inverse correlation with gene expression (creates two separate
# plots)
plots <- MethylMix_PlotModel("MGMT", MethylMixResults, METcancer, GEcancer, METnormal)
plots$MixtureModelPlot
plots$CorrelationPlot


# Plot all functional and differential genes
for (gene in MethylMixResults$MethylationDrivers) {
  MethylMix_PlotModel(gene, MethylMixResults, METcancer, METnormal = METnormal)
}

Rplot01.png