【R コード説明】コヨーテのマイクロバイオーム解析 ~ 1. データの前処理（欠損値の補完）~ #MachineLearning

カナダにおいて、16S rRNA遺伝子アンプリコンシークエンシングを用いて、糞便微生物叢が都市部のコヨーテにおける人為的な食餌（ファストフードの廃棄物など）と不健康の関係、また同様に農村部のコヨーテにおけるタンパク質が豊富な食餌（小型草食動物や偶蹄類など）と健康の関係をどのように媒介するかを調べた研究を見つけました。
また、本研究で活用されたシークエンスデータとメタデータおよびRスクリプトを取得できたので、ここで結果を再現してみることにした。実行コマンドには、説明を付与しますが、これまで扱ったことのない領域についてはAIを活用して記載します。

データの準備

RスクリプトとRDataファイルは下記のリンク先に配架してあります。

fasterq-dumpによってシークエンスデータをダウンロードする

cat SRR_Acc_List.txt | xargs fasterq-dump --split-files --outdir ./download_directory

コード説明（Linux環境）
catコマンドによってAccesion List（IDリスト）を読み込む。
xargsにて、読み込んだ内容をfasterq-dumpに渡す。
--split-filesは、ペアエンド形式への対応のオプション。
--outdirでシークエンスデータの出力先を指定する。

RDataファイルの読み込み

# cy_workspaceを読み込む
load("cy_workspace.RData")

# phyloseqオブジェクトを確認する
phyloseq

# コヨーテのメタデータを確認する
cy.metadata

コード説明
cy_workspace.RDataは、Rの作業環境のことです。作業環境には、ベクトル、行列、データフレーム、リスト、関数などが含まれます。
※Rのセッションを終了するときは、作業環境のイメージを保存することができ、保存したイメージは、次にRを起動したときに自動的に再読み込みされます。

cy_workspaceには、phyloseqオブジェクトとコヨーテのメタデータが含まれています。phyloseqオブジェクトは、マイクロバイオームデータのインポート、保存、分析、グラフィカル表示を容易にするためのクラスとツールのセットです。

具体的に、OTUテーブル、サンプルデータ、分類表、系統樹などのデータが含まれています。また、メタデータには、コヨーテの体重、胴回り、年齢、健康指数、脾臓の重さなどの情報が含まれています。これらのデータを使って、コヨーテの健康状態やマイクロバイオームとの関係を分析することができます。

データの前処理

#### IMPORT AND PROCESS HEALTH METADATA ####
library(MASS) #パッケージの読み込み
library(ggplot2) # パッケージの読み込み

cy.metadata <- read.csv("coyote_metadata.csv") #cy.metadataオブジェクトは、cy_workspace.RDataにあらかじめ含まれているので、このコードは実行不要。

# Convert data to factors, where necessary.
for(i in c("PA.Empty", "PA.Prey", "PA.Veg", "PA.Anthro", "Analysis.Season")){
  cy.metadata[,i] <- as.factor(as.character(cy.metadata[,i]))
  }

# Identify and impute any missing data points using linear regressions on the remaining data.
ggplot(cy.metadata, aes(Mass)) + geom_histogram()
impute.model <- lm(Mass ~ Girth + Cementum.Age + KFI, data=cy.metadata) 

# Model for missing mass.
impute.data <- subset(cy.metadata, is.na(Mass)=="TRUE") 
for(i in 1:nrow(cy.metadata)){
  if(is.na(cy.metadata[i,"Mass"])=="TRUE"){
    cy.metadata[i,"Mass"] <- predict(impute.model, newdata=cy.metadata[i,]) } }

for(i in 1:nrow(cy.metadata)){ # Replace missing normalized spleen mass for coyote 73.
  if(is.na(cy.metadata[i,"Spleen.to.Body.Ratio"])=="TRUE"){
    cy.metadata[i,"Spleen.to.Body.Ratio"] <- cy.metadata[i,"Spleen.Mass"]/cy.metadata[i,"Mass"] } }

