macでインフォマティクス

macでインフォマティクス

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

KrakenUniq

2019 1/17 エラー修正 

 

 メタゲノミクス分類手法は、データセット内の各リードに taxonomic identityをアサインすることを試みる。メタゲノミクスデータにはしばしば何千万ものリードが含まれているため、分類は、通常、長さk(k-mers)の短いワードの正確な一致を使用して行われる。結果には、リードの分類が含まれているが、アライメントされたゲノム内のポジションは含まれない([論文のref.1]のレビューのように)。しかしながら、リードカウントは欺かれる可能性がある。サンプルの抽出、取り扱い、またはシーケンシング中に実験室キットまたは環境から導入された汚染により、大量の偽のidentificationsが引き起こされる可能性がある [ref.2,3]。サンプル投入が少量だと、汚染の問題をさらに複雑にする可能性がある。たとえば、感染症の臨床診断にシークエンシングを使用する場合、投入DNAノ0.1%未満だけが目的の微生物に由来する可能性がある[ref.4, 5]。ゲノムの複雑さの低い領域やデータベースゲノム自体の汚染から、追加の疑似マッチが生じる可能性がある[ref.6]。

 そのような偽陽性のリードはゲノムの小さな部分にのみ一致して起きる。典型的には、種のゲノムが低複雑度領域を含み、その種に合致する唯一のリードがこの領域に入る場合など起きる。真に存在する微生物からのリードは、1つまたは少数の場所に集中するのではなく、ゲノム全体に比較的均一に分布すべきである。ゲノムアライメントはこの情報を明らかにすることができる。しかし、アライメントはリソースを集中的に使用するため、ゲノムすべてのindexを作成する必要があり、すべてのリードをそれらのindexと比較するための比較的遅いアライメントステップが必要になる。いくつかのメタゲノミクス手法では、カバレッジ情報を使用してマッピング定量の精度を向上させているが、これらの手法では、より遅いアラインメント手法を入力として必要としている[ref.7 link]。アセンブリベースの方法は偽陽性を回避するのにも役立つが、これらは非常に豊富な種にのみ有用である[ref.8 link]。

 ここでは、非常に高速なk-merベースの分類と高速なk-mer濃度推定とを組み合わせた新しい方法であるKrakenUniqを提示する。 KrakenUniqはKrakenメタゲノム分類器[ref.9]に基づいており、また、各taxonについて特定されたユニークなk-mersの数を数えるために、効率的な濃度推定アルゴリズムHyperLogLog [ref.10, 11, 12]を使用する方法が追加されている。 KrakenUniqは、各ゲノム中のユニークなk-mersのうち、どれだけがリードでカバーされているかを数えることによって、true-positiveのマッチからfalse-positiveを識別することができる。さらに、KrakenUniqは、メタゲノム分類を改善するための追加の新機能を実装している。(a)複数のデータベースに対して階層的に検索できる。 (b)taxonomy は、系統およびプラスミドのためのノードを含むように拡張され、それによりそれらの検出が可能になる。 (c)database build scriptは、NCBIウイルスゲノムリソース[ref.13]から10万を超えるウイルスの追加を可能にする。KrakenUniqは、Krakenによって提供される情報のスーパーセットを提供し、同時に、同じくらい早く、または少し高速に実行でき、分類時に追加メモリをほとんど使用しない。

 KrakenUniqは、メタゲノム実験で同定されたすべての分類群の効率的なk-merカウント情報を提供するために開発された。 主なワークフローは次のとおりである。リードが処理されると、各k-merにデータベースからtaxonがアサインされる(論文 図1a)。 KrakenUniqは、各分類群のHyperLogLogデータスケッチをインスタンス化し、それにk-mersを追加する(論文 図1bおよび追加ファイル1:HyperLogLogアルゴリズムのセクション1)。リードを分類した後、KrakenUniqはtaxonomic treeを走査し、各taxonの見積りをその親とマージする。 分類レポート出力には、入力データで観測された各taxonの固有のk-mer数とk-merカバレッジデプスが含まれている(論文 図1c)。

(以下省略)

 

f:id:kazumaxneo:20190115204013p:plain図1. Overview of the KrakenUniq algorithm and output.  論文より転載

 

KrakenUniqに関するツイート

 

