Archives

驯化VirSorter: 预测metagenome contigs的prophage

鉴定噬菌体的工具有Phage_Finder, Prophinder , PHAST,PhiSpy,不过要讲的是 VirSorter。VirSorter适合不完整的基因组,单细胞基因组,宏基因组。

VirSorter 运行时间很长问题,主要问题是 HMMER的问题,HMMER支持多线程不理性,即使设置多线程,实际执行的时候基本都是单线程,导致运行时间比较长。

那解决这个问题的方式就是:将线程强制变成进程,根据hmmsearch的特点,将库文件拆分成指定ncpu份, 单独提交可以达到并行目的, 然后将拆分后的结果合并为输出文件即可。

为解决这个问题 Biostack,实现了 hmmsearch-virsorter 做为HMMSEARCH任务提交的中间件,替换掉VirSorter的提交方式,可以顺利进行真实的并行任务提交。

$ hmmsearch-virsorter Program: hmmsearch-virsorter: HMM based annotation. Version: 0.0.1 Contact: ZHANG LEI <zhanglei@logicinformatics.com> Usage: hmmsearch-virsorter [options] <sequence> <tblout> <output> Options: -c INT CPU number, default: [40] -d STR database location, default: [/biostack/database/pfam/Pfam-A.hmm]

现在利用40线程,一个典型的细菌基因组基本3分钟就可以完成前噬菌体鉴定。

生物信息数据工具封装:序列比对之hmmscan-pipe

一、前言

为什么要封装这个hmmscan-pipe呢, hmmscan为HMMER(当前最新版本 3.1b2)程序包的子程序,工作模式是: 蛋白序列对谱序列库(hmmpress构建索引), 常用的氨基酸序列功能注释工具, 常用功能信息数据库有: Pfam、 Superfamily、 dbCAN、 SMART、 TIGRFAM 等。

hmmscan 当前版本可使用 –cpu 指定使用的线程数,但是一般不能有效利用多核心资源,所以最佳实践为: 序列拆分。 使用 fastx-utils partition 分割成指定线程数文件,使用进程对每个文件单独提交,这样就比较适合在集群模式下工作。 hmmscan-pipe依赖hmmscan-utils 支持标准输入流,主要解决两个问题:

hmmscan-utils domtblout, 格式化hmmscan输出格式,使用制表符分隔,采用默认的过滤模式:如果比对片段 >80aa Evalue 阈值使用1e-5, <80aa Evalue 阈值使用1e-4。 hmmscan-utils resolve, 解析HMM匹配区域的交叠问题,去除交叠区域比较大的比对, 实现了文章 A fast and automated solution for accurately resolving protein domain architectures BMC 算法。

hmmscan-utils 程序提供了两子命令程序:

Usage: hmmscan-utils <command> […]

毒力因子注释Protocol:VFDB数据库

一、毒力因子

毒力因子(Virulence factor), 详细介绍参见 维基百科 Virulence_factor 页面, 细菌、病毒、真菌等生成的分子,并产生毒力(主要有侵袭力和毒素等),包括:

1. 在宿主定殖 (colonization),黏附在宿主消化道、呼吸道、生殖道、尿道及眼结膜等处,以免被肠蠕动、黏液分泌、呼吸道纤毛运动等作用所清除 2. 免疫逃避,逃避宿主的免疫应答 3. 免疫抑制,抑制宿主的免疫反应 4. 进入和退出细胞 5. 从宿主获得营养

毒力因子可编码在可移动遗传元件(比如质粒、基因岛、噬菌体等)上并进行水平基因转移(传播),使无害细菌变成危险的病原菌,所以在鉴定毒力因子时一般会考虑: 基因岛、分泌蛋白等。

二、病原菌毒力因子数据库 VFDB

毒力因子数据库VFDB 由中国医学科学院研发,被广泛应用于毒力因子基因鉴定。 VFDB收集了包括30个属( 74个病原菌)的细菌毒力基因序列信息。

VFDB提供了对应的毒力基因核酸和蛋白质序列信息,因此鉴定毒力基因最简单的办法就是序列比对(BLAST),

2.1 数据库预处理

数据预处理需要以下几个步骤:

VFDB的元信息可以通过序列文件以及提供的描述文件获得:

>VFG000676(gb|AAD32411) (lef) anthrax toxin lethal factor precursor [Anthrax toxin (VF0142)] [Bacillus anthracis str. Sterne]`

1、格式化序列文件,只保留毒力基因编号 VFG000676 并获得 VFG000676 -> VF0142 映射关系; […]

生物信息数据工具封装:序列比对之blast-pipe

一、前言

为了生物学家更加容易使用命令行模式的生物信息工具,在数据分析流程集成工具水平下我们设计了 “小封装” 模式,即封装一些常用的几个工作步骤,尽量使用优化后的默认参数,整个小封装依赖实现的 “tool-utils”, 比如 blast-utils、fastx-utils 、 tsv-utils、 sam-utils, 这些小程序一般都是 C语言 实现的高性能应用。

二、blast-pipe介绍

小封装 blast-pipe 目的就是解决序列比对BLAST的提交任务, 命令行接口:

$blast-pipe Program: blast-pipe: blast submit and parse protocol. Version: 0.0.1 Contact: ZHANG LEI <zhanglei@logicinformatics.com> Usage: blast-pipe [options] <sequence> <project> Options: -t STR blast type. blastx|blastp|blastn, default [blastp], for special task, can do like this: ‘blastn -task megablast’ […]

序列功能注释神器:eggNOG-mapper,KEGG/COG/KOG/GO/BiGG 一网打尽

一、引言

对于测序基因组进行 KEGG(Kyoto Encyclopedia of Genes and Genomes)和COG(clusters of orthologous groups,对直系同源基因进行聚类)功能注释,基本成为基因组注释的标配内容, 特别是微生物基因组基因注释,其基因功能注释逻辑基础是 直系同源基因具有相同的功能,最经典的鉴定直系同源基因策略是 BBH(bi-directional best hit)策略,但是通常最直接的直系同源基因很难鉴定,而对同源基因进行聚类并定义一个簇会是更好的策略:每一个簇会包含直系同源基因(伴随物种形成事件出现)和旁系同源基因(伴随拷贝事件出现),每一簇共享同一个功能, KO(KEGG Orholog), COG, eggNOG 等都是基于聚类的方式定义簇,并对簇进行注释。

今天要讲的是eggNOG, eggNOG的出现要从COG说开,下面看看NCBI COG的数据库主要更新历史:

从 1997 年 第一个公布版本,7个完整基因组,720个COG分类, 包含原核基因组和单细胞真核基因组(酵母),2003 年和2014 年进行了版本升级,最后只保留了细菌和古菌,包含了711个基因组以及4,631个COG分类, 26个功能分类。 2013 年构建真核分支COG(KOG, Eukaryotic orthologous groups); 2007 年构建古菌分支COG(arCOG, Archaeal Clusters of Orthologous Genes),2012 年和2014 年arCOG进一步升级,arCOG比较适合用于古菌基因组注释; 2011 年构建Phage分支COG(POG,phage orthologous groups),2013 年进行了升级;

由于计算资源需求,NCBI COG 构建了不同系统分类分支的COG簇,比如arCOG,KOG, POG等,推荐使用这些分支对新测序基因组进行注释,其实eggNOG […]

otutable-utils之计算OTU表汇总信息-otutable_summary和otutab_stats

本文介绍otutable-utils系列的一个小工具:otutab_summary (Github 链接 ), 就是计算OTU表中每个样本的OTU数和tag数,支持未注释的OTU表。

1. 命令行接口 $ otutable_summary Usage: otutable_summary <otutab> version: 0.0.1 2. 命令行接口 otutable_summary OTU_table.txt >OTU_table.summary.txt 3. 其它实现

3.1 otutab_stats

otutab_stats 是 Usearch 的命令行工具,接口也比较简单

usearch -otutab_stats OTU_table.txt -output OTU_table.report.txt

该程序对整个OTU表中样本进行统计分析, 示例结果如下:

1624812 Reads (1.6M) 27 Samples 3565 OTUs 96255 Counts 44360 Count =0 (46.1%) 14975 Count =1 (15.6%) 14774 Count >=10 (15.3%) […]

扩增子序列分析之物种分类: 注释到种水平,也许可以!

一、引言

DADA2 升级至 1.6 版本, 并正式支持 Species assignment, 使用了精确匹配模式判别是否可以分类到种水平,当然种以上分类水平还是要使用RDP算法或者Sintax算法。

DADA2 Species assignment 逻辑很简单:

预处理16S库文件(去除未分类的序列或者分类不明确的序列, 保留属和种水平分类信息), DADA2 提供了rdp_species_assignment_16.fa.gz和silva_species_assignment_v128.fa.gz 两种库。 序列精确比对,DADA2使用了ShortRead进行序列比对; 判断,如果所有有效匹配(可设正向和反向选项)在 属+种 分类都一致,那么是可以直接分配该种分类等级, DADA2 提供了允许分配多个种的选项,不推荐使用这个模式。

代码段:

assignSpecies <- function(seqs, refFasta, allowMultiple=FALSE, tryRC=FALSE, verbose=FALSE) { # Define number of multiple species to return if(is.logical(allowMultiple)) { if(allowMultiple) keep <- Inf else keep <- 1 } else { […]

Linux 命令行小工具之tldr:更好的 man?

一、引言

Linux的伟大成功一个重要的原因就是有各种各样的小工具,比如ls, du、top、du等,可以完成Linux操作系统的各种文件操作、目录操作以及系统管理等,但是,以下几个因素也会成为我们习惯命令行操作的拦路虎:

可是使用的命令行工具太多了。 命令行接口(选项)太多了。

当安装这些小工具不是问题的时候,使用这些小工具就成为最大的问题了,首当其冲就是了解命令行的使用说明,最直接的肯定是 man 工具, 如下图,我们可以在命令行下输入man man就可以其实获得关于man的使用说明。

图一: man使用说明

当然很多程序, 可以通过直接执行命令行工具本身, 比如samtools 或者 bedtools –help 获取帮助信息。

不过今天介绍是一个小工具tldr()。

二、TLDR介绍

tldr-pages 采用一种社区驱动模式, Github仓库地址: https://github.com/tldr-pages/tldr, Web 客户端: https://tldr.ostera.io/。

Linux命令行工具运行命令如下:

图二: tldr使用模式

相对于man页面,tldr具有更多的例子。

三、安装

tldr 提供了各种语言的客户端: 下面列出几种比较常用的安装模式:

Node.js 客户端(需要先安装 Node.js) :

npm install -g tldr

Python 客户端

pip install tldr

GO 客户端:

wget https://github.com/pranavraja/tldr/releases/download/v1/tldr_linux_amd64.tar.gz […]

庖丁解牛,微生物基因组重测序:一条命令完成从sra到vcf

一、引子

前面在 Twitter上看到一个针对细菌基因组重测序流程(鉴定SNP)的使用工具调查,这里我们就庖丁解牛般解读下这些流程工具是如何工作的,流程又是如何设计的。

该调查数据来自 google.docs

从图中可以看出snippy还是比较受欢迎的,Torsten Seemann的流程一向比较稳定,接口一向比较简单。

这次解读的主角是snippy,使用起来确实很方便,比如:

fastq-dump –defline-seq ‘@$sn/$ri’ –defline-qual ‘+’ ERR114970.sra -Z –split-3 |\ seqtk trimfq – | \ snippy –cpus 16 –outdir ERR114970 –ref GCF_000069185.1_ASM6918v1_genomic.gbff \ –prefix ERR114970 -peil /dev/stdin &>log.txt

一条命令可以解决我们从sra到vcf到的问题,当然其中使用的工具很多了,但是都被封装在脚本里了。

二、流程细节 2.1. sra生成interleave格式fastq文件 fastq-dump –defline-seq ‘@$sn/$ri’ –defline-qual ‘+’ ERR114970.sra -Z –split-3

该结果使用-Z 选项将结果输出至标准输出流,格式是Interleave格式,这种模式有很多好处, 比如质量值修剪,序列比对等,后续很多程序都支持,这样比较方便以steaming方式工作。

内部链接: ENA和SRA数据预处理: SRA Toolkit […]

使用U盘安装CentOS7 Linux 操作系统

现在很多服务器都已经放弃了光驱设备,使用U盘安装操作系统成为比较方便的方案,本文介绍如何烧制可引导的U盘系统。

1. 下载ISO镜像

CentOS提供的所有的镜像列表: List of CentOS Mirrors 我们选择:Sohu Inc, Beijing P.R. China

推荐下载 DVD版本: CentOS-7-x86_64-DVD-1708.iso

2. 准备U盘

准备一个8Gb大小的U盘。

3. 安装

安装 win32diskimager-1.0.0, Windows平台下如何安装不需多说。

4. 烧制

打开程序,界面很简单,插入U盘,选中U盘符。

然后直接点击 Writer 确认后完成烧制, USB 2.0的接口估计需要15分钟左右才能烧制完成。

注意: 确保U盘里没有重要文件。

5. 系统安装

一般需要设置下, 让sBIOS系统引导的优先级,选中U盘即可,剩下的和光盘安装一样了。