筆記

請注意，此儲存庫僅用於存檔目的。此儲存庫中的獨立 ANCOM 功能不再維護。強烈建議使用者使用我們的 ANCOMBC R 套件中的ancom或ancombc函數。有關錯誤報告，請存取 ANCOMBC 儲存庫。

ANCOM 使用手冊

當前程式碼在橫斷面和縱向資料集中實作 ANCOM，同時允許使用協變數。需要包含以下庫才能運行 R 程式碼：

library( nlme )
library( tidyverse )
library( compositions )
source( " programs/ancom.R " )

使用說明

資料預處理

我們採用ANCOM-II的方法作為預處理步驟，在進行差異豐度分析之前處理不同類型的零。

用法

feature_table_pre_process(feature_table, meta_data, sample_var, group_var = NULL, out_cut = 0.05, zero_cut = 0.90, lib_cut, neg_lb)

論點

feature_table ：表示觀察到的 OTU/SV 表的資料框或矩陣，其中分類單元在行 ( rownames ) 中，樣本在列 ( colnames ) 中。請注意，這是絕對豐度表，請勿將其轉換為相對豐度表（其中列總計等於 1）。
meta_data ：所有感興趣的變數和協變數的資料框或矩陣。
sample_var ：字元。儲存樣本 ID 的列的名稱。
group_var ：字元。組指示器的名稱。 group_var是檢測結構零點和異常值所必需的。不同零點（結構零點、離群點零點、抽樣零點）的定義請參考ANCOM-II。
out_cut 0 到out_cut之間的數值分數。而混合分佈比例大於1 - out_cut的觀測值將被偵測為異常值。
zero_cut ：0 和zero_cut之間的數值分數。
lib_cut ：數字。文庫大小小於lib_cut的樣本不包含在分析中。
neg_lb ：邏輯。 TRUE 表示分類單元將使用其漸近下界被分類為對應實驗組中的結構零。更具體地說， neg_lb = TRUE表示您正在使用 ANCOM-II 第 3.2 節中所述的兩個標準來偵測結構零點；否則， neg_lb = FALSE將只使用 ANCOM-II 第 3.2 節中的方程式 1 來聲明結構零。

價值

feature_table ：預處理OTU表的資料幀。
meta_data ：預處理元資料的資料框。
structure_zeros ：由 0 和 1 組成的矩陣，其中 1 表示該分類單元被識別為對應組中的結構零。

安康主要功能

用法

ANCOM(feature_table, meta_data, struc_zero, main_var, p_adj_method, alpha, adj_formula, rand_formula, lme_control)

論點

feature_table ：表示 OTU/SV 表的資料框，行中包含分類單元 ( rownames )，列中包含樣本 ( colnames )。它可以是feature_table_pre_process的輸出值。請注意，這是絕對豐度表，請勿將其轉換為相對豐度表（其中列總計等於 1）。
meta_data ：變數的資料框。可以是feature_table_pre_process的輸出值。
struc_zero ：由 0 和 1 組成的矩陣，其中 1 表示該分類單元被辨識為對應組中的結構零。可以是feature_table_pre_process的輸出值。
main_var ：字元。主要感興趣變數的名稱。 ANCOM v2.1 目前支援分類main_var 。
p_adjust_method ：字元。指定調整多重比較 p 值的方法。預設值為“BH”（Benjamini-Hochberg 過程）。
alpha ：顯著性等級。預設值為 0.05。
adj_formula ：表示調整公式的字串（請參閱範例）。
rand_formula ：表示lme中隨機效應公式的字串。有關詳細信息，請參閱?lme 。
lme_control ：指定 lme 擬合控制值的清單。有關詳細信息，請參閱?lmeControl 。

ANCOM 內的測試流程圖

流程圖

價值

out ：包含每個分類群和後續列的W統計量的資料框，這些欄位是 OTU 或分類群在一系列截止值（0.9、0.8、0.7 和 0.6）下是否具有差異豐度的邏輯指標。常用的是detected_0.7 。但是，如果您想對結果更加保守（較小的 FDR），則可以選擇detected_0.8或detected_0.9 ，或者如果您想探索更多發現（較大的功效），則可以使用detected_0.6
fig ：火山圖的ggplot物件。

