NCBI Taxonomy 数据库更新，提供lineage、host信息

我们分析metagenome数据离不开使用NCBI的Taxonomy数据，NCBI Taxonomy 提供了一棵物种树，其实每个节点（Node）都分配了一个数字标识符，可以唯一描述一个系统分类信息。

NCBI Taxonomy 数据库提供了一个 taxdump.tar.gz, 并记录了节点的描述信息（names.dmp）以及树的上下游信息（nodes.dmp），刚刚发布的更新版本提供了额外的lineage信息（rankedlineage.dmp）以及 host 信息。

另外NCBI已经不再给Strain水平分配这种数字标识符，所以NCBI Taxonomy 提供了 typematerial.dmp 文件用于关联种和菌株(strain)的映射关系。

利用新的数据库我们可以很容易对一些短序列分类器进行注释，常用的操作如下：

1、格式化数据库，一般可以使用 tsv-utils

cut -f1,5 fullnamelineage.dmp | sed ‘s/ $//’ >fullnamelineage.db cut -f1,5 taxidlineage.dmp | sed ‘s/ $//’ >taxidlineage.db cut -f1,3 host.dmp >host.db

2、典型使用场景

下面以Kraken为例子，介绍如何格式化为有效信息， kraken的结果：

C E00552:27:HJ2JYALXX:4:1101:5233:1801 435590 203 816:40 435590:21 A:31 435590:13 0:53 […]

taxon-utils之NCBI Taxonomy 标识符lineage注释: taxon-translation

一、引言：

本文介绍一个小工具，taxon-translation, 使用场景:

需要将NCBI Taxonomy数据库的数字标识符转换成整个Lineage信息。

taxon-utils 仓库地址： https://github.com/jameslz/taxon-utils

二：应用 2.1 数据库格式化 taxon-db

taxon-db程序将NCBI的Taxonomy数据库的 node.dmp 和 name.dmp 合并成一个文件。

执行方法：

taxon-db nodes.dmp names.dmp >ncbi.map gzip ncbi.map

ncbi.map 也可以从 taxon-utils 仓库下载， ncbi.map.gz（版本为2017-09-15）链接。

2.1 数据库格式化 taxon-translate

taxon-translate 命令行接口比较简单：

$ taxon-translate Usage: taxon-translate [options] <ncbi-map> <tab> Options: -c INT target col for taxon translate, default: 1 -v print version number […]

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

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

[…]

CentOS 7运维-Windows平台如何使用SSH登录Linux服务器

本文介绍Windows平台下两款登录Linux终端的应用。

一、Putty: 最经典知名的免费SSH客户端

下载地址： 64位 / 32位

经过简单配置后可以工作了。

通过1-7步骤，可以登录到远程服务器端，并使用SSH在终端完成各种操作：下面解读下七个步骤：

1.设置远程服务器IP地址，确保可以访问; 2.设置远程服务器IP端口，默认22， SSH访问; 3.设置会话名称; 4.保存下次可以直接会话列表进行选择； 5.打开终端，当然可以直接双击会话列表中记录登录； 6.输入用户名 7.输入密码

2、MobaXterm: 全能的远程终端登录软件

MobaXterm功能强大、界面漂亮，可以满足大部分的Windows和Linux之间的交互需求，比如SSH登录，SFTP访问等，当然还有一些好玩的小游戏。

2.1 工具下载

选择Home免费版本 MobaXterm Home Edition Windows 平台如何安装跳过。

2.2 推荐理由

构建本地开发平台可以很方便安装 GCC以及各种库文件，方便Windows开发编译各种软件。内置 X server, 可以执行图形界面应用。多标签，分屏 SFTP文件传输直接支持UTF-8编码，不乱码保留用户名和密码一键操作所有活跃终端

2.3 SSH远程登录

通过1-5步骤，可以登录到远程服务器端，并使用SSH在终端完成各种操作：下面解读下五个步骤：

