Archives

序列功能注释神器: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盘即可,剩下的和光盘安装一样了。

如何使用Biostack提供的工具

我们分享的工具集都在 Github (地址) 上, 目前只提供编译后的版本, 如下图:

Biostack Github 仓库

如何下载使用工具,目前所有的工具都是在 CentOS 7编译,可能不兼容版本比较低的系统。

1. 下载工具

使用 git 直接克隆对应仓库, 比如:

git clone https://github.com/jameslz/blast_utils

如果提示没有安装 git, 请先安装git:

sudo yum install git 2. 添加环境变量

如果我们将 blast_utils 克隆到了目录: /project/tools/blast_utils, 只需要执行:

export PATH=/project/tools/blast_utils:$PATH

就可以使用这些程序了, 或者执行下面命令添加到bashrc文件,每次启动bash都会自动加载该目录:

echo export PATH=/project/tools/blast_utils:$PATH >> ~/.bashrc source ~/.bashrc

如果提示权限不够,请执行:

chmod -R 775 /project/tools/blast_utils 3. 命令执行

随意切换到任何目录,我们可以直接使用 […]

使用Usearch进行扩增子数据分析:The Full Story

准备工作:

添加辅助程序路径至环境变量:

export PATH=$PATH:/biostack/tools/pipelines/div_seq-0.2.3/bin/utils export PATH=$PATH:/biostack/labs/seqtk_utils/release export PATH=$PATH:/biostack/labs/tabtk_utils/release export PATH=$PATH:/biostack/tools/microbiome/KronaTools-2.7/bin

准备好原始数据,样品信息表、以及引物信息。 软件需要安装:R、seqtk、usearch、seqtk_utils、tabtk_utils、biom以及fasttree、 Trimal, mafft、KronaTools以及phylommand

1. 数据整体质量评估

参考: 获取fasta/q文件的碱基组成和长度信息: seqtk_info 和 seqtk_fqchk

mkdir quality seqtk_fqchk data/raw.R1.fastq.gz >quality/raw_R1.tsv fqchk_base.R quality/raw_R1.tsv quality/raw_base_R1.pdf R1 fqchk_qual.R quality/raw_R1.tsv quality/raw_qual_R1.pdf R1

可以生成质量和碱基分布图:

2. 数据拆分

参考: 快速对扩增子序列进行数据拆分: seqtk_demultiplex

mkdir demultiplex cut -f1,2 mapping_file.txt | sed ‘s/,/\t/’ >demultiplex/barcode.txt seqtk_demultiplex -b demultiplex/barcode.txt -1 […]

扩增子测序分析之去噪 – DADA2及其平台化

引言

DADA2(Divisive Amplicon Denoising Algorithm)的核心 去噪 然后使用 Amplicon Sequence Variants (ASVs)概念构建类OTU表;

具体算法可以参考下面几篇文章:

Denoising PCR-amplified metagenome data DADA2: High-resolution sample inference from Illumina amplicon data Bioconductor Workflow for Microbiome Data Analysis: from raw reads to community analyses

二、核心算法

DADA2 的核心算法:

1. 错误率模型计算,采用selfConsist无监督学习模型; 2. Rate λ_ji: 通过错误率模型,衡量扩增子序列i是否来在来自j; 3. 丰度P-value值,使用丰度信息衡量:该丰度的扩增子序列是否可以通过错误模型解释,使用了一个经验P值来判断,比如 OMEGA_A = 1e-40 4. 分裂算法( divisive partitioning)形成扩增子簇,簇心定义为denoise序列,簇大小定义denoise序列丰度。 […]

扩增子序列样本序列丰度谱: fastx_uniques、fastx_uniques_persample、 DADA2 和 seqtk_seqtab

本文介绍扩增子序列中的一个分析需求,扩增子序列中样本序列丰度谱,即:Uniques 序列在每个样本的Tag分布, DADA2使用这个概念构建类似OTU表的特征表, QIIME2 也在使用这样特征表概念,使用这个策略,一定要有去噪 这个数据处理过程,去除测序错误造成的序列多样性。

本文只从技术层次上介绍如何产生这样的特征表,具体生物学意义可以参考下面两篇文章:

The 97% threshold is much too low, especially for V4 Exact sequence variants should replace operational taxonomic units in marker-gene data analysis 1. fastx_uniques

fastx_uniques 计算fasta/q序列中的Uniq序列,支持以下参数(详细信息见:fastx_uniques command ):

-threads :线程数; -uc :标准输出格式; -sizeout :输出簇大小信息; -relabel :重新标注序列标识符; -fastaout : fasta文件格式输出; -minuniquesize: 最小序列数/per sample 阈值,如果低于这个数值,序列不输出; -topn N : 输出按降序排列的前N条序列; -strand […]