ggplot(cy.metadata, aes(Snout.to.Base)) + geom_histogram()
impute.model <- lm(Snout.to.Base ~ Girth + KFI + Mass, data=cy.metadata) # Missing body length
impute.data <- subset(cy.metadata, is.na(Snout.to.Base)=="TRUE") 
for(i in 1:nrow(cy.metadata)){
  if(is.na(cy.metadata[i,"Snout.to.Base"])=="TRUE"){
    cy.metadata[i,"Snout.to.Base"] <- predict(impute.model, newdata=cy.metadata[i,]) } }

ggplot(cy.metadata, aes(KFI)) + geom_histogram()
impute.model <- lm(KFI ~ Mass + Girth + Snout.to.Base, data=cy.metadata) # Missing KFI
impute.data <- subset(cy.metadata, is.na(KFI)=="TRUE") 
for(i in 1:nrow(cy.metadata)){
  if(is.na(cy.metadata[i,"KFI"])=="TRUE"){
    cy.metadata[i,"KFI"] <- predict(impute.model, newdata=cy.metadata[i,]) } }

ggplot(cy.metadata, aes(Cementum.Age)) + geom_histogram()
impute.model <- glm.nb(Cementum.Age ~ Mass + Girth + Snout.to.Base, data=cy.metadata)

impute.data <- subset(cy.metadata, is.na(Cementum.Age)=="TRUE") # Missing age.

for(i in 1:nrow(cy.metadata)){
  if(is.na(cy.metadata[i,"Cementum.Age"])=="TRUE"){
    cy.metadata[i,"Cementum.Age"] <- predict(impute.model, newdata=cy.metadata[i,]) } }

rm(impute.model, impute.data)

コード説明

まず、cy.metadata <- read.csv("coyote_metadata.csv")というコードで、coyote_metadata.csvというファイルを読み込んで、cy.metadataという変数に代入する。
※load()による読み込みで、cy.metadataオブジェクトは最初から用意されているので、このコードの実行は、実際のところ不要です。
次に、for(i in c("PA.Empty", "PA.Prey", "PA.Veg", "PA.Anthro", "Analysis.Season"))というループで、cy.metadataの中の5つの列を、因子型（カテゴリカル変数）に変換しています。これは、これらの列が数値ではなく、カテゴリを表す文字列であるためです。例えば、PA.Emptyは、コヨーテの胃が空どうかを示す列で、YesかNoの値をとります。as.factor(as.character(cy.metadata[,i]))という関数で、文字列から因子型に変換します。
次に、ggplot(cy.metadata, aes(Mass)) + geom_histogram()というコードで、cy.metadataの中のMassという列（コヨーテの体重）のヒストグラムを描いています。このヒストグラムを見ると、いくつかの欠損値（NA）があることがわかります。
そこで、impute.model <- lm(Mass ~ Girth + Cementum.Age + KFI, data=cy.metadata)というコードで、Massという列の欠損値を補完するための線形回帰モデルを作成しています。Massを従属変数とし、Girth（コヨーテの胴回り）、Cementum.Age（コヨーテの年齢）、KFI（腎臓脂肪指数）を説明変数としています。
impute.model <-という部分は、線形モデルの結果をimpute.modelという変数に代入することを意味します。
lm(Mass ~ Girth + Cementum.Age + KFI, data=cy.metadata)という部分は、線形モデルを適合する関数の引数です。ここでは、以下のように指定しています。
Mass ~ Girth + Cementum.Age + KFIという式は、線形モデルの式を表します。ここでは、Massを従属変数（目的変数）とし、Girth, Cementum.Age, KFIを独立変数（説明変数）としています。つまり、Mass = β0 + β1 * Girth + β2 * Cementum.Age + β3 * KFI + εというモデルを仮定しています。ここで、β0, β1, β2, β3は回帰係数、εは誤差項です。
次に、impute.data <- subset(cy.metadata, is.na(Mass)=="TRUE")というコードで、cy.metadataの中から、Massが欠損値である行だけを抽出して、impute.dataという変数に代入しています。is.na(Mass)=="TRUE"という条件で、欠損値であるかどうかを判定しています。
次に、for(i in 1:nrow(cy.metadata))というループで、cy.metadataの各行に対して、以下の処理を行っています。
- if(is.na(cy.metadata[i,"Mass"])=="TRUE")という条件で、Massが欠損値であるかどうかを判定しています。
- もし欠損値であれば、cy.metadata[i,"Mass"] <- predict(impute.model, newdata=cy.metadata[i,])というコードで、先ほど作成したモデルを使って、欠損値を予測して補完しています。predict(impute.model, newdata=cy.metadata[i,])という関数で、モデルと新しいデータ（欠損値の行）を渡して、予測値を得ています。
これにより、cy.metadataのMass列の欠損値が補完されます。
次に、for(i in 1:nrow(cy.metadata))というループで、cy.metadataの各行に対して、以下の処理を行っています。
- if(is.na(cy.metadata[i,"Spleen.to.Body.Ratio"])=="TRUE")という条件で、Spleen.to.Body.Ratio（コヨーテの脾臓の重さと体重の比率）が欠損値であるかどうかを判定しています。
- もし欠損値であれば、cy.metadata[i,"Spleen.to.Body.Ratio"] <- cy.metadata[i,"Spleen.Mass"]/cy.metadata[i,"Mass"]というコードで、Spleen.Mass（コヨーテの脾臓の重さ）とMass（コヨーテの体重）の比で計算して補完しています。
これにより、cy.metadataのSpleen.to.Body.Ratio列の欠損値が補完されます。