インストール

ubuntu14.04、miniconda3-4.0.5環境で動作確認を行った。

依存

Jellyfish ver.2には対応していないので注意する。

ハードウエア要件

  • KrakenUniq requires a lot of RAM - ideally 128GB - 512GB. For more memory efficient classification consider using centrifuge.

本体 Github

#condaを使って導入できる(linuxのみ)
conda install -y -c bioconda krakenuniq

> krakenuniq --help

$ krakenuniq --help

Usage: krakenuniq --report-file FILENAME [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

  --gzip-compressed       Input is gzip compressed

  --bzip2-compressed      Input is bzip2 compressed

  --hll-precision INT     Precision for HyperLogLog k-mer cardinality estimation, between 10 and 18 (default: 12)

  --exact                 Compute exact cardinality instead of estimate (slower, requires memory proportional to cardinality!)

  --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

  --output FILENAME       Print output to filename (default: stdout); "off" 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

 

Experimental:

  --uid-mapping           Map using UID database

 

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

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

> krakenuniq-download --help

$ krakenuniq-download --help

Unknown database --help. 

 

krakenuniq-download [<options>] <pattern> <pattern>*

 

ARGUMENTS

 <pattern> can be one of

     'contaminants'     Contaminant sequences from UniVec and EmVec.

     'taxonomy'         NCBI taxonomy mappings from ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/

     'nucleotide'       Download nucleotide sequences using a query specified using --search or --ac.

     'microbial-nt'     Download microbial sequences from nt database.

     'viral-neighbors'  Download viral strain sequences from the NCBI Viral Genome Resource.

                        (Search: "viruses[organism] not cellular organisms[orgn] not wgs[prop] not gbdiv syn[prop] and (nuccore genome samespecies[filter] and (complete[title]) not unverified[title])").

     'genbank/DOMAIN'   Download all complete genomes for DOMAIN from GenBank.

     'refseq/DOMAIN'    Download all complete genomes for DOMAIN from RefSeq.

     'refseq/DOMAIN/ASS_LEVEL'

     'refseq/DOMAIN/ASS_LEVEL/COLUMN=value1(/COLUMN=value2)*' 

        Possible values for DOMAIN: bacteria viral archaea fungi protozoa invertebrate plant vertebrate_mammalian vertebrate_other.

        Possible values for ASS_LEVEL: Any, Complete_Genome, Chromosome, Scaffold and Contig. 

        Possible values for COLUMN: Any column in the NCBI assembly_summary.txt, e.g. species_taxid or assembly_accession.

        Example: 'refseq/vertebrate_mammalian/Any/species_taxid=9606' <- download all human assemblies

 

COMMON OPTIONS

 -o <directory>     Folder to which the files are downloaded. Default: '.'

 --db <directory>   Alternative to -o: Download to <directory>/{library,taxonomy}.

 --threads <# of threads>  Number of processes when downloading (uses xargs). Default: '5'

 --rsync, -R        Download using rsync.

 --overwrite        Redownload and overwrite files with the same name.

 --verbose          Be verbose.

 --dust, -D         Mask low-complexity regions using dustmasker.

 --min-seq-len X    Filter all sequences from the FASTA files that have less than X bp.

 

WHEN USING DATABASE nucleotide or viral-neighbors:

 --search SEARCH_TERM  Download all sequences returned from a NCBI nucleotide search.

                       When used with viral-neighbors, it subsets the viral genomes by the search.

                       E.g. "txid1570291[Organism]" for Ebola virus sequences (taxonomy ID 1570291).

 --ac AC1,AC2          Alternative to --search. Download specified ACs of nucleotide database.

 --mapping-file MAP    Map accessions using the specified mapping file(s) (comma-separated).

                       Possible values: nucl_est, nucl_gb, nucl_gss, nucl_wgs.

                       For viral-neighbors, the default is nucl_gb. Unset by giving it an empty string.

                       Downloaded from ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/.

 --retmode RM          Specify return mode. Possible options: text (default), asn.1, xml.

 --rettype RT          Specify return type for download. Note that no mapping files are generated when

                       rettype is specified. Possible options: fasta (default), gb, gbc, native, acc, seqid, ft, 

                       gbwithparts, fasta_cds_na, fasta_cds_aa. Note that only gb and fasta files are split, while the other formats stay in chunks.

                       The resulting files will have the rettype as extension.

 

WHEN USING DATABASE nt:

 --taxa TAXLIST     Provide list of taxa that are kept in the database. Comma separated, supply either taxIDs or

                    one of the following division names: archaea bacteria fungi parasitic_worms protozoa viral

                    Default: "bacteria,archaea,viral,fungi,protozoa"

 --exclude-environmental-taxa

                    Exclude taxa that are named 'environmental samples'.

 

WHEN USING DATABASE refseq OR genbank:

 --fna <seq types>  Comma-separated list of sequence types, including genomic, rna, rna_from_genomic, cds_from_genomic. Default: genomic.

                    See the assembly project FTP site for available sequences

 -u                 Filter unplaced sequences.

> krakenuniq-build --help

$ krakenuniq-build --help

Usage: krakenuniq-build [task option] [options]

 

Task options (exactly one can be selected -- default is build):

  --download-taxonomy        Download NCBI taxonomic information

  --download-library TYPE    Download partial library (TYPE = one of "refseq/bacteria", "refseq/archaea", "refseq/viral"). 

                             Use krakenuniq-download for more options.

  --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)

  --rebuild                  Create DB from library like --build, but remove

                             existing non-library/taxonomy files before build

  --clean                    Remove unneeded files from a built database

  --shrink NEW_CT            Shrink an existing DB to have only NEW_CT k-mers

  --standard                 Download and create default database, which contains complete genomes 

                             for archaea, bacteria and viruses from RefSeq, as well as viral strains 

                             from NCBI. Specify --taxids-for-genomes and --taxids-for-sequences

                             separately, if desired.

 

  --help                     Print this message

  --version                  Print version information

 

