關於巨集基因組常用的有參分析流程，主要是快速獲得物種組成和功能組成，之前分享了

今天再介紹來自同一作者的另一個軟體，可以一步完成功能和代謝組成。

HUMAnN2: The HMP Unified Metabolic Analysis Network 2，它在巨集基因組研究中非常有用，通過這個分析，不僅能獲得微生物的物種丰度資訊，還能準確有效地獲得微生物代謝途徑和功能模組資訊。

中文版本翻譯日期(2018-05-01)

HUMAnN是基於巨集基因組、巨集轉錄組資料分析微生物通路丰度的有效工具。這一過程稱為功能譜，目的是描述群體成員的代謝潛能。可以回答微生物群體成員可能幹什麼，或在幹什麼的問題

軟體特點：

1. 可對已知和末知生物分析群體功能譜
   
   • MetaPhlAn2和ChocoPhlAn泛基因組資料庫， 可以更快速準確獲得功能譜
   • 物種包括古菌、細菌、真核生物和病毒
2. 可獲得基因組、基因和通路層面的結果
   
   • UniRef資料庫提供基因家族的定義
   • MetaCyc通路基因通路的定義
   • MinPath提供定義的最小通路集
3. 簡單的使用介面(單行命令工作流)
   
   • 使用者只需提供質控的巨集基因組或巨集轉錄組資料
4. 加速序列比對
   
   • 採用Bowtie2加速核酸水平搜尋
   • 採用Diamond加速翻譯蛋白水平搜尋

HUMAnN2工作流程圖

安裝

如果你安裝過python，且有pip安裝工具，可以輕鬆安裝humann2

# 軟體安裝
pip install humann2

# 或可選手動下載安裝
wget //files.pythonhosted.org/packages/43/07/ec41577c3c1f9b578875ade8ed549d14fc2944c13cb7504579d542b62a69/humann2-0.11.1.tar.gz    

# 前面仍不成功，推薦conda安裝更快更好用
conda install humann2

# 測試安裝
humann2_test

# 比如我使用conda安裝程式至/conda/bin目錄，且沒有新增環境變數，可以使用絕對路徑呼叫程式

# 下載資料庫
wd=/conda/bin
$wd/humann2_databases --available
# 5.37GB
$wd/humann2_databases --download chocophlan full /data/humann2
# 5.87GB，解壓後11G
$wd/humann2_databases --download uniref uniref90_diamond /data/humann2
 

依賴關係

# Diaomond http://ab.inf.uni-tuebingen.de/software/diamond/
wget http://github.com/bbuchfink/diamond/releases/download/v0.9.21/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz
sudo ln -fs `pwd`/diamond /usr/local/bin/

## 分析實戰

輸入檔案為fastq，輸出檔案為指定目錄中有各定量表格

cd ~/ath/jt.terpene.meta/clean_data/JT-545
# 可接受壓縮檔案fastq，並自建目錄
$wd/humann2 --input 25/JT-545_25.rmhost.1.fq.gz --output humann2_25 &
$wd/humann2 --input 26/JT-545_26.rmhost.1.fq.gz --output humann2_26 &
$wd/humann2 --input 27/JT-545_27.rmhost.1.fq.gz --output humann2_27 &

輸出檔案

輸出檔案位於輸入目錄中的輸出目錄

1. 基因家族檔案

# Gene Family   $SAMPLENAME_Abundance-RPKs
UNMAPPED        187.0
UniRef50_unknown        150.0
UniRef50_unknown|g__Bacteroides.s__Bacteroides_fragilis 150.0
UniRef50_A6L0N6: Conserved protein found in conjugate transposon    67.0
UniRef50_A6L0N6: Conserved protein found in conjugate transposon|g__Bacteroides.s__Bacteroides_fragilis 57.0
UniRef50_A6L0N6: Conserved protein found in conjugate transposon|g__Bacteroides.s__Bacteroides_finegoldii   5.0
UniRef50_A6L0N6: Conserved protein found in conjugate transposon|g__Bacteroides.s__Bacteroides_stercoris    4.0
UniRef50_A6L0N6: Conserved protein found in conjugate transposon|unclassified   1.0
UniRef50_O83668: Fructose-bisphosphate aldolase 60.0
UniRef50_O83668: Fructose-bisphosphate aldolase|g__Bacteroides.s__Bacteroides_vulgatus  31.0
UniRef50_O83668: Fructose-bisphosphate aldolase|g__Bacteroides.s__Bacteroides_thetaiotaomicron  22.0
UniRef50_O83668: Fructose-bisphosphate aldolase|g__Bacteroides.s__Bacteroides_stercoris 7.0