以下、Snout.to.BaseについてMassの欠損値補完と同じ手順です。

次に、ggplot(cy.metadata, aes(Snout.to.Base)) + geom_histogram()というコードで、cy.metadataの中のSnout.to.Baseという列（コヨーテの鼻先から尾の付け根までの長さ）のヒストグラムを描いています。このヒストグラムを見ると、いくつかの欠損値（NA）があることがわかります。
そこで、impute.model <- lm(Snout.to.Base ~ Girth + KFI + Mass, data=cy.metadata)というコードで、Snout.to.Baseという列の欠損値を補完するための線形回帰モデルを作成しています。このモデルは、Snout.to.Baseを従属変数とし、Girth（コヨーテの胴回り）、KFI（腎臓脂肪指数）、Mass（コヨーテの体重）を説明変数としています。
次に、impute.data <- subset(cy.metadata, is.na(Snout.to.Base)=="TRUE")というコードで、cy.metadataの中から、Snout.to.Baseが欠損値である行だけを抽出して、impute.dataという変数に代入しています。is.na(Snout.to.Base)=="TRUE"という条件で、欠損値であるかどうかを判定しています。
次に、for(i in 1:nrow(cy.metadata))というループで、cy.metadataの各行に対して、以下の処理を行っています。
- if(is.na(cy.metadata[i,"Snout.to.Base"])=="TRUE")という条件で、Snout.to.Baseが欠損値であるかどうかを判定しています。
- もし欠損値であれば、cy.metadata[i,"Snout.to.Base"] <- predict(impute.model, newdata=cy.metadata[i,])というコードで、先ほど作成した線形回帰モデルを使って、欠損値を予測して補完しています。predict(impute.model, newdata=cy.metadata[i,])という関数で、モデルと新しいデータ（欠損値の行）を渡して、予測値を得ています。
これにより、cy.metadataのSnout.to.Base列の欠損値が補完されます。

以下、KFI（腎臓脂肪指数）の欠損値についてMassと同じ手法の補完方法です。