範例

標準分析

檢測不同實驗組間差異豐度的類群（此處為「受試者」）
OTU 表示例：data/moving-pics-table.tsv
元資料範例：data/moving-pics-sample-metadata.tsv

library( readr )
library( tidyverse )

otu_data = read_tsv( " data/moving-pics-table.tsv " , skip = 1 )
otu_id = otu_data $ `feature-id`
otu_data = data.frame ( otu_data [, - 1 ], check.names = FALSE )
rownames( otu_data ) = otu_id

meta_data = read_tsv( " data/moving-pics-sample-metadata.tsv " )[ - 1 , ]
meta_data = meta_data % > % 
  rename( Sample.ID = SampleID )

source( " programs/ancom.R " )

# Step 1: Data preprocessing

feature_table = otu_data ; sample_var = " Sample.ID " ; group_var = NULL
out_cut = 0.05 ; zero_cut = 0.90 ; lib_cut = 1000 ; neg_lb = FALSE
prepro = feature_table_pre_process( feature_table , meta_data , sample_var , group_var , 
                                   out_cut , zero_cut , lib_cut , neg_lb )
feature_table = prepro $ feature_table # Preprocessed feature table
meta_data = prepro $ meta_data # Preprocessed metadata
struc_zero = prepro $ structure_zeros # Structural zero info

# Step 2: ANCOM

main_var = " Subject " ; p_adj_method = " BH " ; alpha = 0.05
adj_formula = NULL ; rand_formula = NULL ; lme_control = NULL

res = ANCOM( feature_table , meta_data , struc_zero , main_var , p_adj_method , 
            alpha , adj_formula , rand_formula , lme_control )

write_csv( res $ out , " outputs/res_moving_pics.csv " )

# Step 3: Volcano Plot

# Number of taxa except structural zeros
n_taxa = ifelse(is.null( struc_zero ), nrow( feature_table ), sum(apply( struc_zero , 1 , sum ) == 0 ))
# Cutoff values for declaring differentially abundant taxa
cut_off = c( 0.9 * ( n_taxa - 1 ), 0.8 * ( n_taxa - 1 ), 0.7 * ( n_taxa - 1 ), 0.6 * ( n_taxa - 1 ))
names( cut_off ) = c( " detected_0.9 " , " detected_0.8 " , " detected_0.7 " , " detected_0.6 " )

# Annotation data
dat_ann = data.frame ( x = min( res $ fig $ data $ x ), y = cut_off [ " detected_0.7 " ], label = " W[0.7] " )

fig = res $ fig +  
  geom_hline( yintercept = cut_off [ " detected_0.7 " ], linetype = " dashed " ) + 
  geom_text( data = dat_ann , aes( x = x , y = y , label = label ), 
            size = 4 , vjust = - 0.5 , hjust = 0 , color = " orange " , parse = TRUE )
fig

針對協變數進行調整

檢測不同實驗組（此處為「受試者」）的差異豐度類群，同時調整其他協變量（此處為「報告的抗生素使用情況」）
OTU 表示例：data/moving-pics-table.tsv
元資料範例：data/moving-pics-sample-metadata.tsv

library( readr )
library( tidyverse )

otu_data = read_tsv( " data/moving-pics-table.tsv " , skip = 1 )
otu_id = otu_data $ `feature-id`
otu_data = data.frame ( otu_data [, - 1 ], check.names = FALSE )
rownames( otu_data ) = otu_id

meta_data = read_tsv( " data/moving-pics-sample-metadata.tsv " )[ - 1 , ]
meta_data = meta_data % > % 
  rename( Sample.ID = SampleID )

source( " programs/ancom.R " )

# Step 1: Data preprocessing

