宏基因組序列物種分類之kraken 1/2和Bracken的使用

2021-02-21 宏基因組

細菌基因組測序完，想看看樣本有沒有被其他的菌汙染?人的轉錄組測序完，想快速看看人、微生物的序列的比例?元/宏基因組測序完，想快速獲得樣本中物種的豐度信息？

REFERENCE

Wood DE, Salzberg SL: Kraken: ultrafast metagenomic sequence classification using exact alignments.Genome Biology 2014, 15:R46.

Kraken 1

Kraken 1是2014年Wood DE在Genome Biology發表的宏基因組序列分類軟體，能夠快速對宏基因樣品中的reads進行分類。

Kraken在序列比對環節基於精確k-mer匹配和精簡資料庫的方法，採取精確匹配，其核心是Kraken有一種特殊資料庫，用以預先計算序列中包含的特殊的Kmer序列。

下面是來自kraken官網關於各分類器的測評結果：

Kraken速度很快，精度較低，適用於做微生物檢測的預處理。通過一些實際的數據測試發現：與Metaphlan2相比，Kraken速度較快，獲得的物種數目較多，相對應的假陽性率也較高。

Usage: kraken [options] <filename(s)>

Options:

--db NAME Name for Kraken DB

(default: none)

--threads NUM Number of threads (default: 1)

--fasta-input Input is FASTA format

--fastq-input Input is FASTQ format

--fastq-output Output in FASTQ format

--gzip-compressed Input is gzip compressed

--bzip2-compressed Input is bzip2 compressed

--quick Quick operation (use first hit or hits)

--min-hits NUM In quick op., number of hits req'd for classification

NOTE: this is ignored if --quick is not specified

--unclassified-out FILENAME

Print unclassified sequences to filename

--classified-out FILENAME

Print classified sequences to filename

--out-fmt FORMAT Format for [un]classified sequence output. supported

options are: {legacy, paired, interleaved}

--output FILENAME Print output to filename (default: stdout); "-" will

suppress normal output

--only-classified-output

Print no Kraken output for unclassified sequences

--preload Loads DB into memory before classification

--paired The two filenames provided are paired-end reads

--check-names Ensure each pair of reads have names that agree

with each other; ignored if --paired is not specified

--help Print this message

--version Print version information

If none of the *-input or *-compressed flags are specified, and the

file is a regular file, automatic format detection is attempted.

$ kraken --threads 40 --db minikraken_20171013_4GB --preload --

paired --fastq-input --gzip-compressed ${B}_1.fastq.gz ${B}_2.fastq.gz | kraken-report --db minikraken_20171013_4GB > "$B"_kraken.tab

$ less -S "$B"_kraken.tab

50.30 27375076 27375076 U 0 unclassified

49.70 27047187 13633 - 1 root

49.35 26857009 348 - 131567 cellular organisms

49.35 26855582 105532 D 2 Bacteria

33.57 18270135 0 - 1783270 FCB group

33.57 18269977 5 - 68336 Bacteroidetes/Chlorobi group

33.57 18269761 107058 P 976 Bacteroidetes

33.34 18144208 69 C 200643 Bacteroidia

33.34 18144105 961382 O 171549 Bacteroidales

30.47 16580771 0 F 815 Bacteroidaceae

30.47 16580771 2380340 G 816 Bacteroides

22.89 12454849 9623614 S 821 Bacteroides vulgatus

5.20 2831235 2831235 - 435590 Bacteroides vulgatus ATCC 8482

0.73 396849 0 S 357276 Bacteroides dorei

0.73 396849 396849 - 997877 Bacteroides dorei CL03T12C01

0.60 326487 326471 S 28116 Bacteroides ovatus

0.00 16 16 - 1379690 Bacteroides ovatus V975

0.42 227354 141795 S 818 Bacteroides thetaiotaomicron

0.16 85559 85559 - 226186 Bacteroides thetaiotaomicron VPI-5482

0.38 206113 206113 S 1796613 Bacteroides caecimuris

0.33 180286 180286 S 47678 Bacteroides caccae

0.29 156976 68554 S 817 Bacteroides fragilis

0.10 55847 55847 - 862962 Bacteroides fragilis 638R

0.03 17859 17859 - 295405 Bacteroides fragilis YCH46

0.03 14716 14716 - 272559 Bacteroides fragilis NCTC 9343

如上是Kraken的結果，可以看出它沒有估算出物種的豐度。 

Bracken

這時可以使用另一款軟體Bracken (Bayesian Reestimation of Abundance with KrakEN)，它是一種從宏基因組數據中高度準確的計算物種豐度的統計方法。