Options:

  --db DBDIR                 Kraken DB directory (mandatory except for --help/--version)

  --threads #                Number of threads (def: 1)

  --new-db NAME              New Kraken DB name (shrink task only; mandatory

                             for shrink task)

  --kmer-len NUM             K-mer length in bp (build/shrink tasks only;

                             def: 31)

  --minimizer-len NUM        Minimizer length in bp (build/shrink tasks only;

                             def: 15)

  --jellyfish-hash-size STR  Pass a specific hash size argument to jellyfish

                             when building database (build task only)

  --jellyfish-bin STR        Use STR as Jellyfish 1 binary.

  --max-db-size SIZE         Shrink the DB before full build, making sure

                             database and index together use <= SIZE gigabytes

                             (build task only)

  --shrink-block-offset NUM  When shrinking, select the k-mer that is NUM

                             positions from the end of a block of k-mers

                             (default: 1)

  --work-on-disk             Perform most operations on disk rather than in

                             RAM (will slow down build in most cases)

  --taxids-for-genomes       Add taxonomy IDs (starting with 1 billion) for genomes.

                             Only works with 3-column seqid2taxid map with third 

                             column being the name

  --taxids-for-sequences     Add taxonomy IDs for sequences, starting with 1 billion.

                             Can be useful to resolve classifications with multiple genomes

                             for one taxonomy ID.

  --min-contig-size NUM      Minimum contig size for inclusion in database.

                             Use with draft genomes to reduce contamination, e.g. with values between 1000 and 10000.

  --library-dir DIR          Use DIR for reference sequences instead of DBDIR/library.

  --taxonomy-dir DIR         Use DIR for taxonomy instead of DBDIR/taxonomy.

 

Experimental:

  --uid-database             Build a UID database (default no)

  --lca-database             Build a LCA database (default yes)

  --no-lca-database          Do not build a LCA database

  --lca-order DIR1           Impose a hierarchical order for setting LCAs.

  --lca-order DIR2           The directories must be specified relative to the libary directory

  ...                        (DBDIR/library). When setting the LCAs, k-mers from sequences in

                             DIR1 will be set first, and only unset k-mers will be set from

                             DIR2, etc, and final from the whole library.

Use this option when including low-confidence draft genomes,

                             e.g use --lca-order Complete_Genome --lca-order Chromosome to

                             prioritize more complete assemblies.

                             Keep in mind that this option takes considerably longer.

> krakenuniq-extract-reads --help