feature_table = otu_data ; sample_var = " Sample.ID " ; group_var = NULL
out_cut = 0.05 ; zero_cut = 0.90 ; lib_cut = 1000 ; neg_lb = FALSE
prepro = feature_table_pre_process( feature_table , meta_data , sample_var , group_var , 
                                   out_cut , zero_cut , lib_cut , neg_lb )
feature_table = prepro $ feature_table # Preprocessed feature table
meta_data = prepro $ meta_data # Preprocessed metadata
struc_zero = prepro $ structure_zeros # Structural zero info

# Step 2: ANCOM

main_var = " Subject " ; p_adj_method = " BH " ; alpha = 0.05
adj_formula = " ReportedAntibioticUsage " ; rand_formula = NULL ; lme_control = NULL

res = ANCOM( feature_table , meta_data , struc_zero , main_var , p_adj_method , 
            alpha , adj_formula , rand_formula , lme_control )

重複測量/縱向分析

隨機截距模型

檢測不同實驗組之間差異豐富的類群（這裡是「交付」），同時考慮隨機效應（這裡是隨機主題效應「studyid」）
OTU 表示例：data/ecam-table-taxa.tsv
元資料範例：data/ecam-sample-metadata.tsv
透過指定group_var來識別結構零。在這裡我們想知道不同層級的delivery之間是否存在一些結構性零點

library( readr )
library( tidyverse )

otu_data = read_tsv( " data/ecam-table-taxa.tsv " , skip = 1 )
otu_id = otu_data $ `feature-id`
otu_data = data.frame ( otu_data [, - 1 ], check.names = FALSE )
rownames( otu_data ) = otu_id

meta_data = read_tsv( " data/ecam-sample-metadata.tsv " )[ - 1 , ]
meta_data = meta_data % > % 
  rename( Sample.ID = `#SampleID` ) % > %
  mutate( month = as.numeric( month ))

source( " programs/ancom.R " )

# Step 1: Data preprocessing

feature_table = otu_data ; sample_var = " Sample.ID " ; group_var = " delivery "
out_cut = 0.05 ; zero_cut = 0.90 ; lib_cut = 0 ; neg_lb = TRUE
prepro = feature_table_pre_process( feature_table , meta_data , sample_var , group_var , 
                                   out_cut , zero_cut , lib_cut , neg_lb )
feature_table = prepro $ feature_table # Preprocessed feature table
meta_data = prepro $ meta_data # Preprocessed metadata
struc_zero = prepro $ structure_zeros # Structural zero info

# Step 2: ANCOM

main_var = " delivery " ; p_adj_method = " BH " ; alpha = 0.05
adj_formula = NULL ; rand_formula = " ~ 1 | studyid "
lme_control = list ( maxIter = 100 , msMaxIter = 100 , opt = " optim " )

res = ANCOM( feature_table , meta_data , struc_zero , main_var , p_adj_method , 
            alpha , adj_formula , rand_formula , lme_control )

write_csv( res $ out , " outputs/res_ecam.csv " )

# Step 3: Volcano Plot

# Number of taxa except structural zeros
n_taxa = ifelse(is.null( struc_zero ), nrow( feature_table ), sum(apply( struc_zero , 1 , sum ) == 0 ))
# Cutoff values for declaring differentially abundant taxa
cut_off = c( 0.9 * ( n_taxa - 1 ), 0.8 * ( n_taxa - 1 ), 0.7 * ( n_taxa - 1 ), 0.6 * ( n_taxa - 1 ))
names( cut_off ) = c( " detected_0.9 " , " detected_0.8 " , " detected_0.7 " , " detected_0.6 " )

# Annotation data
dat_ann = data.frame ( x = min( res $ fig $ data $ x ), y = cut_off [ " detected_0.7 " ], label = " W[0.7] " )

fig = res $ fig + 
  geom_hline( yintercept = cut_off [ " detected_0.7 " ], linetype = " dashed " ) + 
  geom_text( data = dat_ann , aes( x = x , y = y , label = label ), 
            size = 4 , vjust = - 0.5 , hjust = 0 , color = " orange " , parse = TRUE )
fig

針對其他協變數調整的隨機截距模型

