近年來(lái)發(fā)展起來(lái)的宏基因組學(xué)技術(shù)避開(kāi)傳統(tǒng)的微生物分離培養(yǎng)方法直接從環(huán)境樣品中提取總DNA, 通過(guò)構(gòu)建和篩選宏基因組文庫(kù)來(lái)獲得新的功能基因和生物活性物質(zhì), 宏基因組文庫(kù)既包括了可培養(yǎng)的, 又包括了未可培養(yǎng)的微生物遺傳信息, 因此增加了獲得新生物活性物質(zhì)的機(jī)會(huì)。特定生物種基因組研究使人們的認(rèn)識(shí)單元實(shí)現(xiàn)了從單一基因到基因集合的轉(zhuǎn)變，宏基因組研究將使人們擺脫物種界限，揭示更高更復(fù)雜層次上的生命運(yùn)動(dòng)規(guī)律。

在目前的基因結(jié)構(gòu)功能認(rèn)識(shí)和基因操作技術(shù)背景下，細(xì)菌宏基因組成為研究和開(kāi)發(fā)的主要對(duì)象。細(xì)菌宏基因組、細(xì)菌人工染色體文庫(kù)篩選和基因系統(tǒng)學(xué)分析使研究者能更有效地開(kāi)發(fā)細(xì)菌基因資源，更深入地洞察細(xì)菌多樣性。但是宏基因組研究的實(shí)驗(yàn)數(shù)據(jù)中往往含有成千上萬(wàn)種微生物，如何準(zhǔn)確的篩選出不同實(shí)驗(yàn)條件或環(huán)境條件下的關(guān)鍵生物標(biāo)志物，指導(dǎo)后續(xù)的實(shí)驗(yàn)和功能挖掘成為了一個(gè)重要的研究方向。

本代碼通過(guò)基于宏基因組的微生物豐度數(shù)據(jù)（物種豐度或OUT豐度），結(jié)合原始生物樣本的信息，利用非參數(shù)檢驗(yàn)及LDA（Linear Discriminant Analysis，線性判別分析）的算法，篩選出在不同實(shí)驗(yàn)條件下顯著差異的微生物，并估算出其相對(duì)差異的分?jǐn)?shù)（LDA score），以此來(lái)代表該物種在某種實(shí)驗(yàn)環(huán)境下的相對(duì)豐度。

通過(guò)結(jié)合原始樣本信息，既可以得出在不同條件下的顯著的biomarker（生物標(biāo)志物），幫助后續(xù)的數(shù)據(jù)分析。用戶不需要設(shè)置復(fù)雜的模型參數(shù)，只需要指定微生物的豐度矩陣以及對(duì)應(yīng)的樣本信息即可，樣本信息文件中可包含多種實(shí)驗(yàn)分組條件。代碼會(huì)匹配根據(jù)用戶指定的分組計(jì)算出顯著的生物標(biāo)志物，最多同時(shí)可以指定兩個(gè)分組條件一起計(jì)算。計(jì)算完成之后，會(huì)輸出biomarker的分?jǐn)?shù)矩陣，同時(shí)對(duì)分?jǐn)?shù)矩陣使用直方圖進(jìn)行可視化，展示biomarker。

使用方法：

Rscript run_LDA.R [-h] [-v] [-q] [-k float number] [-w float number]

[-l float number] [-g string] [-b string] [-o string]

Abundance File Meta File

參數(shù)說(shuō)明：

usage:?

run_LDA.R [-h] [-v] [-q] [-k float number] [-w float number]

[-l float number] [-g string] [-b string] [-o string]

Abundance File Meta File

positional arguments:Abundance File ???????Species/OTU Abundance File to be used

Meta File ????????????Samples Metadata file

optional arguments:-h, --help ???????????show this help message and exit-v, --version ????????Print software Version [default]-k float number, --kruskal.threshold float numberThe p-value for the Kruskal-Wallis Rank Sum Test(default 0.05).-w float number, --wilcox.threshold float numberThe p-value for the Wilcoxon Rank-Sum when 'blockCol'is present (default 0.05).-l float number, --lda.threshold float numberThe effect size threshold (default 2.0).-g string, --group.var stringColumn name in meta data file indicating groups,usually a factor with two levels (default "GROUP")-b string, --block.var stringOptional column name in meta data file indicatingblocks, usually a factor with two levels (default"GROUP")-o string, --output.prefix string

