Archives

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

一、引言

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

DADA2

DADA2 Species assignment 逻辑很简单:

  1. 预处理16S库文件(去除未分类的序列或者分类不明确的序列, 保留属和种水平分类信息), DADA2 提供了rdp_species_assignment_16.fa.gzsilva_species_assignment_v128.fa.gz 两种库。
  2. 序列精确比对,DADA2使用了ShortRead进行序列比对;
  3. 判断,如果所有有效匹配(可设正向和反向选项)在 属+种 分类都一致,那么是可以直接分配该种分类等级, 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 {
    keep <- as.integer(allowMultiple)
  }
  # Get character vector of sequences
  seqs <- getSequences(seqs)
  # Read in the reference fasta
  refsr <- readFasta(refFasta)
  ids <- as(id(refsr), "character")
  # Crude format check
  if(!length(unlist(strsplit(ids[[1]], "\\s"))) >= 3) {
    if(length(unlist(gregexpr(";", ids[[1]]))) >= 3) {
      stop("Incorrect reference file format for assignSpecies (this looks like a file formatted for assignTaxonomy).")
    } else {
      stop("Incorrect reference file format for assignSpecies.")
    }
  }
  genus <- sapply(strsplit(ids, "\\s"), `[`, 2)
  species <- sapply(strsplit(ids, "\\s"), `[`, 3)
  # Identify the exact hits
  hits <- vector("list", length(seqs))
  lens <- nchar(seqs)
  for(len in unique(lens)) { # Requires all same length sequences
    seqdict <- PDict(seqs[lens==len])
    vhit <- (vcountPDict(seqdict, sread(refsr))>0)
    if(tryRC) vhit <- vhit | (vcountPDict(seqdict, reverseComplement(sread(refsr)))>0)
    hits[lens==len] <- lapply(seq(nrow(vhit)), function(x) vhit[x,])
  }
  # Get genus species return strings
  rval <- cbind(unlist(sapply(hits, mapHits, refs=genus, keep=1)),
                unlist(sapply(hits, mapHits, refs=species, keep=keep)))
  colnames(rval) <- c("Genus", "Species")
  rownames(rval) <- seqs
  if(verbose) cat(sum(!is.na(rval[,"Species"])), "out of", length(seqs), "were assigned to the species level.\n")
  rval
}

该思路可以迁移其它扩增子数据分析平台,比如QIIME2uParse,为了更方面命令行模式下执行该分析,我们重新实现了一个独立版本: species_assign

二、 species_assign程序

2.1 准备数据库

数据库推荐使用silva_species_assignment_v128, 通过下面命令获取序列和分类映射表;

wget https://zenodo.org/record/824551/files/silva_species_assignment_v128.fa.gz
gunzip silva_species_assignment_v128.fa.gz
grep ">" silva_species_assignment_v128.fa | sed 's/>//' | cut -f1,2,3 -d ' ' >silva_species_assignment_v128.map

构建Bowtie索引库:

bowtie-build  example_species_assignment.fa  example_species_assignment
2.2 命令行接口

species_assign 提供了和dada2 assignSpecies 一致的可选项, 即

  1. 正链还是双链匹配模式
  2. 匹配多个物种, 一般不会选择这个参数。

命令行接口:

$ species_assign
Usage: species_assign [options] <map> <bowtie>
Options:
  -b  enable both strand. Default:[plus]
  -m  enable multiple species assignment. Default:[not allowed]
  -v  print version number

可选项:

-b: both strands.
-m: multiple species match.
2.3 程序示例

结合bowtie我们可以轻松完成dada2的物种注释功能:

bowtie  -f -v 0 -a silva_species_assignment_v128 example_seqs.fa 2>/dev/null |\
   species_assign  silva_species_assignment_v128.map -

解读:

bowtie 参数:

-f: 输入格式为fasta;
-v: end-to-end 比对模式, 允许0个错配。
-a:输出所有匹配;
2>/dev/null:将log信息丢掉;

species_assign 程序在 taxon-utils 仓库。

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

Last Update: 2017/11/28 23:41

Comments are closed.