Archives

Metagenome序列分类之蛋白质空间搜索: Kaiju

今天开始会介绍一系列Metagenome(以后使用元基因组)数据分析中的序列分类问题,这类问题一直是研究的热点和难点。

热点:有用,比如直接鉴定环境样本(粪便、痰液、肺积液)病原,算法多样性, 从基于碱基组成统计分类、基于核酸Kmer查找、基于氨基酸序列MEM查找到序列比对,各种算法八仙过海、各显神通。

难点:序列分类准确性、敏感性、内存有效性、计算密集性等都是待解决问题。

Kaiju是一款基于在氨基酸空间进行相似性搜索的序列分类算法。

1. 算法介绍

Kaiju算法发表在了Nature Communications杂志 :Fast and sensitive taxonomic classification for metagenomics with Kaiju

Kaiju算法

图1. Kaiju算法

Kaiju算法可以通过上图进行解读,Kaiju两种模式,基于MUM(MaximUm exact Matches)模式和基于Greedy Score模式。

MUM模式:

1. 对输入的核酸序列进行六框翻译,遇到终止密码后切断生成多条ORF序列,输入序列可为氨基酸序列;
2. 排序,按照ORF序列长度排序;
3. 数据库搜索,对氨基酸序列进行BWT变换和FM索引,搜索每个ORF框并获取最长MEM(最小MEM长度,默认为11),如后
   续的ORF不可能获取更长的MEM,终止搜索,使用具有最长MEM的命中序列对序列进行分类,如搜索到具有多个相同长度
   MEMs,使用LCA回溯。

Greedy模式:

1. 对输入的核酸序列进行六框翻译,遇到终止密码后切断生成多条ORF序列,可输入氨基酸序列;
2. 排序,按照BLOSUM62值排序;
3. 数据库搜索,对氨基酸序列进行BWT变换和FM索引,搜索每个ORF的MEM(Greedy模式,最小MEM长度设置为7)并获取最
   高的Score值(MEM向左侧延伸直至最大允许错配个数或者最左侧,比如官方Web服务设置为5, 最小匹配Score值为75),
   如后续的ORF不可能获取更高的Score值,终止搜索,使用具有最高Score值的命中序列对序列进行分类,如搜索到具有
   多个相同Score值的命中,使用LCA回溯。

几点说明:

  1. BLOSUM62 矩阵, ORF Score 以及后续的MEM左侧延伸后Score值计算都是按照 BLOSUM62计算。
  2. Kaiju的Greedy模式敏感度要比MEM模式高,但是牺牲了分类速度,对于一些新的基因尤为明显,Greedy采用的模式是启发式算法,只对MEM的末端延伸。
  3. LCA算法见下图,从所有叶子节点向上回溯,找到所有命中节点的共同祖先,比如a|b|c叶子节点的LCA节点就是x节点。

LCA

图2. LCA回溯

2. 安装

2.1 索引文件

Kaiju 源代码存放在Github : Fast taxonomic classification of metagenomic sequencing reads using a protein reference database

Kaiju 提供编译后的可执行文件,可以直接下载:

执行Kaiju 需要下载库文件(BWT索引文件),并提供了以下几个库(2017-05-16版本):

1. RefSeq
25M 氨基酸序列 (来自7,065细菌和古菌基因组和9,334病毒基因组) 索引大小:7.7GB

2.proGenomes
20M 氨基酸序列, 以及9,334病毒基因组序列,索引大小(9.1 GB)

3. NCBI BLAST nr
94M 氨基酸序列 (NCBI NR 数据库中的细菌、古菌和病毒) ,索引大小(23 GB)

4. NCBI BLAST nr+euk
103M 氨基酸序列 (NCBI NR 数据库中的细菌、古菌和病毒、真菌以及微真核生物) ,索引大小(28 GB)

几点说明:

1.微真核生物,内容列表可以参考:微真核生物列表 2.progenomes来自 http://progenomes.embl.de/, 目前包含了5,000个种,25,038个注释的基因组, Kaiju使用的是代表性基因组,类似Refseq。

数据下载:

wget http://kaiju.binf.ku.dk/database/kaiju_index.tgz
wget http://kaiju.binf.ku.dk/database/kaiju_index_pg.tgz
wget http://kaiju.binf.ku.dk/database/kaiju_index_nr_euk.tgz
wget http://kaiju.binf.ku.dk/database/kaiju_index_nr_euk.tgz
2.2 可执行程序

本地编译安装:

git clone  https://github.com/bioinformatics-centre/kaiju
cd kaiju/src
make

直接下载:

wget https://github.com/bioinformatics-centre/kaiju/releases/download/v1.5.0/kaiju-1.5.0-linux-x86_64.tar.gz
tar xzvf  kaiju-1.5.0-linux-x86_64.tar.gz
mv  kaiju-v1.5.0-linux-x86_64-static  kaiju-1.5.0

记得添加环境变量:

export PATH=$PWD/kaiju/bin:$PATH
export PATH=$PWD/kaiju-1.5.0/bin:$PATH

下载完数据库和软件, 我们可以愉快的处理我们的数据了。

3. 如何使用

3.1 构建本地数据库

格式要求(官方实例):

>1_1358
MAQQRRGGFKRRKKVDFIAANKIEVVDYKDTELLKRFISERGKILPRRVTGTSAKNQRKVVNAIKRARVMALLPFVAEDQN
>2_44689
MASTQNIVEEVQKMLDTYDTNKDGEITKAEAVEYFKGKKAFNPERSAIYLFQVYDKDNDGKITIKELAGDIDFDKALKEYK
EKQAKSKQQEAEVEEDIEAFILRHNKDDNTDITKDELIQGFKETGAKDPEKSANFILTEMDTNKDGTITVKELRVYYQKVQ
KLLNPDQ
>3_352472
MKTKSSNNIKKIYYISSILVGIYLCWQIIIQIIFLMDNSIAILEAIGMVVFISVYSLAVAINGWILVGRMKKSSKKAQYED
FYKKMILKSKILLSTIIIVIIVVVVQDIVINFILPQNPQPYVYMIISNFIVGIADSFQMIMVIFVMGELSFKNYFKFKRIE
KQKNHIVIGGSSLNSLPVSLPTVKSNESNESNTISINSENNNSKVSTDDTINNVM
>4_91061
MTNPFENDNYTYKVLKNEEGQYSLWPAFLDVPIGWNVVHKEASRNDCLQYVENNWEDLNPKSNQVGKKILVGKR

序列名称由三部分组成, a. “标识符” b. “_” 下划线,c. “NCBI Taxonomy ID”

执行下面命令:

mkbwt -n 10 -l 1   -a ACDEFGHIKLMNPQRSTVWY -o database  database.faa
mkfmi database
rm database.bwt  database.sa

先后执行了:mkbwt(BWT变换)和mkfmi(FM索引)

3.2 序列分类

这里使用先前下载的标准 kaiju 数据库(Refseq)

命令行接口:

$ kaiju
Error: Please specify the location of the nodes.dmp file, using the -t option.

Kaiju 1.5.0
Copyright 2015,2016 Peter Menzel, Anders Krogh
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>

Usage:
   kaiju -t nodes.dmp -f kaiju_db.fmi -i reads.fastq [-j reads2.fastq]

Mandatory arguments:
   -t FILENAME   Name of nodes.dmp file
   -f FILENAME   Name of database (.fmi) file
   -i FILENAME   Name of input file containing reads in FASTA or FASTQ format

Optional arguments:
   -j FILENAME   Name of second input file for paired-end reads
   -o FILENAME   Name of output file. If not specified, output will be printed to STDOUT
   -z INT        Number of parallel threads (default: 1)
   -a STRING     Run mode, either "mem"  or "greedy" (default: mem)
   -e INT        Number of mismatches allowed in Greedy mode (default: 0)
   -m INT        Minimum match length (default: 11)
   -s INT        Minimum match score in Greedy mode (default: 65)
   -x            Enable SEG low complexity filter
   -p            Input sequences are protein sequences
   -v            Enable verbose output

kaiju 命令行接口还算比较简单,这里我们只需要设置一下几个参数:

-t: NCBI的Taxonomy数据库 nodes.dmp 文件;
-f: Kaiju的索引文件;
-i: 输入文件;
-j: 双端序列模式,第二个序列使用 -j参数;
-z: 线程数;
-a: 选择运行模式, mem 还是 greedy;
-m: 最小匹配长度, MEM模式,一般采用默认参数11;
-e: 贪婪模式最大允许错配个数,官方Web使用5, 默认0个;
-s: 最小匹配Score值,默认65,官方Web使用75(配合-e=5使用);
-p: 如输入为氨基酸序列,使用 -p 可选项;
-x: 低复杂度过滤;
-o: 输出文件;

MEM模式:

kaiju -t  /biostack/database/kaiju/nodes.dmp           \
      -f /biostack/database/kaiju/kaiju_db_nr_euk.fmi  \
      -i  10000_R1.fastq  -j  10000_R2.fastq           \
      -o  10000_mem.tsv -a  mem  -z 24  -m 11

Greedy模式:

kaiju -t  /biostack/database/kaiju/nodes.dmp           \
      -f /biostack/database/kaiju/kaiju_db_nr_euk.fmi  \
      -i  10000_R1.fastq  -j  10000_R2.fastq           \
      -o  10000_greedy.tsv -a greedy -z 24 -e 5 -s 75

3.3 后续处理

输出文件格式

标准的Kraken文件格式:

C       ST-E00144:247:HNJNLCCXX:3:2121:14052:14230      286730

三列解读:

第一列: 'C' 或者 'U', 分类标识符,C为分类的, U为未分类的;
第二列: 序列名称;
第三例: NCBI Taxonomy ID;

物种统计:

$kaijuReport
Error: Please specify the location of the names.dmp file with the -n option.

Kaiju 1.5.0
Copyright 2015,2016 Peter Menzel, Anders Krogh
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>

Usage:
   kaijuReport -t nodes.dmp -n names.dmp -i kaiju.out -o kaiju.report

Mandatory arguments:
   -i FILENAME   Name of input file
   -o FILENAME   Name of output file.
   -t FILENAME   Name of nodes.dmp file
   -n FILENAME   Name of names.dmp file.
   -r STRING     Taxonomic rank, must be one of: phylum, class, order, family, genus, species

Optional arguments:
   -m FLOAT      Number in [0, 100], denoting the minimum required percentage for the taxon \
                 to be reported (default: 0.0)
   -c INT        Integer number > 0, denoting the minimum required number of reads for the \
                 taxon to be reported (default: 0)
   -u            Unclassified reads are not counted for the total reads when calculating \
                 percentages for classified reads.
   -p            Print full taxon path.
   -l            Print taxon path containing only ranks specified by a comma-separated list,
                 for example: superkingdom,phylum,order,class,family,genus,species
   -v            Enable verbose output.

Only one of the options -m and -c may be used at a time.

这里统计序列分类有几点需要注意:

  1. 输出结果为指定分类水平的相对丰度,这里需要注意是分母的问题, 可以包含未分类序列也可以不包含未分类序列, 通过 -u 切换。
  2. 可以将特定分类水平 reads数小于指定值(-c 可选项),或者相对丰度小于特定值(-m可选项),归入单独分类, 这样在后期统计或者可视化分析还是比较有好处。

执行:

kaijuReport -t /biostack/database/kaiju/nodes.dmp \
            -n /biostack/database/kaiju/names.dmp \
            -i 10000_mem.tsv \
            -o 10000_mem.report

这样可以轻松获取样本种群结构的:counts数和相对丰度信息。

Krona可视化:

kaiju2krona
Error: Error: Please specify the name of the output file, using the -o option.

Kaiju 1.5.0
Copyright 2015,2016 Peter Menzel, Anders Krogh
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>

Usage:
   kaiju2krona -t nodes.dmp -n names.dmp -i kaiju.out -o kaiju2krona.out

Mandatory arguments:
   -i FILENAME   Name of input file
   -o FILENAME   Name of output file.
   -t FILENAME   Name of nodes.dmp file
   -n FILENAME   Name of names.dmp file

Optional arguments:
   -v            Enable verbose output.
   -u            Include count for unclassified reads in output.

kaiju2krona的作用就是将Kaiju输出的结果转化成KtImport

执行:

kaiju2krona  -u                           \
   -t /biostack/database/kaiju/nodes.dmp  \
   -n /biostack/database/kaiju/names.dmp  \
   -i 10000_mem.tsv  \
   -o  10000_krona.tsv

ktImportText  10000_krona.tsv -o  10000_krona.html

通过-u设置包含为分类的序列信息。

lineage注释

cat  10000_mem.tsv | grep -Pv "^U" \
     |cut -f2,3  \
     | taxon-translate  -c 2  /biostack/database/taxonomy/ncbi.map  -

这样我们可以获取每条序列的分类路径,比如:

ST-E00144:247:HNJNLCCXX:3:2108:31538:65230      186802  d:Bacteria;p:Firmicutes;c:Clostridia;o:Clostridiales

通过本文介绍可以大致了解Kaiju的算法原理以及如何对序列进行分类和构建本地数据库,剩下的就是实践。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有。

Last Update: 2017-11-22 5:44 PM

Comments are closed.