多组学分析MOFA用法

 多组学分析是结合多种组学，如转录组、基因组等，构建基因表达网络，深层次理解分子的调控和因果关系，使用无监督的机器学习方法，并且可以设置使用GPU进行训练

本文是对“Drug-perturbation-based stratification of blood cancer”的复现，旨在说明多组学分析的进行方法

CLL：慢性淋巴细胞白血病

1.1 载入包和数据

#if (!require("BiocManager", quietly = TRUE))
#install.packages("BiocManager")
#BiocManager::install("MOFA2")
#BiocManager::install("MOFAdata")
library(MOFA2)
library(MOFAdata)
library(data.table)
library(ggplot2)
library(tidyverse)
# Load the example data
utils::data("CLL_data")   
# Create the MOFA object
MOFAobject <- create_mofa(CLL_data)
plot_data_overview(MOFAobject)
CLL_metadata <- fread("ftp://ftp.ebi.ac.uk/pub/databases/mofa/cll_vignette/sample_metadata.txt")

cll_data数据示例

Drugs：药物反应的生存率值

methylation：甲基化的M值，M值是一种衡量CpG位点甲基化水平的指标，计算公式为log2(M/U)，其中M代表甲基化信号强度，U代表非甲基化信号强度。甲基化信号强度是一种衡量DNA甲基化水平的方法，通过使用特定的探针或引物来检测甲基化的CpG位点的荧光信号的强度来计算的。甲基化信号强度可以反映甲基化的程度和分布，适用于样本内的特征分析。

mRNA: 标准化的表达值

Mutations: 突变状态

cll_metadata

TTT:Time to Treatment，指从诊断到开始治疗的时间

TTD: Time to Death,指从诊断到死亡的时间

trisomy12: 12号染色体三体,trisomy12是CLL/SLL中最常见的染色体异常之一，约占16%的病例。通常表现为CD38阳性、ZAP70阳性、IgHV未突变、CD49d高表达等特征。

1.2模型参数设置

data_opts <- get_default_data_options(MOFAobject)
model_opts <- get_default_model_options(MOFAobject)
model_opts$num_factors <- 15
train_opts <- get_default_training_options(MOFAobject)
train_opts$convergence_mode <- "slow"
train_opts$seed <- 42


MOFAobject <- prepare_mofa(MOFAobject,
                           data_options = data_opts,
                           model_options = model_opts,
                           training_options = train_opts
)

分为model、train分别设置参数，一般使用默认值即可

1.3 运行模型

# Run the model
outfile = file.path(getwd(),"model.hdf5")
MOFAobject.trained <- run_mofa(MOFAobject, outfile)
#MOFAobject <- readRDS(url("http://ftp.ebi.ac.uk/pub/databases/mofa/cll_vignette/MOFA2_CLL.rds"))

MOFAobject.trained

 exectations包含factor和weight

factor是对变量的不同线性组合，共分成15个，是对患者的划分。

1.4 可视化分析

# Plot the factors
-
plot_factor_cor(MOFAobject.trained)
plot_data_scatter(MOFAobject.trained, 
                  view = "mRNA",
                  factor = 1,  
                  features = 4,
                  sign = "positive",
                  color_by = "IGHV"
) + labs(y="RNA expression")

我们有大量具有较大正权重和负权重的基因。 正值较大的基因在有 IGHV 突变的样本中表达较多，而负值较大的基因在没有 IGHV 突变的样本中表达较多。 我们来验证一下。 函数plot_data_scatter 生成具有最大正权重的前 4 个基因的因子 1 值（x 轴）与表达值（y 轴）的散点图。 样品按 IGHV 状态着色：