• 檔名：OUTPUTDIR/$OUTPU{T}_{D}IR/$SAMPLENAME_genefamilies.tsv
• 群體中每個基因家族的丰度。基因家族是一組進化上相關的編碼蛋白序列，通常具有相似功能。
• 基因家族的丰度在群體水平分級顯著，顯示已知和未知物種的貢獻度。
• 使用MetaPhlAn2軟體和ChocoPhlAn資料庫，檢索核酸翻譯的蛋白資料庫
• 基因家族的丰度採用RPK(每Kb的reads)以標準化不同的基因長度；RPK單位代表基因或轉錄本在群體中的拷貝數；RPK值可進一步求和標準化，用於不同樣品測序深度的比較。
• 如果輸入檔案是基因表，不會建立基因家族檔案
• UNMAPPED是兩步核酸和蛋白搜尋後，仍無法比對的reads數量。

• UniRef50_unknown 代表可比對ChocoPhlAn，但沒有註釋

  1. 通路丰度檔案

#Pathway   $SAMPLENAME_Abundance
UNMAPPED    140.0
UNINTEGRATED    87.0
UNINTEGRATED|g__Bacteroides.s__Bacteroides_caccae   23.0
UNINTEGRATED|g__Bacteroides.s__Bacteroides_finegoldii   20.0
UNINTEGRATED|unclassified   12.0
PWY0-1301: melibiose degradation    57.5
PWY0-1301: melibiose degradation|g__Bacteroides.s__Bacteroides_caccae   32.5
PWY0-1301: melibiose degradation|g__Bacteroides.s__Bacteroides_finegoldii   4.5
PWY0-1301: melibiose degradation|unclassified   3.0
PWY-5484: glycolysis II (from fructose-6P)  54.7
PWY-5484: glycolysis II (from fructose-6P)|g__Bacteroides.s__Bacteroides_caccae 16.7
PWY-5484: glycolysis II (from fructose-6P)|g__Bacteroides.s__Bacteroides_fi

• 檔名：OUTPUT_DIR/$SAMPLENAME_pathabundance.tsv
• 代表群體中通路的丰度
• 即有群體水平，又有物種水平丰度
• 通路按丰度大小排序，物種組分也按丰度大小排序，全為0的通路不輸出
• 通路的比例是是完整拷貝的丰度，如線性通路Gene1-4，分別為10，5，5，5。則按5計算。
• 與基因不同，通路的丰度並一定是群體組分的總合。A物種[5, 5, 10, 10]為5個拷貝，B物種[10, 10, 5, 5]為5個拷貝，而總體有15個拷貝； 詳細的計算說明見英文幫助原文
• 對於單個基因必須的反應步驟為零丰度時，可進行所需最低丰度填充。
• MetaCyc預設定義最簡通路解析群體觀測的代謝通路；
• 使用者也可以自定義通路資料庫

• 非線性基因拷貝數、無法比對序列處理方法請參考英文原文

  1. 通路覆蓋度檔案

# Pathway   $SAMPLENAME_Coverage
UNMAPPED    1.0
UNINTEGRATED    1.0
UNINTEGRATED|g__Bacteroides.s__Bacteroides_caccae   1.0
UNINTEGRATED|g__Bacteroides.s__Bacteroides_finegoldii   1.0
UNINTEGRATED|unclassified   1.0
PWY0-1301: melibiose degradation    1.0
PWY0-1301: melibiose degradation|g__Bacteroides.s__Bacteroides_caccae   1.0
PWY0-1301: melibiose degradation|g__Bacteroides.s__Bacteroides_finegoldii   1.0
PWY0-1301: melibiose degradation|unclassified   1.0
PWY-5484: glycolysis II (from fructose-6P)  1.0
PWY-5484: glycolysis II (from fructose-6P)|g__Bacteroides.s__Bacteroides_caccae 0.7
PWY-5484: glycolysis II (from fructose-6P)|g__Bacteroides.s__Bacteroides_finegoldii 0.7
PWY-5484: glycolysis II (from fructose-6P)|unclassified 0.3

