2021/8/17 論文引用
2022/02/17 krona追記
2022/02/24 krona関係のコマンド修正
2023/12/01 説明追加
ヒトの腸内細菌叢の特徴を明らかにするために、微生物の分離とDNA配列の決定の両方が行われてきた。また，最新のバイオインフォマティクスツールを用いて，新規に構築されたゲノム（Metagenome-Assembled Genomes-MAGs）も大きな貢献をしてきた。その結果，最近，200,000以上の非冗長な参照ゲノムからなる Unified Human Gastrointestinal Genome（UHGG）コレクション（ref.7）が発表され，この分野における大きな節目となった。これらの研究は，ヒトの腸内で見られる膨大な種類のゲノムを特定し，確固たる基礎を築いた。しかし、これらの研究では、健康な人の腸内に存在するゲノムの世界的な普及状況、すなわち、その出現頻度に関する情報を提供していない。このような知識は、世界中の健康なヒト腸内マイクロバイオームを反映したヒト腸内関連原核生物ゲノムのコレクションを設定するために不可欠である。これは、ヒト消化管内部のマイクロバイオーム研究における比較研究に使用することを目的としたカスタムデータベースを構築するために特に重要である。
地域的には，腸内細菌叢は環境によって大きく形成されていること[ref.8]，その組成がさまざまな疾患や障害と関連していること[ref.9,10,11,12]が研究によって明らかにされており，したがって，現在，腸内細菌叢の治療的介入が導入されている段階にある[ref.13,14]。しかし，健康なヒトの腸内細菌叢に関する世界的な基準がないことがボトルネックとなっている[ref.15]。このことが，世界規模での腸内細菌叢の理解と，大規模な介入戦略の導入の両方を妨げている。
本研究の目的は，すべてのヒト腸内細菌叢研究のための普遍的な基準として，健康なヒトに関連する腸内細菌の単一かつ包括的なゲノムコレクション（HumGut）を作成することであった。NCBIのRefSeqゲノムとともに、上述のUHGGコレクションを利用した。HumGutを構築するための戦略を論文図1に示す。
HumGutのゲノムは、世界中で収集された健康なヒトの腸内メタゲノムに含まれる量によってランク付けされている。最もよく遭遇するゲノム（つまりリストの上位）は、dereplicationの際に代表的な分類として選択され、最も関連性の高いリストを確保している。比較的単純なコンセプトのように思われるかもしれないが、この研究は、増え続ける原核生物のゲノムのために、一般に公開されているヒトの腸内メタゲノムを迅速にスクリーニングできるバイオインフォマティクスのツールが開発されたことによって初めて可能になった。
5,700本以上の健康なヒト腸内メタゲノムをスクリーニングし、RefSeqと最近発表されたUHGGコレクションから得られた490,000本以上の一般に入手可能な原核生物ゲノムが含まれていることを確認した。その結果、381,000個以上のゲノムが得られ、これらのゲノムは健常人のメタゲノムに含まれる頻度に基づいてスコアリングされ、ランク付けされた。次に、これらのゲノムを97.5%の配列同一性の分解能でクラスター化し、クラスターを代表するもの（合計30,691個）をHumGutコレクションとして保存した。Kraken2を用いた分類では、メタゲノム中のリードの平均94.5%を分類できたが、UHGGでは86%、標準的なKraken2データベースでは44%であった。UHGGと同様に95%の配列同一性で複製されたゲノムからなる粗いHumGutコレクションでは、リードの88.25%を分類した。HumGutは、標準的なKraken2データベースの半分のサイズで、UHGGのサイズと同程度だが、両者を凌駕している。
HumGutコレクションには、97.5%の配列同一性の分解能でクラスタリングされた30,000以上のゲノムが含まれており、ヒト腸内の有病率でランク付けされている。IBD患者からのメタゲノムも同様にこのコレクションにマッピングされることを示し、このリファレンスがHumGutの取得に使用されたメタゲノムリファレンスセット以外の研究にも適していることを示している。すべてのデータとメタデータ、および有用なコードは、http://arken.nmbu.no/~larssn/humgut/で公開されている。
ゲノムデータの公開サイト
http://arken.nmbu.no/~larssn/humgut/
Github
#kraken2とbrackenのインストール
#conda (link)
mamba create -n kraken2 -y
conda activate kraken2
mamba install -c bioconda -y kraken2
#bracken (link)
mamba install -c bioconda -y bracken
データベースのビルド
Githubの説明に則ってkraken2のデータベースを構築する流れを確認する。
１、KRAKEN2データベースを構築するフォルダを作成し、その中に上記HPからダウンロードしたGTDBかNCBIのnames.dmpとnodes.dmpファイルを配置する。
mkdir -p KRAKEN2/taxonomy
mv gtdb_names.dmp KRAKEN2/taxonomy/names.dmp
mv gtdb_nodes.dmp KRAKEN2/taxonomy/nodes.dmp
２、ヒトの腸からのデータであれば、ホストゲノムがどの程度検出されるか調べられるように、ヒトゲノムもデータベースに含めることが推奨されている。kraken2のコマンドを使ってダウンロードできる。
kraken2-build --download-library human --db KRAKEN2/
- --download-library Download partial library (TYPE = one of "archaea", "bacteria", "plasmid", "viral", "human", "fungi", "plant", "protozoa", "nr", "nt", "UniVec", "UniVec_Core")
KRAKEN2/libraryに関連ファイルと共にダウンロードされる。テストした時は20分ほどかかった。
３、上記HPからダウンロードした全ゲノムのtar ballを解凍する。
tar -xvf HumGut.tar.gz
fna/にgzipped fastaファイルが展開される。
４、１つのファイルに結合する。kraken2は圧縮ゲノムfastaからのデータベースのビルドに対応していないので、zcatを使って結合する。
zcat fna/*.fna.gz > HumGut975_library.fna
結合後のfastaファイルは63GBとなった (2021年8月現在)。
5、kraken2-buildコマンドで HumGut975_library.fnaファイルを追加する。
kraken2-build --add-to-library HumGut975_library.fna --db KRAKEN2/
- --add-to-library Add FILE to library
5のステップの後、上記で作成した HumGut975_library.fnaは削除してO.K。
補足
Githubでは、５のステップの前に、kraken2の大きなメモリ使用量を減らすため、クラスタリングの配列同一性の閾値を変更してゲノムの数を減らす手順が説明されています。
６−１、kraken2のデータベースをビルドする。
#先にgtdb_nodesとgtdb_namesのdumpファイルをtaxonomyに配置しておく
mkdir KRAKEN2/taxonomy
mv gtdb_nodes.dmp KRAKEN2/taxonomy/nodes.dmp
mv gtdb_names.dmp KRAKEN2/taxonomy/names.dmp
#build
kraken2-build --build --threads 40 --db KRAKEN2/
このステップが一番時間がかかる。試した時は２−３時間かかった。ビルドが完了すると、KRAKEN2フォルダ内にhash.k2d、opts.k2d、taxo.k2d、seqid2taxid.mapのファイルができる。
６−２、Kraken 2データベースの中身を確認する。
kraken2-inspect --db KRAKEN2 | head -n 20
kraken2-inspectの出力形式は、kraken2の--reportオプションで生成されるレポートと同じ。
７、kraken2データベースができたら、これを拡張してbrackenデータベースも構築する。40スレッド、k-mer長35（kraken2と同じ）、read length 100。
bracken-build -d KRAKEN2 -t 40 -k 35 -l 100
KRAKEN2フォルダにdatabase.kraken, database.100mers.kraken, database.100mers.kmer_distribの3つのファイルが追加される。brackenはkrakenの曖昧なk-merマッチングを改善するもので、よく似た複数マッチの場合に、各ゲノムのユニーク領域からのユニークマッチの数を使って、ベイズ推定して種レベルの確率的な存在量を推定する。krakenや他のtaxonomic プロファイリングツールの精度を上げる目的で使われる。
レポジトリではkraken2/brackenの代わりにkrakenUniqのデータベースを構築する流れも説明されています。興味がある方は確認してみて下さい。
1, kraken2のラン（manual）
ペアエンドfastq。gzip圧縮なら--gzip-compressedを付ける。
kraken2 --db KRAKEN2 --threads 40 --gzip-compressed --paired --report report.txt --use-names --classified-out seqs#.fq input_R1.fq.gz input_R2.fq.gz > output.txt
- --minimum-hit-groups <NUM> Minimum number of hit groups (overlapping k-mers
sharing the same minimizer) needed to make a call (default: 2)
- --unclassified-out <FILENAME> Print unclassified sequences to filename
- --classified-out <FILENAME> Print classified sequences to filename
- --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
- --use-names Print scientific names instead of just taxids
全データを報告する標準出力（> output.txt）はここでは使わない。不要なら> /dev/nullでもO.K。
krakenのレポートファイルは以下のフォーマットになっている（マニュアルより）。
- Percentage of reads covered by the clade rooted at this taxon
- Number of reads covered by the clade rooted at this taxon
- Number of reads assigned directly to this taxon
- A rank code, indicating (U)nclassified, (D)omain, (K)ingdom, (P)hylum, (C)lass, (O)rder, (F)amily, (G)enus, or (S)pecies. All other ranks are simply '-'.
- NCBI taxonomy ID
- indented scientific name
2, Brackenのラン
kraken report ファイル（krakenの--report出力の方）を指定する。種レベルの推定（-l S）。
bracken -d KRAKEN2 -i kraken_report.txt -r 100 -l S -o bracken
- -i location of the built Kraken 1.0 or Kraken 2.0 database
- -r read length
- -t Default = 10. For species classification, any species with <= 10 (or otherwise specified) reads will not receive any additional reads from higher taxonomy levels when distributing reads for abundance estimation. If another classification level is specified, thresholding will occur at that level.
- -l Default = 'S'. This specifies that abundance estimation will calculate estimated reads for each species. Other possible options are K (kingdom level), P (phylum), C (class), O (order), F (family), and G (genus).
kraken_report_braken_species.txtとbrackenが出力される。kronaではkraken_report_braken_species.txtを使う。
３、kraken2のreportを扱うためのユーティリティ；Kraken Toolsが公開されている。
git clone https://github.com/jenniferlu717/KrakenTools.git
python KrakenTools/kreport2krona.py -r report_bracken_species.txt -o MYSAMPLE.bracken
ktImportText MYSAMPLE.bracken -o MYSAMPLE.bracken.html
#複数ファイル
ktImportText MYSAMPLE.bracken* -o all.html
kronaで視覚化された MYSAMPLE.bracken.htmlが得られる。
参考；kronaのインストールとランについて
ここではkreport2kronaスクリプトとktImportText を使って視覚化したが、それにはkronaをインストールしておく必要がある。
#インストールとデータベースの準備
mamba install -c bioconda -y krona
rm -rf $HOME/miniconda3/envs/kraken2/opt/krona/taxonomy
mkdir -p krona/taxonomy
ktUpdateTaxonomy.sh krona/taxonomy/
cp -r krona/taxonomy/ $HOME/miniconda3/envs/kraken2/opt/krona/
#kronaのラン
ktImportTaxonomy -q 2 -t 3 report_bracken_species.txt -o krona.html
#全サンプル(1-5)
ktImportTaxonomy -q 2 -t 3 *_bracken_species.txt -o krona.html
HumGutのメタデータファイルはGithubとHPから入手できます。
引用
HumGut: a comprehensive human gut prokaryotic genomes collection filtered by metagenome data
Pranvera Hiseni, Knut Rudi, Robert C Wilson, Finn Terje Hegge, Lars Snipen
Microbiome. 2021 Jul 31;9(1):165
2021 8/17
HumGut: a comprehensive human gut prokaryotic genomes collection filtered by metagenome data
Pranvera Hiseni, Knut Rudi, Robert C. Wilson, Finn Terje Hegge & Lars Snipen
Microbiome volume 9, Article number: 165 (2021)
#help
> kraken2-build
$ kraken2-build
Must select a task option.
Usage: kraken2-build [task option] [options]
Task options (exactly one must be selected):
--download-taxonomy Download NCBI taxonomic information
--download-library TYPE Download partial library
(TYPE = one of "archaea", "bacteria", "plasmid",
"viral", "human", "fungi", "plant", "protozoa",
"nr", "nt", "UniVec", "UniVec_Core")
--special TYPE Download and build a special database
(TYPE = one of "greengenes", "silva", "rdp")
--add-to-library FILE Add FILE to library
--build Create DB from library
(requires taxonomy d/l'ed and at least one file
in library)
--clean Remove unneeded files from a built database
--standard Download and build default database
--help Print this message
--version Print version information
Options:
--db NAME Kraken 2 DB name (mandatory except for
--help/--version)
--threads # Number of threads (def: 1)
--kmer-len NUM K-mer length in bp/aa (build task only;
def: 35 nt, 15 aa)
--minimizer-len NUM Minimizer length in bp/aa (build task only;
def: 31 nt, 12 aa)
--minimizer-spaces NUM Number of characters in minimizer that are
ignored in comparisons (build task only;
def: 7 nt, 0 aa)
--protein Build a protein database for translated search
--no-masking Used with --standard/--download-library/
--add-to-library to avoid masking low-complexity
sequences prior to building; masking requires
dustmasker or segmasker to be installed in PATH,
which some users might not have.
--max-db-size NUM Maximum number of bytes for Kraken 2 hash table;
if the estimator determines more would normally be
needed, the reference library will be downsampled
to fit. (Used with --build/--standard/--special)
--use-ftp Use FTP for downloading instead of RSYNC; used with
--download-library/--download-taxonomy/--standard.
--skip-maps Avoids downloading accession number to taxid maps,
used with --download-taxonomy.
--load-factor FRAC Proportion of the hash table to be populated
(build task only; def: 0.7, must be
between 0 and 1).
--fast-build Do not require database to be deterministically
built when using multiple threads. This is faster,
but does introduce variability in minimizer/LCA
pairs. Used with --build and --standard options.
> kraken2-inspect -h
# kraken2-inspect -h
Usage: kraken2-inspect [options]
Options:
--db NAME Name for Kraken 2 DB
(default: none)
--threads NUM Number of threads to use
--skip-counts Only print database summary statistics
--use-mpa-style Format output like Kraken 1's kraken-mpa-report
--report-zero-counts Report counts for ALL taxa, even if
counts are zero
--help Print this message
--version Print version information
> kraken2
# 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.
> bracken -h
# bracken -h
/root/miniconda3/bin/bracken: illegal option -- h
Usage: bracken -d MY_DB -i INPUT -o OUTPUT -w OUTREPORT -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
OUTREPORT New Kraken REPORT output file with Bracken read estimates
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,S1,etc] (default: S)
THRESHOLD number of reads required PRIOR to abundance estimation to perform reestimation (default: 0)
> bracken-build -h
/home/kazu/miniconda3/envs/kraken2/bin/bracken-build: illegal option -- h
Usage: bracken_build -k KMER_LEN -l READ_LEN -d MY_DB -x K_INSTALLATION -t THREADS
KMER_LEN kmer length used to build the kraken database (default: 35)
THREADS the number of threads to use when running kraken classification and the bracken scripts
READ_LEN read length to get all classifications for (default: 100)
MY_DB location of Kraken database
K_INSTALLATION location of the installed kraken/kraken-build scripts (default assumes scripts can be run from the user path)
**Note that this script will try to use kraken2 as default. If kraken2 is not installed, kraken will be used instead
関連
* kraken1とkraken2のデータベースに互換性はないので注意。
参考
kraken2のオリジナルのビルド済みデータベース