1.设置新会话； 2.选择SSH登录； 3.设置远程服务器IP地址，确保可以访问; 4.设置远程服务器用户名； 5.设置远程服务器IP端口，默认22；

[…]

blast-utils之表格文件处理：blast-utils小工具集合

BLAST类序列比对工具在生物学数据分析中应用最为广泛，比如16S全长序列比对(MEGABLAST)，氨基酸功能注释(BLASTP)等，另外还有各种快速比对工具，比如： DIAMOND、GHOSTX、RAPSearch2等，这些工具默认都支持制表符分隔的输出文件（BLAST程序指定的 -outfmt 输出格式）。

格式如下：

列中文描述英文描述 1查询序列表标识符Query 2目标序列标识符Target 3相似度identity 4比对长度Alignment length 5错配数Number of mismatches 6gap数Number of gap opens 7查询序列比对起始位置Start position in query 8查询序列比对终止位置End position in query 9目标序列比对起始位置Start position in target 10目标序列比对终止位置End position in target 11E值E-value 12Bit score值Bit score

BLAST比对结果通常一条查询序列可以分几个片段和目标序列比对上，其中每个片段的比对结果都是一个HSP, 比对到的目标序列称之为一个hit，比对结果可以包含多个hit, 我们对序列注释有时候采用的策略是 best hit 或者 best_hsp, 这样就需要对比对结果进行各种处理。

blast_utils （github 仓库）小工具集合，主要使用场景包括：

根据E值、Bit […]

16S扩增子组成谱预测功能谱：本地快速执行otutable_norm和otutable_melt以及otutable_kann

当前可以使用16S扩增子组成谱预测功能谱的工具主要有：PICRUSt、Tax4fun 和 Piphilin，本文主要使用PICRUSt的预测文件，进行重新实现。

PICRUSt 采用了一种ASR（ancestral state reconstruction，可参考1985年的一篇文章：Phylogenies and the Comparative Method）方法重构物种的16S拷贝数和基因组的基因功能信息(使用KEGG、COG描述)，基于系统进化的逻辑是同一个clade分支下的基因组（物种）具有相似的形状，这样就可以使用已测序基因组预测未测序基因组的功能空间，主要依赖参数为系统进化树分支距离。

图来自PICRUSt文章

原理图分拆解读：

1. 使用Greengene 16S序列构建系统进化树，去除未测序基因组； 2. 预测所有祖先状态的性状（16S拷贝数、基因功能）； 3. 预测未知基因组节点的性状（16S拷贝数、基因功能）； 4. 标准化，使用预测的拷贝数对OTU丰度进行标准化； 5. 基因功能空间预测（KO、COG）丰度表；

其实关键算法在1、2、3部分。这篇文章我们直接使用PICRUSt根据 Greengene 13.5 进行的预测结果进行一些数据分析。

需要的文件为双端合并后的序列数据（Usearch sample=’xxx’ 注释文件），我们会在下面一步一步演示如何操作。准备工作:

下载PICRUSt预计算拷贝数和KO表: 16S_13_5_precalculated.tab.gz和 ko_13_5_precalculated.tab.gz 下载Greengene 97OTU代表序列: 97_otus.fasta.gz

如何使用这些程序，请移步：如何使用biostack 提供的命令行工具

1. 构建OTU表

虽然有文章讨论过使用closed- 和 open-reference 构建 OTUs 有多不准确，如下：

Accuracy of […]

fasta/q序列抽取和排序: seqtk_reoder

一、使用场景

给定fasta/q文件和一个序列列表，按照给定的序列顺序抽取序列，如不存在则忽略，先前的 seqtk_subseq（如何根据列表快速提取fasta/q序列子集/补集:seqtk_subseq ）抽取出序列按照标识符在序列文件中出现的顺序输出，有时我们需要按照给定的标识符输出序列。

seqtk_reoder使用reoder作为名字输出就是指序列输出顺序按照指定顺序输出。

二、工具介绍

通过seqtk_info 查看测试序列信息：

$seqtk_info shear.fastq #sequence base min_len max_len avg_len 47525 23905075 503 503 503.00