prefix for output file(default "lda_result")

操作步驟：

1、打開(kāi)命令行界面，輸入“Rscript run_LDA.R -h”查看軟件詳細(xì)的幫助文檔，確定該程序所需的輸入文件和關(guān)鍵參數(shù)的設(shè)置。

2、用戶根據(jù)幫助文檔中的參數(shù)說(shuō)明內(nèi)容，對(duì)參數(shù)進(jìn)行設(shè)置。這里，必須輸入?yún)?shù)有2個(gè)，分別是Abundance File，表示物種的豐度矩陣文件，以物種或OTU為行，樣本為列，保存為csv文件；Meta File 表示各個(gè)樣本的meta信息文件，第一列為樣本名，后面每一列為樣本的某一種分組或?qū)嶒?yàn)類型，也保存為csv文件。需要確保Abundance File中出現(xiàn)的樣本名，在Meta File 中必須能夠找到。-g，表示選取Meta矩陣中的某一列作為分組依據(jù)，進(jìn)行LDA分析。其余參數(shù)為可選參數(shù)，在用戶沒(méi)有顯式指定的情況下，將采用默認(rèn)值運(yùn)行。具體說(shuō)明如下：-h，輸出軟件幫助文檔；-v，輸出軟件版本號(hào)；-k，kruskal.threshold，kruskal 非參數(shù)檢驗(yàn)的p值，只有小于此閾值的微生物會(huì)被保留下來(lái)；-w, wilcox.threshold, wilcox非參數(shù)檢驗(yàn)的的p值，只有小于此閾值的微生物會(huì)被保留下來(lái)；-l, lda.threshold, lda score的回歸系數(shù)閾值，絕對(duì)值大于此閾值的會(huì)被保留；-b，第二個(gè)分組變量的列名；-o，結(jié)果輸出文件的前綴。參數(shù)符號(hào)的string，number提示該參數(shù)的數(shù)據(jù)類型。

3、完成參數(shù)提交后，按下確認(rèn)鍵，整個(gè)程序即正式開(kāi)始進(jìn)入執(zhí)行。每步執(zhí)行內(nèi)容都會(huì)給出提示。程序執(zhí)行完畢后，界面會(huì)顯示” Finished LDA "結(jié)束語(yǔ)。

結(jié)果展示：

共輸出一個(gè)表格csv格式文件，一個(gè)pdf格式的圖片文件

1、lda_result_plot.pdf

該圖表示在用戶指定的分組中顯著差異的大于設(shè)定閾值的biomarker。LDA score分值越大，差異越大，表示在該組中豐度顯著高于其他組。不同的顏色用于區(qū)分不同的分組。

2、lda_result.csv

該表格表示各個(gè)差異物種的具體LDA分?jǐn)?shù)值，共有三列，第一列是具體的物種或OUT名稱，第二列是具體的LDA Score，是對(duì)原始LDA score去log10之后的數(shù)值, 第三列是該物種顯著富集的分組。

特別說(shuō)明：本代碼經(jīng)申請(qǐng)軟件著作權(quán)，僅轉(zhuǎn)讓使用權(quán)，不轉(zhuǎn)讓所有權(quán)

如需代碼及示例數(shù)據(jù)等文件，請(qǐng)掃碼聊天框回復(fù) “代碼”領(lǐng)?。?/span>

寫在文末：

如果您近期想做生信方面的文章而苦于沒(méi)有思路，或者不知道如何來(lái)入手生信分析，或者兌具體的某一個(gè)圖有作圖需求，都可以掃碼咨詢小云，我們有專業(yè)的技術(shù)團(tuán)隊(duì)，生信熱點(diǎn)思路設(shè)計(jì)、生信分析、熱點(diǎn)方向生信挖掘等，如有需要，可掃碼下方二維碼了解詳情