• 檔名：$OUTPUT_DIR/$SAMPLENAME_pathabundance.tsv
• 檔案提供了一種有(1)和無(0)的群體通路計演算法，而不是相對丰度
• 每個反應有置信得分
• 計算通路覆蓋與相對丰度一致
• 通路丰度在群體水平和物種水平計算
• 群體水平比物種水平更可信
• 只輸出非零丰度的通路

• 通路覆蓋度與通路丰度順序相同

  1. 中間臨時檔案

• Bowtie2比對結果($DIR/$SAMPLENAME_bowtie2_aligned.sam)

• 簡化Bowtie2比對結果$DIR/$SAMPLENAME_bowtie2_aligned.tsv
• Bowtie2資料庫索引$DIR/$SAMPLENAME_bowtie2_index*
• Bowtie2末比對的reads$DIR/$SAMPLENAME_bowtie2_unaligned.fa
• 自定義ChocoPhlAn資料庫$DIR/$SAMPLENAME_custom_chocophlan_database.ffn
• MetaPhlAn2的bowtie2比對結果$DIR/$SAMPLENAME_metaphlan_bowtie2.txt
• MetaPhlAn2 bugs list$DIR/$SAMPLENAME_metaphlan_bugs_list.tsv
• 翻譯後的比對結果$DIR/$SAMPLENAME_$TRANSLATEDALIGN_aligned.tsv
• 翻譯後仍末比對序列$DIR/$SAMPLENAME_$TRANSLATEDALIGN_unaligned.fa
• 日誌檔案$DIR/$SAMPLENAME.log

配置軟體

# 顯示引數
$wd/humann2_config --print
# 修改引數格式
$wd/humann2_config --update $SECTION $NAME $VALUE
# 如修改執行緒數
$wd/humann2_config --update run_modes threads 12

HUMAnN2小工具

humann2_barplot

Basic usage: $ humann2_barplot --input $TABLE.tsv --feature $FEATURE --outfile $FIGURE
$TABLE.tsv = a stratified HUMAnN2 output file
$FEATURE = Feature from the table to plot (defaults to first feature)
$FIGURE = Where to save the figure
Run with -h to see additional command line options

可選擇某個Feature進行柱狀圖視覺化。–help引數可檢視相關排序、標準化選項。

合併多樣品結果為表格

此步非常重要，我們無法多少個樣品，humann2結果僅為一列。多樣品需經本步合併為矩陣，方便下游統計分析和差異比較。

Basic usage: $ humann2_join_tables --input $INPUT_DIR --output $TABLE
$INPUT_DIR = a directory containing gene/pathway tables (tsv or biom format)
$TABLE = the file to write the new single gene table (biom format if input is biom format)
Optional: --file_name $STR will only join gene tables with $STR in file name
Run with -h to see additional command line options

其它小工具

• 構建自定義資料庫 humann2_build_custom_database
• 檢視屬水平基因家族與通路 humann2_gene_families_genus_level
• 增加物種分類(降低可信度) humann2_infer_taxonomy
• 合併MetaPhlAn2分類結果 humann2_reduce_table
• 對樣品、通路進行合併/重分組操作humann2_regroup_table
• 對Feature進行重新命名 humann2_rename_table
• 表格標準化 humann2_renorm_table
• 巨集轉錄組標準化 humann2_rna_dna_norm
• 將結果分為分層、末分層兩個檔案 humann2_split_stratified_table
• 將多樣品表分類單樣品 humann2_split_table
• 挖掘菌株水平差異 humann2_strain_profiler
• 輸出組成通路的基因丰度 humann2_unpack_pathways

其它說明