$ krakenuniq-extract-reads --help

/home/uesaka/.pyenv/versions/miniconda3-4.2.12/bin/krakenuniq-extract-reads version [unknown] calling Getopt::Std::getopts (version 1.12 [paranoid]),

running under Perl version 5.26.2.

 

Usage: krakenuniq-extract-reads [-OPTIONS [-MORE_OPTIONS]] [--] [PROGRAM_ARG1 ...]

 

The following single-character options are accepted:

With arguments: -t

Boolean (without arguments): -f -v -i -a -p

 

Options may be merged together.  -- stops processing of options.

Space is not required between options and their arguments.

  [Now continuing due to backward compatibility and excessive paranoia.

   See 'perldoc Getopt::Std' about $Getopt::Std::STANDARD_HELP_VERSION.]

 

krakenuniq-extract-reads: Extract all reads from FASTQ file that are matched to a specied taxon by KrakenUniq

 

Usage: krakenuniq-extract-reads [OPTIONS] <taxon> <kraken> <fasta/fastq>

 

<taxon>         taxonomy ID, possibly multiple separated by ','

<kraken>        kraken result file

<fasta/fastq>   fasta/fastq file, possibly gzipped

 

Options:

  -a  input is FASTA file (default: FASTQ)

  -f  output in FASTA format

  -i  invert: print all reads not matching taxon

  -t TAXDB Include children of taxonomy IDs, using TAXDB to find them

  -v  verbose

  -p  paired-end reads: use a '%' in fasta/q file name as placeholder for 1 and 2

 

Example:

    krakenuniq-extract-reads -p 9606 result.kraken input_%.fq 

    outputs all reads of input_1.fq and input_2.fq that have the taxonomy ID 9606 to STDOUT.

> krakenuniq-report --help

$ krakenuniq-report --help

Usage: krakenuniq-report --db KRAKEN_DB_NAME [OPTIONS] <kraken output file(s)>

 

OPTIONS:

  --show-zeros    Show full taxonomy table.

  --taxon-counts  Input files are in the format '<taxon ID><tab><count>' instead of Kraken output.

  --taxon-list    Input files is list of taxon IDs instead of Kraken output.

  -h              This message.

  

This script should only be used when post-processing raw KrakenUniq output, and k-mer counts and coverages are not needed. For most use-cases, krakenuniq --report-file is better than krakenuniq-report. 

> krakenuniq-mpa-report --help

$ krakenuniq-mpa-report --help

Usage: krakenuniq-mpa-report [--db KRAKEN_DB_NAME] [options] <kraken output file(s)>

 

Options:

  --db NAME             Name of Kraken database

                        (default: none)

  --show-zeros          Display taxa even if they lack a read in any sample

  --header-line         Display a header line indicating sample IDs

                        (sample IDs are the filenames)

  --intermediate-ranks  Display taxa not at the standard ranks with x__ prefix

> krakenuniq-filter --help

$ krakenuniq-filter --help

Usage: krakenuniq-filter [--db KRAKEN_DB_NAME] [--threshold NUM] <kraken output file(s)>

 

Threshold must be between 0 and 1.

> krakenuniq-translate --help

$ krakenuniq-translate --help

Usage: krakenuniq-translate [--db KRAKEN_DB_NAME] [--mpa-format] <kraken output file(s)>

 

 

データベースのビルド

1、taxonomyのダウンロード (350MBくらい)

krakenuniq-download --db DBDIR taxonomy

$ ls -alh taxonomy/

total 345M

drwxr-xr-x 2 uesaka user 4.0K Jan 14 10:27 .

drwxr-xr-x 4 uesaka user 4.0K Jan 14 10:35 ..

-rw-r--r-- 1 uesaka user 167M Jan 14 10:23 names.dmp

-rw-r--r-- 1 uesaka user 133M Jan 14 10:23 nodes.dmp

-rw-r--r-- 1 uesaka user  45M Jan 14 10:23 taxdump.tar.gz

-rw-r--r-- 1 uesaka user   29 Jan 14 10:27 timestamp

2、データベースのダウンロード

#Choose one of the following: bacteria, viral, archaea, fungi, protozoa, invertebrate, plant, vertebrate_mammalian, vertebrate_other


