概要

本稿では、Metaboanalyst の R パッケージを用いたメタボロームデータ解析の流れを紹介する。

目的

本稿では、IBD multiomics database に登録されているメタボロームデータの一部をデモデータとして用い、CD (クローン病患者) 群と nonIBD (非炎症性腸疾患患者) 群の糞便中に含まれる菌叢代謝物を比較する。今回の代謝物は LC-MS スペクトルから定量されたが、NMR や MS で測定した場合も本稿の解析が可能である。

R コマンド

ライブラリとデータの読み込み

まず MetaboAnalystR を読み込む。

Rコンソール

library("MetaboAnalystR")

次に mSetオブジェクトを作る。MetaboanalystR では、mSet に次々と属性を加えていくことで解析が進行する。

Rコンソール

mSet<-InitDataObjects("pktable", "stat", FALSE)

メタボロームデータを csv ファイルで読み込む。

Rコンソール

mSet<-Read.TextData(mSet, "MetabolomeData.csv", "colu", "disc");

csv ファイルのフォーマットは以下のとおり。
一行目に被験者の ID、二行目に被験者の診断名、三行目以降に代謝物の LC-MS スペクトル強度が記載されている。
詳しくは公式HP を参考にすると良い。

MetabolomeData.csv

Metabolite,MSM9VZHX,MSM9VZJ3,MSMA26AL,PSMA265T,MSM9VZHF,HSMA33JH
diagnosis,nonIBD,CD,CD,CD,nonIBD,CD
1-3-7-trimethylurate,264995.0,7268.0,165010.0,183495.0,4448779.0,171060.0
1-methylguanine,81422.0,156158.0,103739.0,53031.0,126482.0,50352.0
1-methylguanosine,8369.0,17360.0,9266.0,2041.0,6577.0,1725.0
1-methylhistamine,118680.0,178849.0,1271.0,0.0,12265.0,569462.0
1-methylnicotinamide,505048.0,212683.0,1786.0,0.0,4796.0,82993.0

欠損値のフィルタリングと推定

RemoveMissingPercent では、全サンプルの何 % で欠損値だった代謝物をフィルタリングするのかを決める。以下の例では 50% 以上のサンプルで取れてこなかった代謝物が除去される。
ImputeVarでは、欠損値を推定する。method ="min" の場合、全サンプルの最低値に置換される。

Rコンソール

mSet<-SanityCheckData(mSet) # データの状態をチェック
mSet<-RemoveMissingPercent(mSet, percent=0.5)
mSet<-ImputeVar(mSet, method="min")
mSet<-PreparePrenormData(mSet)

サンプルの標準化

SumNorm（組成比変換）、LogNorm（対数変換）、AutoNorm（平均値との差を標準偏差で割る）という3種類のオプションがある。標準化の方法によって解析結果は異なるため、注意深く検討する必要がある。

Rコンソール

mSet<-Normalization(mSet, "SumNorm", "LogNorm", "AutoNorm", ratio=FALSE, ratioNum=20)
mSet<-PlotNormSummary(mSet, "norm_0_", "png", 72, width=NA)
mSet<-PlotSampleNormSummary(mSet, "snorm_0_", "png", 72, width=NA)

標準化結果は、norm_0_dpi72.png と snorm_0_dpi72.png にされる。標準化後の分布が正規分布に近いほど、t 検定の妥当性が上がる。

PCA（主成分分析）

以下のコマンドで PCA を実施する。

Rコンソール

mSet<-PCA.Anal(mSet)
mSet<-PlotPCAPairSummary(mSet, "pca_pair_0_", "png", 72, width=NA, 5)
mSet<-PlotPCAScree(mSet, "pca_scree_0_", "png", 72, width=NA, 5)
mSet<-PlotPCA2DScore(mSet, "pca_score2d_0_", "png", 72, width=NA, 1,2,0.95,0,0)
mSet<-PlotPCALoading(mSet, "pca_loading_0_", "png", 72, width=NA, 1,2);
mSet<-PlotPCABiplot(mSet, "pca_biplot_0_", "png", 72, width=NA, 1,2)
mSet<-PlotPCA3DLoading(mSet, "pca_loading3d_0_", "json", 1,2,3)

結果は、pca_pair_0_dpi72.png に表示される。