一共47525序列；

2.1 seqtk_split

命令行接口：

$ seqtk_reorder Usage: seqtk_reorder <fasta/q> <list> version: 0.0.1

seqtk_reorder 接受两个参数，序列文件和列表文件，支持stream流模式；

实例

$cat shear.idx | seqtk_reorder shear.fasta – $cat shear.idx | seqtk_reorder shear.fastq.gz – $cat shear.idx | seqtk_reorder […]

fasta/q序列拆分三剑客: seqtk_split、seqtk_partition 和 seqtk_shred

一、使用场景

本部分内容主要解决以下几个问题：

对指定fasta/q序列文件拆分成每份N条序列问题，比如执行InterProScan。对指定fasta/q序列文件拆分成指定N份问题，比如我们有40核心，程序不支持多线程（或者支持的不理想），那拆分40份多进程执行，比如HMMER 对指定fasta/q序列文件序列拆分，每条序列保留最大l长度，并保证分割非序列存在N个碱基交叠。使用的场景多，比如 454 Newbler拼装软件，不支持很长度序列，如果我们对不同拼装版本的组装结果混合组装，将序列撕碎（保留比较大的overlap）再进行组装，就没有问题了。二、工具介绍

以下三个程序分别解决以上三个问题。

通过seqtk_info 查看测试序列信息：

$seqtk_info shear.fastq #sequence base min_len max_len avg_len 47525 23905075 503 503 503.00

一共47525序列；

2.1 seqtk_split

命令行接口：

$ seqtk_split Usage: seqtk_split [options] <fasta/q> <size> <prefix> Options: -t fastx file suffix: default: [fasta]; -v print version number

seqtk_split 接受三个参数，序列文件、序列大小（条数）和输出文件前缀。可以通过 -t 指定输入文件类型，比如 […]

去除嵌合体： uchime_ref和uchime3_denovo

嵌合体从序列组成来说是指一条序列由两条或者更多条生物学序列连接而成(见下面示意图)，扩增子测序实验中，主要由于PCR阶段不完全延伸形成，形成的嵌合体会在随后PCR过程被放大。具体解释可参考 chimera_formation。

嵌合体示意图:

下面讨论如何去除这些嵌合体，以及什么场景采用什么策略。

USEARCH系统去除嵌合体策略主要有：

使用嵌合体参考库参考：uchime2_ref 从头去除嵌合体参考：uchime3_denovo OTU构建集成: 参考： uparse 和 unoise3

请参考: UCHIME2 算法： UCHIME2: Improved chimera detection for amplicon sequences

1. 使用嵌合体参考库

命令行接口：

usearch -uchime_ref reads.fasta -db silva.udb -uchimeout out.txt -strand plus -threads 20 -mode balanced

主要参数解释：

-db: 参考库，推荐SILVA； -strand：目前只支持 plus，所以序列最好是调整好方向的； -uchimeout：嵌合体输出文件； -threads：线程数； -mode：可以使用的模型，支持 high_confidence、specific、balanced、sensitive等几种模式。\ […]

切除末端和Mask掉氨基酸序列中的终止密码： seqtk_strip_stop_codon

本文介绍一个小工具 seqtk_strip_stop_codon，实现该下程序的起因是 InterProScan 程序，使用InterProScan对氨基酸序列进行功能注释（也是 Gene Ontology 注释的一套方案），对输入序列要求比较严格，不能存在终止密码 “*” 或者 “.”，可通过切除序列末端和Mask(将 “*” 或者 “.” 替换成 “X”)序列中间的终止密码子进行格式化。

程序接口：

$ seqtk_strip_stop_codon Usage: seqtk_strip_stop_codon <fasta> version: 0.0.1

实例

seqtk_strip_stop_codon Trinity.pep

Transdecode 预测的ORF序列末端包含终止密码，如果要使用InterProScan需要使用seqtk_strip_stop_codon 进行处理。

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

Last […]

March 2021
M	T	W	T	F	S	S
« Mar
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30	31