1. 選擇基因家族精度的水平
2. 選擇UniRfe90還是50? 推薦90
3. 選擇翻譯搜尋模式
   
   1. Bypass translated search: 無法比對泛基因組的儲存，再比對蛋白，結果中沒有末分類的層級
   2. Filtered translated search，是EC-filtered protein database(是UniRef中Level4層面國際生物化學聯合會酶委員分類)的預設方案
   3. Comprehensive translated search 使用最綜合的蛋白資料庫，但速度比2慢5倍。預設為方案2
4. 雙端序列：HUMAnN2不考慮雙端序列，全當作單端
5. PICRUSt輸出可繼續用本軟體分析，需要下載KEGG完整資料庫，拆分預測巨集基因組為每個樣品單個檔案，再執行humann2，再合併。詳見官網
6. 使用KEEG資料庫，目前新版收費，免費版本最新為v56，可同HUMAnN1一起下載
7. 結果可使用QIIME的core_diversity_analyses.py進行多樣性分析
8. HUMAnN2分析巨集轉錄組

常見問題

1. 如何輸出分析過程更多資訊？新增--verbose引數
2. 如何使用多核？--threads $CORES或修改預設設定
3. 如何清空臨時檔案？--remove-temp-output
4. 指定ChocoPhlAn資料庫位置？--nucleotide-database $DIR
5. 指定UniRef資料庫位置？--protein-database $DIR
6. 使用Metaphlan2結果繼續分析？--taxonomic-profile bugs_list.tsv
7. 修改樣本名？--output-basename $NAME
8. 去除分層的結果？--remove-stratified-output
9. 使用unipathways databases？--pathways unipathway
10. 輸出biom格式結果--output-format biom
11. 修改相似度閾值--identity-threshold <50.0>
12. 修改metaphlan2引數--metaphlan-options="-t rel_ab"

報錯與解決

1. CRITICAL ERROR: Can not call software version for diamond

diamond沒有在環境變數，下載解壓並確保新增到環境變數

1. The database file for MetaPhlAn does not exist at /mnt/bai/yongxin/software/metaphlan2/db_v20/mpa_v20_m200.pkl . Please provide the location with –metaphlan-options

沒有找到metaphlan2的資料庫，是metaphlan2新版本目錄更改了位置，永久方法是建一箇舊位置的硬鏈。

進入metaphlan2安裝目錄

mkdir db_v20
ln `pwd`/databases/mpa_v20_m200.* db_v20/

1. CRITICAL ERROR: Error executing: /mnt/bai/public/bin/diamond blastx –query /mnt/bai/yongxin/ath/jt.terpene.meta/clean_data/JT-545/humann2/JT-545_25.rmhost.1_humann2_temp/JT-545_25.rmhost.1_bowtie2_unaligned.fa –evalue 1.0 –threads 1 –max-target-seqs 20 –outfmt 6 –db /data/humann2/uniref/uniref90_annotated.1.1 –out /mnt/bai/yongxin/ath/jt.terpene.meta/clean_data/JT-545/humann2/JT-545_25.rmhost.1_humann2_temp/tmp2_lUDg/diamond_m8_TL5Rl_ –tmpdir /mnt/bai/yongxin/ath/jt.terpene.meta/clean_data/JT-545/humann2/JT-545_25.rmhost.1_humann2_temp/tmp2_lUDg

/mnt/bai/public/bin/diamond是個目錄，不知為什麼系統會找這個目錄當前程式，系統我也裝在了 /usr/local/bin/diamond 中。修改此目錄為程式

巨集基因組相關學習資源

如果基礎知識體系不完善，自學存在困難的小夥伴，急時上車也是不錯的選擇。成為實驗中不可或缺的人，趕快點選閱讀原文報名我們的培訓，加速你入行！http://www.ehbio.com/Training

猜你喜歡

寫在後面

為鼓勵讀者交流、快速解決科研困難，我們建立了“巨集基因組”專業討論群，目前己有國內外1500+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註“姓名-單位-研究方向-職稱/年級”。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。

學習擴增子、巨集基因組科研思路和分析實戰，關注“巨集基因組”