macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

ヒト腸内細菌のゲノムコレクション HumGut

2021 ヒト腸内 (human gut) download database Microbiome (Journal)

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のオリジナルのビルド済みデータベース