#REgseqの全bacteria、archaea
krakenuniq-download --db DBDIR --threads 20 --dust refseq/bacteria refseq/archaea

#plant(3種のみ コンタミチェックに使う)
krakenuniq-download --db DBDIR_plant --threads 20 --dust refseq/plant

#fungi (10種)
krakenuniq-download --db DBDIR_fungi --threads 20 --dust refseq/invertebrate

#protozoa (3種) 
krakenuniq-download --db DBDIR_protozoa --threads 20 --dust refseq/protozoa

#invertebrate
krakenuniq-download --db DBDIR_invertebrate --dust refseq/invertebrate


## Contaminant sequences from UniVec and EmVec, plus the human reference genome
krakenuniq-download --db DBDIR refseq/vertegrate_mammalian/Chromosome/species_taxid=9606

## All viral genomes from RefSeq plus viral 'neighbors' in NCBI Nucleotide
krakenuniq-download --db DBDIR_virus --threads 20 refseq/viral/Any viral-neighbors

## All microbial (including eukaryotes) sequences in the NCBI nt database (並列処理に対応していない?)
krakenuniq-download --db DBDIR --dust microbial-nt

#
krakenuniq-download --db DBDIR --threads 20 \
--taxa "archaea,bacteria,viral,fungi,protozoa,helminths" \
--dust --exclude-environmental-taxa nt

#DB/にnt~fnaファイルがダウンロードされる。
  • --dust    Mask low-complexity regions using dustmasker.
  •  --db <directory>     Alternative to -o: Download to <directory>/{library,taxonomy}.
  • --threads <NUM>    Number of processes when downloading (uses xargs). Default: '5'
  • --taxa <taxlist>     Provide list of taxa that are kept in the database. Comma separated, supply either taxIDs or
  • --exclude-environmental-taxa   Exclude taxa that are named 'environmental samples'.

(refseq/のplant、fhungi、protozoa、bacteria、archaea、viral/Any、 viral-neighborsのみ動作確認を行った)

 

 3、データベースビルド(時間がかかる)

krakenuniq-build --db DBDIR --kmer-len 31 --threads 40 --taxids-for-genomes --taxids-for-sequences
  • --kmer-len <NUM>     K-mer length in bp (build/shrink tasks only; def: 31)

  • --taxids-for-genomes       Add taxonomy IDs (starting with 1 billion) for genomes.
  • --taxids-for-sequences    Add taxonomy IDs for sequences, starting with 1 billion.

 

 

実行方法

レポートファイル作成まではkrakenと同じ流れで進める。

#メモリに余裕があれば、データベースをmemoryにpreloadしてインメモリデータベースにしておくことでかなり早くなる
krakenuniq --db DBDIR/ --preload --threads 20

#k-merアサインを実行
krakenuniq --threads 20 --db DBDIR/ --fastq-input input1.fq input2.fq > kmer-assignment
  • --preload    Loads DB into memory before classification 

 > head -n 5 kmer-assignment

U M02077:12:000000000-A8737:1:1101:26651:6691 0 300 0:270

C M02077:12:000000000-A8737:1:1101:7099:6692 51589 300 0:124 2237:5 0:36 2237:12 0:31 51589:8 0:31 2237:23

U M02077:12:000000000-A8737:1:1101:27349:6692 0 299 0:269

U M02077:12:000000000-A8737:1:1101:25832:6693 0 300 0:270

リードごとにtaxonomyアサインされる(raw KrakenUniq output)。

 

レポート作成

krakenuniq-report --db DBDIR/ kmer-assignment > report
  • --show-zeros    Show full taxonomy table

  

 

すでにkrakenのデータベースがあるなら、以下のコマンドでtaxDBを作成するだけでランできる。 

cd KRAKEN_DB/
#DBディレクトリで次のコマンドを打ち、DB/にtaxDBを追加
build_taxdb taxonomy/names.dmp taxonomy/nodes.dmp > taxDB

krakenuniq --report-file FILENAME --preload --db KRAKEN_DB/\
--threads 20 --fastq-input input.fq

 

引用

KrakenUniq: confident and fast metagenomics classification using unique k-mer counts

F. P. Breitwieser, D. N. Baker and S. L. Salzberg

Genome Biology 201819:198

 

関連