LEfSe分析的在線+本地運行的超詳細教程

LEfSe（Linear discriminant analysis Effect Size）一種廣泛用於宏基因組生物高維數據的biomaker鑑別和解釋的工具，通過分類比較、生物學一致性測試和效應大小估計，描述和驗證兩個或多個微生物群落之間的差異物種、基因或途徑（Segata et al, 2011）。本篇以某16S高通量測序數據所得的微生物群落數據為例，展示使用LEfSe尋找組間顯著差異的微生物類群，以及確定它們的重要程度。下文中涉及的示例數據、腳本代碼、示例結果等，已存放至百度盤，可能會用得到。https://pan.baidu.com/s/1ctwyoTqT1ofygMUnTj7CPA

LEfSe通過將具有統計學意義的標準檢驗與表徵生物學一致性和效應相關性的附加檢驗結合起來，確定最有可能解釋類別之間差異的特徵（如物種、進化枝、基因或功能，統稱為biomaker）。類別比較方法通常用於預測biomaker，這些標記由違反類別之間無差異零假設的特徵組成；並通過檢測特徵子集的豐度模式，估計顯著變化的程度；效應大小提供了每個特徵對類別差異貢獻大小的估計。

以下是對LEfSe工作原理的簡單概括，細節部分可參考Segata等（2011）的論文。

Data：輸入數據包含m個樣本（列）的集合，每個樣本由n個變量（行，通常進行行歸一化，紅色代表高值，綠色代表低值）組成；每個樣本都對應了一些類別（即樣本分組），代表正在研究的主要生物學假設，可以為兩組或多組，也可以包含反映類內分組的子類別標籤。

Step1：通過Kruskall-Wallis檢驗分析所有變量，檢驗不同類別中的值是否存在差異分布。

Step2：對於Step1中違反零假設的變量，即Kruskall-Wallis檢驗顯著的結果，繼續使用成對的Wilcoxon檢驗檢驗不同類別內的子類別之間所有成對比較是否與類別水平趨勢顯著一致。

Step3：所得的向量子集用於構建線性判別分析（LDA）模型，利用類間的相對差異對變量進行降維和分類，輸出包含一列與類別有關的特徵列表，這些特徵與類別中的子類別分組保持一致，並根據它們區分類別的效應大小進行排名。

 

因此LEfSe是用於對不同生物學方面的相關性進行排名以及進行進一步研究和分析的有價值的工具。它強調統計意義和生物相關性，能夠在組與組之間尋找具有統計學差異的biomaker。在該方法中引入先驗知識作為監督，解決了傳統上與高維數據挖掘有關的挑戰。

接下來，展示LEfSe工具的使用，以及對結果部分的說明。

LEfSe的輸入文件要求將微生物豐度表與分組信息合併，整理為如下所示的樣式。

每一列代表一個樣本，對於文件的前幾行：

class（必須存在），各樣本的主分組名稱；

sub_class（可預設），各樣本的次分組名稱；

subject_id（可預設），代表編號。

之後的各行，為各微生物類群的層級劃分，首先是分類水平的最高級，依次往下推，不斷的增加分類水平，不同的分類水平需要用「|」隔開。在各樣本中對應的數值，為相應微生物類群在該樣本中的總相對豐度。

對於LEfSe識別的輸入格式，可以手動使用Excel整理微生物豐度表（就是有點麻煩），也可自寫代碼完成。如下命令行展示了一個從OTU表到LEfSe輸入格式的轉換過程，代碼可見網盤附件「OTU錶轉化為LEfSe輸入格式的腳本」。

#OTU 表到 LEfSe 輸入格式的轉換過程

##測試文件說明
#otu_table.txt，OTU 豐度表，絕對豐度相對豐度都可以
#otu_table_tax_assignments.txt，OTU 注釋文件，我用 silva 注釋的（大家隨意）
#group.txt，樣本分組文件
#1-summarize_trans.py、2-lefse_trans.py，兩個自備的腳本
 
##首先藉助 qiime 工具，從 OTU 水平的豐度表獲得「界門綱目科屬」水平的微生物類群豐度表
#可使用 conda 安裝 qiime 環境
#conda create -n qiime1 qiime
 
#激活 qiime 環境
source activate qiime1
 
#在 otu_table.tsv 開頭添加一行「# Constructed from biom file」，以便將 otu_table.tsv 轉為 qiime 可識別樣式
cp otu_table.txt otu_table.tsv
sed -i '1i\# Constructed from biom file' otu_table.tsv
 
#otu_table.tsv 轉換為 biom 格式
biom convert -i otu_table.tsv -o otu_table.biom --table-type="OTU table" --to-json
 
#添加 otu 注釋信息至 biom 格式的 otu 表（otu_table.biom ）的最後一列，並將列名稱命名為 taxonomy
biom add-metadata -i otu_table.biom --observation-metadata-fp otu_table_tax_assignments.txt -o otu_table.silva.biom --sc-separated taxonomy --observation-header OTUID,taxonomy 
 