各軸の寄与度は、pca_scree_0_dpi72.png の青線に表示される。緑色は累積寄与度を表す。

PLS および PLS-DA

PLS は、PCA で代謝物の違いが見られなかった場合に実施されることがある。

Rコンソール

# PLS
mSet<-PLSR.Anal(mSet, reg=TRUE)
mSet<-PlotPLSPairSummary(mSet, "pls_pair_0_", "png", 72, width=NA, 5)
mSet<-PlotPLS2DScore(mSet, "pls_score2d_0_", "png", 72, width=NA, 1,2,0.95,0,0)
mSet<-PlotPLS3DScoreImg(mSet, "pls_score3d_0_", "png", 72, width=NA, 1,2,3, 40)
mSet<-PlotPLSLoading(mSet, "pls_loading_0_", "png", 72, width=NA, 1, 2);
mSet<-PlotPLS3DLoading(mSet, "pls_loading3d_0_", "json", 1,2,3)
# PLS-DA
mSet<-PLSDA.CV(mSet, "L",5, "Q2")
mSet<-PlotPLS.Classification(mSet, "pls_cv_0_", "png", 72, width=NA)
mSet<-PlotPLS.Imp(mSet, "pls_imp_0_", "png", 72, width=NA, "vip", "Comp. 1", 15,FALSE)

PLS の結果は、pls_score2d_0_dpi72.png や pls_score3d_0_dpi72.png に表示される。各軸は、CD 群か nonIBD 群かという群情報からの寄与も受けるため、PCA よりも分類されやすい。

PLS-DA では、CD 群か nonIBD 群かを区別する際の重要度を表す VIP score が計算される。pls_imp_0_dpi72.png から、どちらの群でどのような代謝物が増加あるいは減少しているのかが分かる。

Fold change

Fold change (FC) は、2 群間で代謝物を比較した際の、差の大きさを表す。

Rコンソール

mSet<-FC.Anal.unpaired(mSet, 2.0, 0)
mSet<-PlotFC(mSet, "fc_0_", "png", 72, width=NA)

結果は fc_0_dpi72に表示される。 $Log_2(FC)$ が 1 以上の代謝物はピンク色でプロットされる。

t 検定

標準化した結果、正規分布に近づけば t 検定を実施することができる。

Rコンソール

mSet<-Ttests.Anal(mSet, F, 0.05, FALSE, TRUE, FALSE)
mSet<-PlotTT(mSet, "tt_0_", "png", 72, width=NA)

結果は tt_0_dpi72.png に表示される。

Random Forest 分析

PLS-DA と同様に、CD 群と nonIBD 群を区別する際に重要な代謝物を選定することができる。

Rコンソール

mSet<-RF.Anal(mSet, 500,7,1)
mSet<-PlotRF.Classify(mSet, "rf_cls_0_", "png", 72, width=NA)
mSet<-PlotRF.VIP(mSet, "rf_imp_0_", "png", 72, width=NA)
mSet<-PlotRF.Outlier(mSet, "rf_outlier_0_", "png", 72, width=NA)

結果はrf_imp_0_dpi72.pngに表示される。MeanDecreaseAccuracyは、その代謝物を除いたときに、CD 群と nonIBD 群の分類精度がどれくらい落ちるのかを表す指標である。VIP score と同様に、この値が大きいほど 2 群に分類する際に重要な代謝物といえる。

Streamlit を用いた欠損値処理と標準化条件の検討

メタボローム解析の結果は、欠損値の処理方法および標準化条件によって異なることがある。
ここでは、欠損値処理と標準化条件の検討方法を提案する。

網羅的解析の実施

以下の script を実行し、欠損値処理と標準化を変えながら網羅的に解析する。
具体的には、RemoveMissingPercent に 0.1, 0.5, 0.9 のいずれかを代入し、Normalization に SumNorm, LogNorm, AutoNorm または NULL を代入する。
今回の場合、RemoveMissingPercent で 3 通り、Normalization で $2^3$ 通り、合計 24 通りのパラメータで解析する。
実行すると、各代謝物に対して Fold change (FC)、 Q value (t検定)、 VIP score (PLS-DA)、MDA (Random Forest 解析) が計算される。
成果物は results.csv に出力される。

Metaboanalyst.R

library("MetaboAnalystR")
rm( list=ls(all=TRUE) )

df <- data.frame()