$ bracken -h

Usage: bracken -d MY_DB -i INPUT -o OUTPUT -r READ_LEN -l LEVEL -t THRESHOLD

MY_DB location of Kraken database

INPUT Kraken REPORT file to use for abundance estimation

OUTPUT file name for Bracken default output

READ_LEN read length to get all classifications for (default: 100)

LEVEL level to estimate abundance at [options: D,P,C,O,F,G,S] (default: S)

THRESHOLD number of reads required PRIOR to abundance estimation to perform reestimation (default: 0)

$ bracken -d minikraken_20171013_4GB/ -i $B\_kraken.tab -t 10 -o $B.out

#獲得如下結果：

name taxonomy_id taxonomy_lvl kraken_assigned_reads added_reads new_est_reads fraction_total_reads

Uncultured phage WW-nAnB strain 2 1449896 S 388 113 501 0.00007

Uncultured phage WW-nAnB strain 3 1449897 S 349 131 480 0.00007

Aureimonas sp. AU20 1349819 S 108 39 147 0.00002

Phaeobacter piscinae 1580596 S 21 4 25 0.00000

Sinorhizobium sp. CCBAU 05631 794846 S 25 12 37 0.00001

Mucilaginibacter sp. BJC16-A31 1234841 S 31 0 31 0.00000

Arcanobacterium phocae 131112 S 34 0 34 0.00000

Kineococcus radiotolerans 131568 S 82 3 85 0.00001

Actinomyces radingae 131110 S 314 5 319 0.00005

Sediminicola sp. YIK13 1453352 S 19 0 19 0.00000

Methylotenera mobilis 359408 S 42 1 43 0.00001

Stenotrophomonas rhizophila 216778 S 76 52 128 0.00002

Acholeplasma oculi 35623 S 12 0 12 0.00000

Dictyoglomus turgidum 513050 S 12 0 12 0.00000

Chelatococcus sp.