檢測不同實驗組之間差異豐富的類群（這裡是“交付”），同時考慮隨機效應（這裡是隨機主題效應“studyid”）並調整其他協變量（這裡是“月份”）
OTU 表示例：data/ecam-table-taxa.tsv
元資料範例：data/ecam-sample-metadata.tsv
透過指定group_var來識別結構零。在這裡我們想知道不同層級的delivery之間是否存在一些結構性零點

library( readr )
library( tidyverse )

otu_data = read_tsv( " data/ecam-table-taxa.tsv " , skip = 1 )
otu_id = otu_data $ `feature-id`
otu_data = data.frame ( otu_data [, - 1 ], check.names = FALSE )
rownames( otu_data ) = otu_id

meta_data = read_tsv( " data/ecam-sample-metadata.tsv " )[ - 1 , ]
meta_data = meta_data % > % 
  rename( Sample.ID = `#SampleID` ) % > %
  mutate( month = as.numeric( month ))

source( " programs/ancom.R " )

# Step 1: Data preprocessing

feature_table = otu_data ; sample_var = " Sample.ID " ; group_var = " delivery "
out_cut = 0.05 ; zero_cut = 0.90 ; lib_cut = 0 ; neg_lb = TRUE
prepro = feature_table_pre_process( feature_table , meta_data , sample_var , group_var , 
                                   out_cut , zero_cut , lib_cut , neg_lb )
feature_table = prepro $ feature_table # Preprocessed feature table
meta_data = prepro $ meta_data # Preprocessed metadata
struc_zero = prepro $ structure_zeros # Structural zero info

# Step 2: ANCOM

main_var = " delivery " ; p_adj_method = " BH " ; alpha = 0.05
adj_formula = " month " ; rand_formula = " ~ 1 | studyid "
lme_control = list ( maxIter = 100 , msMaxIter = 100 , opt = " optim " )

res = ANCOM( feature_table , meta_data , struc_zero , main_var , p_adj_method , 
            alpha , adj_formula , rand_formula , lme_control )

隨機係數/斜率模型

檢測不同實驗組之間差異豐富的類群（這裡是“交付”），同時考慮隨機效應（這裡是隨機時間效應“月”和隨機受試者效應“studyid”）並調整其他協變量（這裡是“月”））
OTU 表示例：data/ecam-table-taxa.tsv
元資料範例：data/ecam-sample-metadata.tsv
透過指定group_var來識別結構零。在這裡我們想知道不同層級的delivery之間是否存在一些結構性零點

library( readr )
library( tidyverse )

otu_data = read_tsv( " data/ecam-table-taxa.tsv " , skip = 1 )
otu_id = otu_data $ `feature-id`
otu_data = data.frame ( otu_data [, - 1 ], check.names = FALSE )
rownames( otu_data ) = otu_id

meta_data = read_tsv( " data/ecam-sample-metadata.tsv " )[ - 1 , ]
meta_data = meta_data % > % 
  rename( Sample.ID = `#SampleID` ) % > %
  mutate( month = as.numeric( month ))

source( " programs/ancom.R " )

# Step 1: Data preprocessing

feature_table = otu_data ; sample_var = " Sample.ID " ; group_var = " delivery "
out_cut = 0.05 ; zero_cut = 0.90 ; lib_cut = 0 ; neg_lb = TRUE
prepro = feature_table_pre_process( feature_table , meta_data , sample_var , group_var , 
                                   out_cut , zero_cut , lib_cut , neg_lb )
feature_table = prepro $ feature_table # Preprocessed feature table
meta_data = prepro $ meta_data # Preprocessed metadata
struc_zero = prepro $ structure_zeros # Structural zero info

# Step 2: ANCOM

main_var = " delivery " ; p_adj_method = " BH " ; alpha = 0.05
adj_formula = " month " ; rand_formula = " ~ month | studyid "
lme_control = list ( maxIter = 100 , msMaxIter = 100 , opt = " optim " )

res = ANCOM( feature_table , meta_data , struc_zero , main_var , p_adj_method , 
            alpha , adj_formula , rand_formula , lme_control )