#使用 qiime 實現對 otu 表的「界門綱目科屬」統計
summarize_taxa.py -i otu_table.silva.biom -o taxonomy
 
#退出 qiime 環境
source deactivate qiime1
 
##接下來，使用自寫腳本合併「界門綱目科屬」統計的類群，並轉換為 LEfSe 的輸入格式
#合併「界門綱目科屬」統計的類群，並排序，結束後 taxonomy 路徑下獲得「otu_table.silva_all.txt」
python2 1-summarize_trans.py -i taxonomy --prefix otu_table.silva
 
#過濾去除注釋為「Other」的微生物類群
grep -v 'Other' taxonomy/otu_table.silva_all.txt > otu_table.silva_filt.txt
 
#將「otu_table.silva_filt.txt」轉化為 LEfSe 的輸入格式，並合併分組（group.txt 文件）
python2 2-lefse_trans.py otu_table.silva_filt.txt group.txt otu_table_for_lefse.txt
 
#最後的「otu_table_for_lefse.txt」就是了
#該腳本只支持添加 class 組，若存在 sub_class 組可後續手動添加即可
#就是過程有點繁瑣......主要懶得自寫代碼統計「界門綱目科屬」，就用 qiime 來做，就煩 qiime 必須識別 biom 文件，還得來迴轉格式
#後續哪天閒了重新整個不藉助 qiime 的

 

下文開始展示LefSe在線運行過程，使用的示例數據集下載自（我也放網盤附件裡了，lefse_test.txt）：

https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/lefse/input/hmp_small_aerobiosis.txt

別問我為什麼沒用上面自己轉化的數據，因為下文是先寫的，上文代碼是後寫的……

由於我沒找到關於該數據集的文字描述，所以下述幾行文字是我根據它的分組名稱猜的……

oxygen_availability，測量了樣本來源的微生物群落的氧利用能力，並據此將樣本歸為3個類別：高氧利用（High_O2）、中氧利用（Mid_O2）以及低氧利用（Low_O2）。

body_site，根據樣本來源，將樣本歸為6個類別：ear、oral、gut、nasal、vagina、skin。

subject_id，代表編號。

Bacteria等，為各類微生物類群的名稱，及其在各樣本中的相對豐度信息。

LefSe提供了在線分析方法，作為Galaxy平臺的一個分析模塊，方便了我們使用。

Galaxy平臺：http://huttenhower.sph.harvard.edu/galaxy

下文的在線運行結果，已存放在了網盤附件中，見 「LEfSe在線運行示例的結果文件」。

   

首先，把已經生成好的LEfSe配置文件上傳到雲平臺上。

 

上傳成功後，進入LEfSe分析界面，設置參數並運行LEfSe分析，以及可視化。

在線版的LEfSe由執行以下步驟的六個模塊組成，接下來我們一步步分析。

   

A模塊中的可選項主要包括指定樣本的分組類別行、子類別行等，以及是否需要標準化數據。

執行後，LEfSe將預處理輸入文件的格式，獲得LEfSe內置識別文件格式。

   

這一步是LEfSe的計算環節。

使用A模塊中提供的預處理後的文件，設置程序運行參數執行分析。計算過程概要可見上文「LEfSe原理概述」中的描述。

主要參數中，包括了Kruskall-Wallis檢驗的p值閾值、Wilcoxon檢驗的p值閾值、以及線性判別得分（LDA得分）的閾值，用作識別顯著的biomaker的標準。

該步結束後，就獲得LEfSe的計算結果了。

結果下載後可用excel打開，共5列信息，包含了每種微生物類群名稱及其對類間差異貢獻效應的大小。

第1列，每種微生物類群名稱；

第2列，每種微生物類群在各分組類別中豐度平均值中最大值的log10，如果平均豐度小於10的按照10來計算（即該列的最小值為1）；

第3列，若該微生物類群存在組間顯著差異，且線性判別得分高於設定閾值，則該列展示其所富集的分組類別名稱；若不存在顯著差異，則該列為空；

第4列，若該微生物類群存在組間顯著差異，且線性判別得分高於設定閾值，則該列展示其線性判別得分值，用以評估其對觀測到的組間差異的效應大小，該值越高代表該微生物類群越重要；若不存在顯著差異，則該列為空；

第5列，統計檢驗的p值，若顯著，則將這些微生物作為解釋類別之間差異的biomaker來看待；若不顯著，則該列值為「-」；

該步的結果可進一步為可視化模塊（C、D、E、F）提供輸入。

   

以圖形方式展示了發現的生物標誌物（模塊B的輸出）及其效應大小。

查看結果文件，以柱形圖的形式展示了顯著的微生物類群，顏色代表了其富集的分組。

數值為線性判別得分（已經過log10轉化），代表了該類群的效應大小，值越高代表了其對觀測到的組間差異貢獻越大，表明該微生物類群是更重要的biomaker。

   