次に、ggplot(cy.metadata, aes(KFI)) + geom_histogram()というコードで、cy.metadataの中のKFIという列（腎臓脂肪指数）のヒストグラムを描いています。このヒストグラムを見ると、いくつかの欠損値（NA）があることがわかります。
そこで、impute.model <- lm(KFI ~ Mass + Girth + Snout.to.Base, data=cy.metadata)というコードで、KFIという列の欠損値を補完するための線形回帰モデルを作成しています。このモデルは、KFIを従属変数とし、Mass（コヨーテの体重）、Girth（コヨーテの胴回り）、Snout.to.Base（コヨーテの鼻先から尾の付け根までの長さ）を説明変数としています。
次に、impute.data <- subset(cy.metadata, is.na(KFI)=="TRUE")というコードで、cy.metadataの中から、KFIが欠損値である行だけを抽出して、impute.dataという変数に代入しています。is.na(KFI)=="TRUE"という条件で、欠損値であるかどうかを判定しています。
次に、for(i in 1:nrow(cy.metadata))というループで、cy.metadataの各行に対して、以下の処理を行っています。
- if(is.na(cy.metadata[i,"KFI"])=="TRUE")という条件で、KFIが欠損値であるかどうかを判定しています。
- もし欠損値であれば、cy.metadata[i,"KFI"] <- predict(impute.model, newdata=cy.metadata[i,])というコードで、先ほど作成した線形回帰モデルを使って、欠損値を予測して補完しています。predict(impute.model, newdata=cy.metadata[i,])という関数で、モデルと新しいデータ（欠損値の行）を渡して、予測値を得ています。
これにより、cy.metadataのKFI列の欠損値が補完されます。

以下、コヨーテの年齢についての欠損値の補完です。

次に、ggplot(cy.metadata, aes(Cementum.Age)) + geom_histogram()というコードで、cy.metadataの中のCementum.Ageという列（コヨーテの年齢）のヒストグラムを描いています。このヒストグラムを見ると、いくつかの欠損値（NA）があることがわかります。
そこで、impute.model <- glm.nb(Cementum.Age ~ Mass + Girth + Snout.to.Base, data=cy.metadata)というコードで、Cementum.Ageという列の欠損値を補完するための負の二項回帰モデルを作成しています。このモデルは、Cementum.Ageを従属変数とし、Mass（コヨーテの体重）、Girth（コヨーテの胴回り）、Snout.to.Base（コヨーテの鼻先から尾の付け根までの長さ）を説明変数としています。これらの変数は、Cementum.Ageと相関があると考えられるため、選択されています。負の二項回帰モデルは、従属変数が負の整数で、過分散（分散が平均より大きい）である場合に適したモデルです。
次に、impute.data <- subset(cy.metadata, is.na(Cementum.Age)=="TRUE")というコードで、cy.metadataの中から、Cementum.Ageが欠損値である行だけを抽出して、impute.dataという変数に代入しています。is.na(Cementum.Age)=="TRUE"という条件で、欠損値であるかどうかを判定しています。
次に、for(i in 1:nrow(cy.metadata))というループで、cy.metadataの各行に対して、以下の処理を行っています。
- if(is.na(cy.metadata[i,"Cementum.Age"])=="TRUE")という条件で、Cementum.Ageが欠損値であるかどうかを判定しています。
- もし欠損値であれば、cy.metadata[i,"Cementum.Age"] <- predict(impute.model, newdata=cy.metadata[i,])というコードで、先ほど作成した負の二項回帰モデルを使って、欠損値を予測して補完しています。predict(impute.model, newdata=cy.metadata[i,])という関数で、モデルと新しいデータ（欠損値の行）を渡して、予測値を得ています。
これにより、cy.metadataのCementum.Age列の欠損値が補完されます。
最後に、rm(impute.model, impute.data)というコードで、不要になった変数を削除しています。これは、メモリの節約や混乱の防止のために行われます。

〇 Massのヒストグラム

ggplot(cy.metadata, aes(Mass)) + geom_histogram()

〇 Snout.to.Base（鼻先から尾の付け根までの長さ）のヒストグラム

ggplot(cy.metadata, aes(Snout.to.Base)) + geom_histogram()

KFI（腎臓脂肪指数）のヒストグラム

ggplot(cy.metadata, aes(KFI)) + geom_histogram()

コヨーテの年齢のヒストグラム

ggplot(cy.metadata, aes(Cementum.Age)) + geom_histogram()

【R コード説明】コヨーテのマイクロバイオーム解析 ~ 1. データの前処理（欠損値の補完）~