今天开始会介绍一系列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