以圖形方式展示了由分層要素名稱指定的分類樹中發現的biomaker（模塊B的輸出）。

查看結果文件，該結果就是常見的那種「物種進化樹」樣式。

圖中每一個節點代表一種給定的微生物類群，對於不顯著的微生物使用黃色節點標註，顯著的則賦予其它顏色，顏色代表了其富集的分組，即表明這些分支在某些分組中具有更高的豐度。

 

如果覺得該圖中，黃色（不顯著）的點過於密集影響觀測，即覺得重疊嚴重，則可以通過編輯模塊B的輸出結果，將其中的一部分不顯著的微生物類群去除。之後再將編輯後的表讀入Galaxy，作為LEfSe模塊D的輸入，如此出圖中黃色的點就少了，圖片就會「整齊」很多。

例如，我刪掉了模塊B輸出結果中的一些內容，然後將修改過後的文件重新上傳至Galaxy，同上述，通過左側界面的「Get Data」的「Upload File」上傳，上傳時注意選擇文件格式為「lefse_internal_res」。

然後在「LEfSe」中選擇「D） Plot Cladogram」，讀取新上傳後的過濾後的文件，重新作圖。示例結果如下，節點變少了，圖片「整齊」很多。（不過我好像刪的數據有點多……大家視情況來吧）

   

該步將模塊A、B的輸出作為輸入，將某特定變量（手動指定）的行值繪製為一個包含類和子類結構的豐度直方圖。

選擇展示的變量可以為已識別的biomaker或非biomaker，即可以為顯著的也可以為非顯著的，通過在繪圖選項中指定。

例如選擇了其中某微生物類群進行觀測，結果圖中以柱形圖展示了該微生物在各樣本中的豐度信息，顏色代表了樣本的分組狀況。

   

與模塊E相比，該步可一次將所有變量的行值統一繪製為包含類和子類結構的豐度直方圖，結果以zip壓縮包形式保存。

可以只輸出已識別為biomaker的變量，也可以輸出所有變量，通過在繪圖選項中指定。

例如上述只選擇了識別為biomaker（即顯著）的變量進行展示，可在結果中下載打包好的壓縮包。下圖是其中之一的展示圖。

然後是本地版LEfSe工具的使用，可選項更多，更加靈活。

LEfSe目前也已打包在conda中，方便了我們安裝使用。如下繼續使用上述示例數據，展示本地版LEfSe工具的使用，shell命令行。

本地運行的輸出結果示例，也放在了網盤附件中，見 「LEfSe本地運行示例的結果文件」。由於我在本地運行和在線版運行時的設定參數不同，所以兩個示例輸出是不一致的，所以不要疑問為啥二者不一致。

各結果文件的說明，參考上文在線版的輸出示例即可。

#conda 安裝 LEfSe
conda install -c biobakery lefse
#或者
conda install -c bioconda lefse
 
##對應於上述 A-F 6 個模塊，本地版的命令行操作示例如下
 
#A，設置 LEfSe 的數據格式，詳情 format_input.py -h
#-c，指定 class 的行（必須指定）；-s，指定 sub_class 的行（可預設）；-u，指定 subject_id 的行（可預設）；-o，設置歸一化值，默認 -1 即不執行標準化
#註：版本問題，有時 format_input.py 命令找不到，可能為 lefse-format_input.py
format_input.py lefse_test.txt A_lefse_test.in -c 1 -s 2 -u 3 -o 1000000
 
#B，LEfSe 分析，詳情 run_lefse.py -h
#-l 2.0，設定 LDA 得分的對數值的最低閾值為 2
run_lefse.py A_lefse_test.in B_lefse_test.res -l 2.0
 
#備註：對於 B 的輸出，我們可以選擇從中刪一些不必要（不顯著）的數據，以增強 C、D、E 作圖時的美感
 
#C，繪製 LEfSe 得分值，詳情 plot_res.py -h
#註：版本問題，有時 plot_res.py 命令找不到，可能為 lefse-plot_res.py
plot_res.py B_lefse_test.res C_lefse_test.lda.pdf --format pdf --dpi 150 --width 16
 
#D，繪製進化分支圖，詳情 plot_cladogram.py -h
#註：版本問題，有時 plot_cladogram.py 命令找不到，可能為 lefse-plot_cladogram.py
plot_cladogram.py B_lefse_test.res D_lefse_test.cladogram.pdf --format pdf --dpi 150
 
#E，單張圖的展示略，直接使用 F 繪製所有的圖
 
#F，繪製差異特徵，詳情 plot_features.py -h
#註：版本問題，有時 plot_features.py 命令找不到，可能為 lefse-plot_features.py
mkdir -p F_out
plot_features.py A_lefse_test.in B_lefse_test.res F_out/lefse_test --format pdf --dpi 200

     

https://bitbucket.org/biobakery/biobakery/wiki/lefse

Segata N, Izard J, Waldron L, et al. Metagenomic biomarker discovery and explanation. Genome biology, 2011, 12(6).