file <- paste0('MetabolomeData.csv')

for (p in c(0.1, 0.5, 0.9)){
    for (n1 in c('NULL','SumNorm')){
        for (n2 in c('NULL','LogNorm')){    
            for (n3 in c('NULL','AutoNorm')){
                mSet<-InitDataObjects("pktable", "stat", FALSE)
                mSet<-Read.TextData(mSet, file, "colu", "disc");
                mSet<-SanityCheckData(mSet)

                mSet<-RemoveMissingPercent(mSet, percent=p)
                mSet<-ImputeVar(mSet, method="min")
                mSet<-PreparePrenormData(mSet)

                mSet<-Normalization(mSet, n1, n2, n3, ratio=FALSE, ratioNum=20)

                mSet<-Volcano.Anal(mSet, FALSE, 2.0, 0, 0.75,F, 0.1, TRUE, "raw")
                mSet<-PLSR.Anal(mSet, reg=TRUE)
                mSet<-PLSDA.CV(mSet, "L",5, "Q2")
                mSet<-RF.Anal(mSet, 1000,7,1)

                results = cbind (as.data.frame(mSet$analSet$volcano$fc.log)[mSet$dataSet$orig.var.nms,],
                    as.data.frame(mSet$analSet$volcano$p.log)[mSet$dataSet$orig.var.nms,],
                    as.data.frame(mSet$analSet$plsda$vip.mat[,'Comp. 1'])[mSet$dataSet$orig.var.nms,],
                    as.data.frame(mSet$analSet$rf.sigmat)[mSet$dataSet$orig.var.nms,],
                    as.data.frame(mSet$dataSet$orig.var.nms))
                colnames(results) <- c('FC', 'Q_value','VIP','RF_MDA', 'metabolites')

                results <- transform(results,Normalization_method = rep(paste(n1,n2,n3), length(mSet$dataSet$orig.var.nms)))
                results <- transform(results,RemoveMissingPercent = rep(p, length(mSet$dataSet$orig.var.nms)))
                rownames(results) <- NULL

                df <- rbind(df, results)

                rm(mSet)
            }   
        }   
    }
}

write.csv(df, paste0('/results.csv'), fileEncoding = "CP932")

Streamlit を用いた条件検討

Streamlit を用いて、24 通りの解析結果をまとめて表示する。
streamlit run ファイル名.py で、ブラウザにアプリケーションが表示される。Streamlit については、こちらを参考にされるとよい。

st_condition_search.py


import os
import streamlit as st
import altair as alt
import pandas as pd

cwd = os.getcwd()
path = ['results.csv']
file = os.path.join(cwd, *path)

df = pd.read_table(file, sep=',' ,header=0)

axis_list = ['FC', 'Q_value', 'VIP', 'RF_MDA']
metabo_list = sorted(set(list(df['metabolites'])))
RemoveMissingPercent_list = sorted(set(list(df['RemoveMissingPercent'])))
Normalization_method_list = sorted(set(list(df['Normalization_method'])))

x_axis = st.selectbox('Select x axis',  axis_list)
y_axis = st.selectbox('Select y axis',  axis_list)
metabolites = st.multiselect('Select metabolites',metabo_list)
RemoveMissingPercent = st.multiselect('Select RemoveMissingPercent', RemoveMissingPercent_list)
Normalization_method = st.multiselect('Select normalization method', Normalization_method_list)

def select_all_elements(l,ALL):
    if not l:
        l = ALL
    else:
        pass
    return(l)

metabolites = select_all_elements(metabolites,metabo_list)
RemoveMissingPercent = select_all_elements(RemoveMissingPercent,RemoveMissingPercent_list)
Normalization_method = select_all_elements(Normalization_method,Normalization_method_list)

df = df[(df['metabolites'].isin(metabolites))&
        (df['RemoveMissingPercent'].isin(RemoveMissingPercent))&
        (df['Normalization_method'].isin(Normalization_method))] 

plot = alt.Chart(df).mark_circle(size=40).encode(
    x=x_axis,
    y=y_axis,
    tooltip=[x_axis, y_axis, 'metabolites', 'Normalization_method', 'RemoveMissingPercent'],
    color='Normalization_method',
    opacity='RemoveMissingPercent'
).properties(
    width=700,
    height=500
).interactive()

st.write(plot)