LEfSe概述
LEfSe透過將具有統計學意義的標準檢驗與表徵生物學一致性和效應相關性的附加檢驗結合起來，確定最有可能解釋類別之間差異的特徵（如物種、進化枝、基因或功能，統稱為biomaker）。類別比較方法通常用於預測biomaker，這些標記由違反類別之間無差異零假設的特徵組成；並透過檢測特徵子集的丰度模式，估計顯著變化的程度；效應大小提供了每個特徵對類別差異貢獻大小的估計。
以下是對LEfSe工作原理的簡單概括，細節部分可參考Segata等（2011）的論文。
Data：輸入資料包含m個樣本（列）的集合，每個樣本由n個變數（行，通常進行行歸一化，紅色代表高值，綠色代表低值）組成；每個樣本都對應了一些類別（即樣本分組），代表正在研究的主要生物學假設，可以為兩組或多組，也可以包含反映類內分組的子類別標籤。
Step3：所得的向量子集用於構建線性判別分析（LDA）模型，利用類間的相對差異對變數進行降維和分類，輸出包含一列與類別有關的特徵列表，這些特徵與類別中的子類別分組保持一致，並根據它們區分類別的效應大小進行排名。
因此LEfSe是用於對不同生物學方面的相關性進行排名以及進行進一步研究和分析的有價值的工具。它強調統計意義和生物相關性，能夠在組與組之間尋找具有統計學差異的biomaker。在該方法中引入先驗知識作為監督，解決了傳統上與高維資料探勘有關的挑戰。
接下來，展示LEfSe工具的使用，以及對結果部分的說明。
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線上執行
LefSe提供了線上分析方法，作為Galaxy平臺的一個分析模組，方便了我們使用。
Galaxy平臺：http://huttenhower.sph.harvard.edu/galaxy
下文的線上執行結果，已存放在了網盤附件中，見“LEfSe線上執行示例的結果檔案”。
上傳待分析的資料
首先，把已經生成好的LEfSe配置檔案上傳到雲平臺上。
上傳成功後，進入LEfSe分析介面，設定引數並執行LEfSe分析，以及視覺化。
線上版的LEfSe由執行以下步驟的六個模組組成，接下來我們一步步分析。
（A）設定LEfSe的資料格式
A模組中的可選項主要包括指定樣本的分組類別行、子類別行等，以及是否需要標準化資料。
執行後，LEfSe將預處理輸入檔案的格式，獲得LEfSe內建識別檔案格式。
（B）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）提供輸入。
（C）繪製LEfSe得分值
以圖形方式展示了發現的生物標誌物（模組B的輸出）及其效應大小。
檢視結果檔案，以柱形圖的形式展示了顯著的微生物類群，顏色代表了其富集的分組。
數值為線性判別得分（已經過log10轉化），代表了該類群的效應大小，值越高代表了其對觀測到的組間差異貢獻越大，表明該微生物類群是更重要的biomaker。
（D）繪製進化分支圖
以圖形方式展示了由分層要素名稱指定的分類樹中發現的biomaker（模組B的輸出）。
檢視結果檔案，該結果就是常見的那種“物種進化樹”樣式。
圖中每一個節點代表一種給定的微生物類群，對於不顯著的微生物使用黃色節點標註，顯著的則賦予其它顏色，顏色代表了其富集的分組，即表明這些分支在某些分組中具有更高的丰度。
如果覺得該圖中，黃色（不顯著）的點過於密集影響觀測，即覺得重疊嚴重，則可以透過編輯模組B的輸出結果，將其中的一部分不顯著的微生物類群去除。之後再將編輯後的表讀入Galaxy，作為LEfSe模組D的輸入，如此出圖中黃色的點就少了，圖片就會“整齊”很多。
例如，我刪掉了模組B輸出結果中的一些內容，然後將修改過後的檔案重新上傳至Galaxy，同上述，透過左側介面的“Get Data”的“Upload File”上傳，上傳時注意選擇檔案格式為“lefse_internal_res”。
然後在“LEfSe”中選擇“D） Plot Cladogram”，讀取新上傳後的過濾後的檔案，重新作圖。示例結果如下，節點變少了，圖片“整齊”很多。（不過我好像刪的資料有點多……大家視情況來吧）
（E）繪製單個變數的分佈
該步將模組A、B的輸出作為輸入，將某特定變數（手動指定）的行值繪製為一個包含類和子類結構的丰度直方圖。
選擇展示的變數可以為已識別的biomaker或非biomaker，即可以為顯著的也可以為非顯著的，透過在繪圖選項中指定。
例如選擇了其中某微生物類群進行觀測，結果圖中以柱形圖展示了該微生物在各樣本中的丰度資訊，顏色代表了樣本的分組狀況。
（F）繪製差異特徵
與模組E相比，該步可一次將所有變數的行值統一繪製為包含類和子類結構的丰度直方圖，結果以zip壓縮包形式儲存。
可以只輸出已識別為biomaker的變數，也可以輸出所有變數，透過在繪圖選項中指定。
例如上述只選擇了識別為biomaker（即顯著）的變數進行展示，可在結果中下載打包好的壓縮包。下圖是其中之一的展示圖。
LEfSe本地執行
然後是本地版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).
