今天开始会介绍一系列Metagenome（以后使用元基因组）数据分析中的序列分类问题，这类问题一直是研究的热点和难点。
热点：有用，比如直接鉴定环境样本（粪便、痰液、肺积液）病原，算法多样性， 从基于碱基组成统计分类、基于核酸
Kmer查找、基于氨基酸序列MEM查找到序列比对，各种算法八仙过海、各显神通。难点：序列分类准确性、敏感性、内存有效性、计算密集性等都是待解决问题。
Kaiju是一款基于在氨基酸空间进行相似性搜索的序列分类算法。
1. 算法介绍
Kaiju算法发表在了Nature Communications杂志 ：Fast and sensitive taxonomic classification for metagenomics with 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回溯。
几点说明：
- BLOSUM62 矩阵， ORF Score 以及后续的MEM左侧延伸后Score值计算都是按照 BLOSUM62计算。
Kaiju的Greedy模式敏感度要比MEM模式高，但是牺牲了分类速度，对于一些新的基因尤为明显，Greedy采用的模式是启发式算法，只对MEM的末端延伸。- LCA算法见下图，从所有叶子节点向上回溯，找到所有命中节点的共同祖先，比如
a|b|c叶子节点的LCA节点就是x节点。
图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.
这里统计序列分类有几点需要注意：
- 输出结果为指定分类水平的相对丰度，这里需要注意是
分母的问题, 可以包含未分类序列也可以不包含未分类序列， 通过-u切换。- 可以将特定分类水平 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