NOTE: Kraken 2 is the newest version of Kraken (See Kraken 2's Webpage for details). Kraken 1 will continue to be available via the Kraken 1 Github page, but it is no longer being supported.

Kraken 2

與Kraken 1相比，Kraken 2有了很大的改進：

更快速的構建資料庫

資料庫的佔用存儲空間更少

更快的分類速度

還能支持 16S Databases包括Greengenes, SILVA,和 RDP。

$ kraken2

Need to specify input filenames!

Usage: kraken2 [options] <filename(s)>

Options:

--db NAME Name for Kraken 2 DB

(default: none)

--threads NUM Number of threads (default: 1)

--quick Quick operation (use first hit or hits)

--unclassified-out FILENAME

Print unclassified sequences to filename

--classified-out FILENAME

Print classified sequences to filename

--output FILENAME Print output to filename (default: stdout); "-" will

suppress normal output

--confidence FLOAT Confidence score threshold (default: 0.0); must be

in [0, 1].

--minimum-base-quality NUM

Minimum base quality used in classification (def: 0,

only effective with FASTQ input).

--report FILENAME Print a report with aggregrate counts/clade to file

--use-mpa-style With --report, format report output like Kraken 1's

kraken-mpa-report

--report-zero-counts With --report, report counts for ALL taxa, even if

counts are zero

--memory-mapping Avoids loading database into RAM

--paired The filenames provided have paired-end reads

--use-names Print scientific names instead of just taxids

--gzip-compressed Input files are compressed with gzip

--bzip2-compressed Input files are compressed with bzip2

--help Print this message

--version Print version information

If none of the *-compressed flags are specified, and the filename provided

is a regular file, automatic format detection is attempted.

$ kraken2\

--db minikraken2_v2_8GB_201904_UPDATE/ \

--threads 20 \

--report report \

--gzip-compressed --paired \

${B}_1.fastq.gz ${B}_2.fastq.gz

bracken也支持kraken2的結果。

相關閱讀

腸道菌群：16S測序分析流程解讀

腸道菌群：宏基因組測序分析流程解讀

腸道菌群：宏轉錄組測序分析流程解讀

腸道菌群：宏病毒組測序分析流程解讀

猜你喜歡

10000+：菌群分析 寶寶與貓狗 梅毒狂想曲 提DNA發Nature Cell專刊 腸道指揮大腦

系列教程：微生物組入門 Biostar 微生物組  宏基因組

專業技能：學術圖表 高分文章 生信寶典 不可或缺的人

一文讀懂：宏基因組 寄生蟲益處 進化樹

必備技能：提問 搜索  Endnote

文獻閱讀 熱心腸 SemanticScholar Geenmedical

擴增子分析：圖表解讀 分析流程 統計繪圖

16S功能預測   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在線工具：16S預測培養基 生信繪圖

科研經驗：雲筆記  雲協作 公眾號

編程模板: Shell  R Perl

生物科普:  腸道細菌 人體上的生命 生命大躍進  細胞暗戰 人體奧秘  

寫在後面

為鼓勵讀者交流、快速解決科研困難，我們建立了「宏基因組」專業討論群，目前己有國內外5000+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註「姓名-單位-研究方向-職稱/年級」。PI請明示身份，另有海內外微生物相關PI群供大佬合作交流。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍未解決群內討論，問題不私聊，幫助同行。

學習16S擴增子、宏基因組科研思路和分析實戰，關注「宏基因組」

點擊閱讀原文，跳轉最新文章目錄閱讀

相關焦點

  • 多快好省的宏基因組研究技巧
    不通過拼接組裝，直接基於kraken2等kmer，或MetaPhlAn2等標記基因的參考基因組方法進行種屬豐度分類。結合其到菌株的物種分類和豐度數據可較16s方案下的PICRUST更加準確的預測基因構成。
  • 宏基因組的一些坑和解決方案
    資料庫的選擇和完善程度我們使用的是kraken2以及谷禾自建的基於Reseq 99版本的微生物基因組資料庫。左側是使用了完整所有序列的PCA圖，右側是每個樣本隨機抽取100萬reads的PCA圖。可以看到隨機抽樣100萬reads之後的菌群分布情況和完整序列的完全相同。
  • Microbiome:CAMISIM模擬宏基因組和微生物群落
    Github: https://github.com/CAMI-challenge/CAMISIM看了一下使用方法：CAMISIM可以採用兩種方法進行模擬：1是根據分類學文件，從NCBI上完整基因組中進行查找和模擬；結果儘可能和輸入文件相接近。2是de novo模擬。
  • Nature子刊:HUMAnN2實現宏基因組和宏轉錄組種水平功能組成分析
    HUMAnN2是HUMAnN1的升級版，可實現快速宏基因組、宏轉錄組的物種和功能定量，同時提供功能通路內物種組成信息，軟體在文章發表前已經在圈內廣泛使用，關於它的使用教程請參閱公眾號之前關於此軟體的教程和相關流程的文章。
  • 快速看懂腸道菌群宏基因組測序分析報告
    物種注釋metaphlan2主頁：http://segatalab.cibio.unitn.it/tools/metaphlan2/從Clean read出發，使用metaphlan2軟體分析，獲得不同分類層級的物種豐度表。
  • 宏基因組測序和16S rDNA測序有啥區別?
    編者按　16S rDNA測序及宏基因組測序都是微生物研究和應用的重要方法，那麼問題來了：兩者到底有什麼區別呢？什麼情況下需要做16S測序？什麼情況下需要做宏基因組測序？什麼情況下需要二者結合使用呢？那麼在開始宏基因組測序專題前，小編需要給大家解決一個非常重要的問題——16S測序和宏基因組測序的主要區別是什麼？
  • NAR | GMrepo：人類腸道宏基因組資料庫
    對收錄的樣品進行質控，分類注釋和豐度計算，最後將這些信息網頁可視化，以便用戶查詢和重利用人類腸道微生物數據。目前可查詢的表型包括年齡、性別、國家、BMI和近期抗生素使用情況等，用戶還可以通過 GMrepo 獲得預先計算好的物種豐度、表型內和表型間的流行度以及菌株共現網絡等信息。
  • 基於「三+二」宏基因組測序的抗性基因和可移動元件的精確研究
    2、宏基因組組裝：分別使用組裝流程OPERA-MS（圖1）和其他組裝器進行宏基因組組裝，並對組裝集進行質量評估。3、人工構建GIS20腸道微生物群落及其測序：將20種已知的細菌菌株的DNA混合在一起，菌株DNA的豐度為0.1%到30%，再進行三代測序。比較不同工具對該數據的組裝效果。
  • 從宏基因組數據中挖掘、分析和整合病毒信號
    雖然噬菌體代表與健康有關的病毒社區的顯著部分，得到噬菌體有更深的了解是由於病毒分離和純化的困難仍然具有挑戰性。然而，使用宏基因組測序，急劇通過從生態和臨床樣品的宏基因組回收的噬菌體的基因組促進了病毒研究。使用宏基因組數據揭示噬菌體的分類和功能組成以及噬菌體 - 宿主相互作用。此外，可以實現對微生物組的更好的生態學理解以及對其對人類健康的影響的更深入的見解。
  • nanopore宏基因組分析培訓班（第2期）開始報名了
    我們第六期納米孔測序培訓班將專門關注納米孔測序在宏基因組中的應用。本次培訓班，我們將系統介紹納米孔測序在宏基因組中的應用，通過本次培訓班你將完成以下五大主題：1、熟悉納米孔測序；2、掌握納米孔宏基因組測序物種分類；3、宏基因組數據分結果可視化；4、掌握納米孔宏基因組測序基因組拼接；5、從納米孔宏基因組拼接中得到細菌完成圖。
  • 宏基因組方法學研究取得進展
    宏基因組方法學研究取得進展 來源：北京生命科學研究院   發布者：張薦轅   日期：2016-03-23   今日/總瀏覽：1/5100
  • 宏基因組02. HUMAnN2 --宏基因組代謝通路分析
    本系列課程前情回顧宏基因組分析00. 概述和分析環境的部署宏基因組分析01.
  • 利用宏基因組序列來確定蛋白質結構
    在這項研究中，Baker實驗室的蛋白質結構預測伺服器Rosetta分析了JGI集成微生物基因組（IMG）系統上的宏基因組序列。對於蛋白質而言，外觀很重要。當然，這並不是指顏值，而是三維結構。蛋白質是由長的胺基酸鏈組成的，但一維的胺基酸序列似乎沒有意義。只有了解三維結構，研究人員才能弄清蛋白質的結構如何決定它的功能。
  • DADA2中文教程v1.8
    通過注釋序列得到物種信息，下遊我們將使用目前十分流行的R包phyloseq做下遊微生物群落相關分析，盡請期待。這裡我們的16s序列注釋，採用下載Silva參考資料庫進行訓練和注釋，DADA2包為此目的提供了樸素貝葉斯分類器方法實現物種注釋。
  • Nanopore宏基因組測序揭示汙水處理過程中可移動抗生素抗性基因組
    三、研究結果1.質粒和ICEs攜帶的ARGs在WWTPs的抗性組中佔主導地位WWTPs的進水、汙泥、出水的抗性組表現在9個環境宏基因組樣本的1791個攜帶ARG的長reads和316個contigs，基於它們的HGT潛力分為2大類：細胞間可移動組和染色體組
  • 2019微生物組—宏基因組分析專題培訓第三期
    宏基因組簡介發展史、常用技術適用範圍、分析思路22宏基因組有參質控FastQC、Trimmomatic、 MultiQC、KneadData質控、parallel並行計算23物種和功能組成MetaPhlAn2物種組成、HUMAnN2功能組成、功能關聯驅動物種24物種和功能差異比較和可視化LEfSe、STAMP、R語言統計25發表前準備圖片排版、數據釋放、代碼整理(可選)26網絡繪製Cytoscape
  • 利用深度變體自動編碼器改進宏基因組的組裝
    利用深度變體自動編碼器改進宏基因組的組裝 作者：小柯機器人 發布時間：2021/1/5 16:19:03 丹麥哥本哈根大學Simon Rasmussen課題組的最新研究利用深度變體自動編碼器改進了宏基因組的組裝。
  • 人類微生物組研究設計、樣本採集和生物信息分析指南
    表 1: 大腸埃希菌細菌層級分類分類層級分類名稱門變形菌門綱丙型變形菌綱目腸桿菌目科腸桿菌科屬埃希氏桿菌屬種埃希氏菌株EIEC112ac株2.3 操作分類單元和擴增子序列變異操作分類單元（operational taxonomic units，OTUs）的構建對於標記基因（擴增子）數據分析非常重要[17]。
  • 動物所揭開大熊貓消化竹子纖維素和半纖維素之謎
    大熊貓（Ailuropoda melanoleuca）在分類上屬食肉目，具有食肉動物典型的消化系統。但它卻特化以竹類為食，無疑是典型的素食者。竹子是一種高纖維和低營養的食物，已有的研究表明，大熊貓除消化竹子中90%以上的蛋白質和脂肪等營養物質外，還能利用約8%的纖維素和27%的半纖維素。然而，2010年公布的大熊貓基因組卻未能發現大熊貓具有編碼纖維素和半纖維素等消化酶的基因。
  • 三代nanopore宏基因組測序數據分析，北京，11月7-9日
    本次培訓班，我們將系統介紹納米孔測序在宏基因組中的應用，通過本次培訓班你將完成以下四大主題：1、熟悉納米孔測序；2、病原微生物快速鑑定；3、宏基因組數據分結果可視化；4、納米孔建庫測序實驗。課程特色1、真正的小班授課，提供助教，提供一對一指導，零基礎可學習。2.