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的数据库主要更新历史:

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

由于计算资源需求,NCBI COG 构建了不同系统分类分支的COG簇,比如arCOG,KOG, POG等,推荐使用这些分支对新测序基因组进行注释,其实eggNOG 尤其是4.x版本也使用了clade特异的聚类模式。

从2003年至2014年NCBI COG一直未更新,EMBL EggNOG(evolutionary genealogy of genes: Non-supervised Orthologous Groups)继承了NCBI COG的衣钵,极大的扩展了基因组信息。当前版本为 4.5.1 版本, 把包含了2,031个基因组, 其中 352病毒基因组, 190k个直系同源家族。

eggNOG 数据库包含了丰富的注释信息,除了COG/KOG/NOG的分类和注释信息外,还包含了KEGG/GO/SMART/PFAM信息。

最新版本的EggNOG 还提供了自动化注释工具eggnog-mapper,可很方便的完成基因组的功能注释,注释信息可以关联COG/KOG/KEGG/GO/BiGG等。

eggNOG在不同的系统分类水平都进行了构建直系同源簇,因此当我们可以确定基因组系统分类时,可以选择系统进化上最近的分支进行注释,以减少计算量, 比如注释的真菌组可以选择fuNOG, 古菌基因组可选择arNOG, 每一个分类水平都提供了HMM谱, 可使用HMMER快速比对。
最新版本eggNOG构建流程如下图:

eggNOG

eggNOG 大致通过以下几个步骤完成:

  1. 选择基因组
  2. all-against-all SW序列比对,形成SIMAP库
  3. 物种分类水平选择
  4. 对已选择分类水平构建直系同源簇,图无监督聚类算法
  5. 根据物种系统分类,建立直系同源簇的等级关系, 根据算法对OG进行合并和分拆使得不同分类水平的OG保持一致。
  6. 对直系同源簇进行功能注释, 包括了KEG/GO/COG/SMART等
  7. 对直系同源簇构建系统进化关系,确立直系同源和旁系同源关系, 并构建HMM谱

了解eggNOG数据库的构建思路有助于我们理解eggnog-mapper是如何通过实现的。

下面我们着重介绍eggnog-mapper 工具, 看看如何使用eggNOG资源实现序列(核酸和氨基酸都可以)的功能注释。

二、eggNOG-mapper介绍

2.1 eggNOG-mapper算法:

image.png-462.8kB

图一、eggNOG-mapper工作流

通过A/B/C/D 四个步骤完成eggNOG注释和功能关联。

下面组条解读:

1. 数据库搜索,可以通过HMMER和DIAMOND获取种子序列,DIAMOND直接与序列库进行相似性搜索,通过 HMMER 与特定 level 的HMMER谱搜索获取最佳HMM,然后使用 phmmer 获取最佳种子序列。

2. 解析直系同源基因 Orthologs, 存在几种模式:
              a. one-to-one, 没有发生duplication 事件;
              b. one-to-many, 物种形成后发生的duplication 事件;
              c. many-to-many;

3. 去除远源Orthologs,选择系统分类上比较近的Orthologs
4. 对Orholog类型进行处理,功能转移(注释)
2.2 eggNOG-mapper数据准备:

当前版本: 1.0.3

下载库文件:

emapperdb-4.5.1

hmmdb_levels/                                      29-May-2016 08:05       -
OG_fasta.tar.gz                                    13-Oct-2017 09:18      6G
eggnog.db.gz                                       17-Oct-2017 12:00      4G
eggnog_proteins.dmnd.gz                            13-Oct-2017 09:18      2G
og2level.tsv.gz   

注: 如果使用的是最新的Diamond版本,需要重新对序列库文件进行格式化,按照下面方式进行操作:

gunzip eggnog_proteins.dmnd
diamond getseq --db eggnog_proteins.dmnd  >eggnog_proteins.pep
diamond makedb --db eggnog_proteins.dmnd  --in  eggnog_proteins.pep

更多使用帮助请参考: Basic-Usage

2.3 安装和使用:

安装比较简单,直接下载解压后即可使用, 可参考 InstallationBasic Usage

git clone https://github.com/jhcepas/eggnog-mapper
export PATH=$PWD/eggnog-mapper:$PATH
emapper.py -i test/polb.fa --output polb_bact -d bact
emapper.py -i test/p53.fa --output p53_maNOG  -m diamond

emapper 的主要命令行接口:

$ emapper.py  -h
usage: emapper.py [-h] [--guessdb] [--database] [--dbtype {hmmdb,seqdb}]
                  [--data_dir] [--qtype {hmm,seq}] [--tax_scope]
                  [--target_orthologs {one2one,many2one,one2many,many2many,all}]
                  [--excluded_taxa]
                  [--go_evidence {experimental,non-electronic}]
                  [--hmm_maxhits] [--hmm_evalue] [--hmm_score]
                  [--hmm_maxseqlen] [--hmm_qcov] [--Z] [--dmnd_db DMND_DB]
                  [--matrix {BLOSUM62,BLOSUM90,BLOSUM80,BLOSUM50,BLOSUM45,PAM250,PAM70,PAM30}]
                  [--gapopen GAPOPEN] [--gapextend GAPEXTEND]
                  [--seed_ortholog_evalue] [--seed_ortholog_score] [--output]
                  [--resume] [--override] [--no_refine] [--no_annot]
                  [--no_search] [--report_orthologs] [--scratch_dir]
                  [--output_dir] [--temp_dir] [--no_file_comments]
                  [--keep_mapping_files] [-m {hmmer,diamond}] [-i]
                  [--translate] [--servermode] [--usemem] [--cpu]
                  [--annotate_hits_table] [--version]

optional arguments:
  -h, --help            show this help message and exit
  --version

Target HMM Database Options:
  --guessdb             guess eggnog db based on the provided taxid
  --database , -d       specify the target database for sequence searches.
                        Choose among: euk,bact,arch, host:port, or a local
                        hmmpressed database
  --dbtype {hmmdb,seqdb}
  --data_dir            Directory to use for DATA_PATH.
  --qtype {hmm,seq}

Annotation Options:
  --tax_scope           Fix the taxonomic scope used for annotation, so only
                        orthologs from a particular clade are used for
                        functional transfer. By default, this is automatically
                        adjusted for every query sequence.
  --target_orthologs {one2one,many2one,one2many,many2many,all}
                        defines what type of orthologs should be used for
                        functional transfer
  --excluded_taxa       (for debugging and benchmark purposes)
  --go_evidence {experimental,non-electronic}
                        Defines what type of GO terms should be used for
                        annotation:experimental = Use only terms inferred from
                        experimental evidencenon-electronic = Use only non-
                        electronically curated terms

HMM search_options:
  --hmm_maxhits         Max number of hits to report. Default=1
  --hmm_evalue          E-value threshold. Default=0.001
  --hmm_score           Bit score threshold. Default=20
  --hmm_maxseqlen       Ignore query sequences larger than `maxseqlen`.
                        Default=5000
  --hmm_qcov            min query coverage (from 0 to 1). Default=(disabled)
  --Z                   Fixed database size used in phmmer/hmmscan (allows
                        comparing e-values among databases).
                        Default=40,000,000

diamond search_options:
  --dmnd_db DMND_DB     Path to DIAMOND-compatible database
  --matrix {BLOSUM62,BLOSUM90,BLOSUM80,BLOSUM50,BLOSUM45,PAM250,PAM70,PAM30}
                        Scoring matrix
  --gapopen GAPOPEN     Gap open penalty
  --gapextend GAPEXTEND
                        Gap extend penalty

Seed ortholog search option:
  --seed_ortholog_evalue
                        Min E-value expected when searching for seed eggNOG
                        ortholog. Applies to phmmer/diamond searches. Queries
                        not having a significant seed orthologs will not be
                        annotated. Default=0.001
  --seed_ortholog_score
                        Min bit score expected when searching for seed eggNOG
                        ortholog. Applies to phmmer/diamond searches. Queries
                        not having a significant seed orthologs will not be
                        annotated. Default=60

Output options:
  --output , -o         base name for output files
  --resume              Resumes a previous execution skipping reported hits in
                        the output file.
  --override            Overwrites output files if they exist.
  --no_refine           Skip hit refinement, reporting only HMM hits.
  --no_annot            Skip functional annotation, reporting only hits
  --no_search           Skip HMM search mapping. Use existing hits file
  --report_orthologs    The list of orthologs used for functional transferred
                        are dumped into a separate file
  --scratch_dir         Write output files in a temporary scratch dir, move
                        them to final the final output dir when finished.
                        Speed up large computations using network file
                        systems.
  --output_dir          Where output files should be written
  --temp_dir            Where temporary files are created. Better if this is a
                        local disk.
  --no_file_comments    No header lines nor stats are included in the output
                        files
  --keep_mapping_files  Do not delete temporary mapping files used for
                        annotation (i.e. HMMER and DIAMOND search outputs)

Execution options:
  -m {hmmer,diamond}    Default:hmmer
  -i                    Input FASTA file containing query sequences
  --translate           Assume sequences are genes instead of proteins
  --servermode          Loads target database in memory and keeps running in
                        server mode, so another instance of eggnog-mapper can
                        connect to this sever. Auto turns on the --usemem flag
  --usemem              If a local hmmpressed database is provided as target
                        using --db, this flag will allocate the whole database
                        in memory using hmmpgmd. Database will be unloaded
                        after execution.
  --cpu
  --annotate_hits_table
                        Annotatate TSV formatted table of query->hits. 4
                        fields required: query, hit, evalue, score. Implies
                        --no_search and --no_refine.

主要可选参数解读:

-i: 输入文件
-o: 输出文件前缀
-m: 使用HMMER策略还是DIAMOND策略,默认使用HMMER,如果eggNOG存在很近的序列,HMMER\
    鉴定的种子子直系同源序列的准确性会较低,选用DIAMOND比较普适,尤其是metagenome \
    序列
--cpu:使用的线程数
--translate:如使用的核酸序列,选择HMMER策略时需要先翻译成氨基酸序列,DIAMOND模式 \
             会自动选择blastx序列比对模式,如果输入的是微生物基因组核酸序列,推荐使\
             用Prodigal预测基因,如果是拼装的转录本,可以考虑使用TRANSDECODE预测 \
             ORF

--usemem: 将 emapper.db 读入内存,他然后执行emapper, 默认为磁盘数据库查找
--servermode:使用服务器模式,即将emapper.db调入内存,供其它进程使用,分配监听端口。
--output_dir:输出文件目录
--report_orthologs: 列出所有进行功能转移的直系同源基因
--no_refine: 只针对HMMER模式,需要结合 --no_annot 使用 
--no_annot: 只汇总鉴定的最佳seed序列以及相应的E值和Bitscore值
--no_search:可直接基于--no_annot的结果进行后续功能注释
--target_orthologs: one2one,many2one,one2many,many2many,all可选。
--data_dir: 数据库目录
--tax_scope: 指定选择的直系同源基因的物种分类范围,默认为自动判断。

从实现上看,emapper使用SQLite数据库存储注释信息,序列比对(HMMER和Diamond)和序列注释是独立的两步,可以分开执行,注释部分可直接使用序列比对的中间结果。

了解这些参数,执行起来就比较简单,比如:

HMMER模式

emapper.py -i test/polb.fa --output polb_bact -d bact

diamond模式

emapper.py -i test/p53.fa --output p53_maNOG  -m diamond

Server模式:

python emapper.py -d bact --cpu 10 --servermode
python emapper.py -d bact:localhost:51500 -i test/polb.fa -o polb

2.4 结果解读:

针对生成文件: polb.emapper.annotations

结果解读:

第一列:查询序列名称;
第二列:eggNOG种子序列;
第三列:eggNOG种子序列 evalue;
第四列:eggNOG种子序列 bit score;
第五列:预测基因名称;
第六列:GO_terms, 预测的GO,分号分隔;
第七列:KEGG_KO: 预测的KO,分号分隔;
第八列:BiGG_Reactions: BiGG代谢反应预测,分号分隔;
第九列:eggNOG Taxonomic Scope信息;
第十列:匹配的OGs;
第十一列:best_OG|evalue|score: Best matching Orthologous Groups (only in HMM mode)
第十二列:COG功能分类;
第十三列:eggNOG功能描述;

注:种子可通过两种方式获取,一种基于HMMER HMM谱 相似性搜索,一种基于DIAMOND序列相似性搜索。
针对Metagenome, 首选使用DIAMOND于整个eggNOG库进行相似性搜索。

细节可参考:Results Interpretation

三、功能关联

程序见: https://github.com/jameslz/ann-utils

GO/KO/BiGG注释

我们根据上述表中的第六列/第七列/第八列 实现了对 KEGG/BiGG/GO 通用的注释提取程序,用于生成类似 goslim 输出结果。

$ emapper-ann
Usage: emapper-ann <text>
Options:
  -d STR   parse annotation, [KEGG, BiGG, GO], default: [KEGG]
  -v       print version number

示例:

emapper-ann -d KEGG   nr.emapper.annotations

KOG/COG注释

有时候我们只关注COG或者KOG分类,这样我们就可以根据 第九列选取对应的KOG/COG匹配信息。

$ emapper-eggnog
Usage: emapper-eggnog [options] <emapper>
Options:
  -c     catalog to report [COG, KOG, XOG], XOG for [KC]OG; default: [COG]
  -v     print version number

示例:

emapper-eggnog -c COG   nr.emapper.annotations

到目前为止, 我们可以根据eggnog-mapper完成KEGG/BiGG/COG/KOG/GO的功能注释。

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

Last Update: 2017-12-01 11:09 AM

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%)

       504  OTUs found in all samples (14.1%)
       739  OTUs found in 90% of samples (20.7%)
      1808  OTUs found in 50% of samples (50.7%)

Sample sizes: min 48273, lo 59502, med 60943, mean 60178.2, hi 61976, max 62871

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

Last Upate: 2017-11-29 10:18 PM

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

一、引言

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

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

一、引言

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

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

当安装这些小工具不是问题的时候,使用这些小工具就成为最大的问题了,首当其冲就是了解命令行的使用说明,最直接的肯定是 man 工具, 如下图,我们可以在命令行下输入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
tar xzvf  tldr_linux_amd64.tar.gz
mv  tldr_linux_amd64 tldr-1.0
export PATH=$PATH:$PWD/tldr-1.0

四、使用

命令行接口(以Node.js客户端为例):

$ tldr

  Usage: tldr command [options]
  Simplified and community-driven man pages
  Options:
    -h, --help               output usage information
    -V, --version            output the version number
    -l, --list               List all commands for the chosen platform in the cache
    -a, --list-all           List all commands in the cache
    -1, --single-column      List single command per line (use with options -l or -a)
    -r, --random             Show a random command
    -e, --random-example     Show a random example
    -f, --render [file]      Render a specific markdown [file]
    -m, --markdown           Output in markdown format
    -o, --os [type]          Override the operating system [linux, osx, sunos]
    --linux                  Override the operating system with Linux
    --osx                    Override the operating system with OSX
    --sunos                  Override the operating system with SunOS
    -t, --theme [theme]      Color theme (simple, base16, ocean)
    -s, --search [keywords]  Search pages using keywords
    -u, --update             Update the local cache
    -c, --clear-cache        Clear the local cache

Examples

    $ tldr tar
    $ tldr du --os=linux
    $ tldr --search "create symbolic link to file"
    $ tldr --list
    $ tldr --list-all
    $ tldr --random
    $ tldr --random-example

To control the cache

    $ tldr --update
    $ tldr --clear-cache

To render a local file (for testing)

    $ tldr --render /path/to/file.md

下面介绍几种功能性可选项:

-l: 列出缓存中指定操作系统的所有命令;
-a: 列出缓存中所有命令;
-m: 输出格式为markdown 格式;
--linux/--osx/--sunos:重载操作系统类型
-t: 颜色主题模式;
-u: 升级本地缓存;
-c:清除本地缓存;

示例:

$ tldr tar

  tar

  Archiving utility.
  Often combined with a compression method, such as gzip or bzip.

  - Create an archive from files:
    tar cf target.tar file1 file2 file3

  - Create a gzipped archive:
    tar czf target.tar.gz file1 file2 file3

  - Extract an archive in a target folder:
    tar xf source.tar -C folder

  - Extract a gzipped archive in the current directory:
    tar xzf source.tar.gz

  - Extract a bzipped archive in the current directory:
    tar xjf source.tar.bz2

  - Create a compressed archive, using archive suffix to determine the compression program:
    tar caf target.tar.xz file1 file2 file3

  - List the contents of a tar file:
    tar tvf source.tar

很方便的一个小工具,如果忘记一个命令具体使用,请使用:tldr command

所支持的不同平台的命令行使用用例,请下载:

tldr-book

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

Last Update: 2017/11/27 21:34

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.dmpname.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

可以通过 -c 指定列进行翻译。

测试文件(test.tsv):

read_1  68      Lysobacter enzymogenes

第一列为序列名称, 第二列为NCBI Taxon ID, 第三列为科学命名。

对第二列进行翻译:

taxon-translate -c 2  ncbi.map.gz test.tsv

我们可以获得:

read_1  69      Lysobacter enzymogenes  \
d:Bacteria;p:Proteobacteria;c:Gammaproteobacteria;  \
o:Xanthomonadales;f:Xanthomonadaceae;g:Lysobacter;s:Lysobacter enzymogenes

很直接,很简单。

该程序每次都要将ncbi.map.gz 读入内存, 因此可以先将ncbi.map.gz对于对应的数据结构序列化,然后直接内存映射,这样可以对taxon-translate进一步优化。

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

Last Update: 2017-11-22 5:56 PM

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

今天开始会介绍一系列Metagenome(以后使用元基因组)数据分析中的序列分类问题,这类问题一直是研究的热点和难点。

热点:有用,比如直接鉴定环境样本(粪便、痰液、肺积液)病原,算法多样性, 从基于碱基组成统计分类、基于核酸Kmer查找、基于氨基酸序列MEM查找到序列比对,各种算法八仙过海、各显神通。

难点:序列分类准确性、敏感性、内存有效性、计算密集性等都是待解决问题。

Kaiju是一款基于在氨基酸空间进行相似性搜索的序列分类算法。

1. 算法介绍

Kaiju算法发表在了Nature Communications杂志 :Fast and sensitive taxonomic classification for metagenomics with Kaiju

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回溯。

几点说明:

  1. BLOSUM62 矩阵, ORF Score 以及后续的MEM左侧延伸后Score值计算都是按照 BLOSUM62计算。
  2. Kaiju的Greedy模式敏感度要比MEM模式高,但是牺牲了分类速度,对于一些新的基因尤为明显,Greedy采用的模式是启发式算法,只对MEM的末端延伸。
  3. LCA算法见下图,从所有叶子节点向上回溯,找到所有命中节点的共同祖先,比如a|b|c叶子节点的LCA节点就是x节点。

LCA

图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.

这里统计序列分类有几点需要注意:

  1. 输出结果为指定分类水平的相对丰度,这里需要注意是分母的问题, 可以包含未分类序列也可以不包含未分类序列, 通过 -u 切换。
  2. 可以将特定分类水平 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

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

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

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

下载地址: 64位 / 32位

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

image.png-91.3kB

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

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

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

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

MobaXterm界面

2.1 工具下载

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

2.2 推荐理由

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

2.3 SSH远程登录

Mobaxterm

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

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

第一次登录会提示输入用户密码,并会提示是否保存密码,选择后下次自动登录。

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

BLAST类序列比对工具在生物学数据分析中应用最为广泛,比如16S全长序列比对(MEGABLAST), 氨基酸功能注释(BLASTP)等, 另外还有各种快速比对工具,比如: DIAMONDGHOSTXRAPSearch2等,这些工具默认都支持制表符分隔的输出文件(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_utilsgithub 仓库) 小工具集合, 主要使用场景包括:

  1. 根据E值、Bit score值、相似度筛选所有满足条件的结果, blast_hits
  2. 根据E值、Bit score值、相似度筛选满足条件的最佳HSP, blast_hsp
  3. 选取每条查询序列的Top Hsp, best_hsp
  4. 对目标序列进行注释, blast_annotation

下面对每个程序使用进行介绍:

1. blast_hits

命令行接口

$ blast_hits
Usage: blast_hits  [options]  <blast|->
Options:
  -b DOUBLE  MIN bit score, default: [60]
  -e DOUBLE  MAX E-value, default: [0.001]
  -i DOUBLE  MIN identity, default: [0]
  -v print version number

根据上述描述, blast_hits 接受三个可选参数, bit score 、E-value 和 identity,就是输出所有满足条件的比对结果。

实例:

blast_hits -b 60  -e 0.001 -i 60  S7711.tsv

2. blast_hsp

命令行接口

$ blast_hsp
Usage: blast_hits  [options]  <blast|->
Options:
  -b DOUBLE  MIN bit score, default: [60]
  -e DOUBLE  MAX E-value, default: [0.001]
  -i DOUBLE  MIN identity, default: [0]
  -v print version number

根据上述描述, blast_hsp 接受三个可选参数, bit score 、E-value 和 identity,就是输出满足条件的第一个hsp。

实例:

blast_hsp -b 60  -e 0.001 -i 60  S7711.tsv

3. best_hsp

命令行接口

$best_hsp
Usage: best_hsp <blast|->
version: 0.0.2

$best_hsp 接受一个参数, 最输入文件不进行限制,只输出第一个hsp。

实例:

best_hsp S7711.tsv

4. blast_annotation

因为我们有时需要对目标序列进行各种关联,需要使用目标序列的标识符信息,所以:

LDEE01000001.1.2    sp|P33330|SERC_YEAST

可能会比:

LDEE01000001.1.2        sp|P33330|SERC_YEAST Phosphoserine aminotransferase ...

更方便;

blast_annotation 做的事情就是对sp|P33330|SERC_YEAST 进行注释,变成 sp|P33330|SERC_YEAST Phosphoserine aminotransferase ... 模式。

程序也比较简单, 需要提供序列->注释的映射文件。

实例:

blast_annotation  annotation.txt  S7711.tsv  >S7711_ann.tsv

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

Last update:2017-11-15 7:07 PM

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

一、引子

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

SNP Pipeline

该调查数据来自 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 和 seqtk_sra

2.2 去除低质量序列

seqtk trimfq - 

steaming 模式切除低质量序列。

内部链接: NGS数据质量过滤:Usearch fastq-filter、seqtk trimfq和seqtk_filter

2.3. 序列比对

至此,snippy 预处理工作已经做完,现在可以进入snippy的工作模式,这里只记录主要步骤。

gbk转gff和fasta:

使用`BioPerl`从GBK文件中提取fasta和gff文件用于构建序列比对索引和snpEff库。

构建基因组文件序列

samtools faidx reference/ref.fa

构建基因组索引:

bwa index reference/ref.fa

序列比对:

bwa mem  -v 2 -M -R '@RG    ID:ERR114970    SM:ERR114970' \
    -p -t 16 reference/ref.fa /dev/stdin | \
    samtools view -@ 16 -F 3844 -q 60 -S -b -u -T reference/ref.fa - |\
    samtools sort -O bam -o ERR114970.bam -@ 16 -

构建bam文件索引:

samtools index ERR114970.bam

计算比对深度:

samtools depth -q 20 ERR114970.bam | bgzip > ERR114970.depth.gz

这一步完成了序列比对工作,使用的都是比较常规的操作。

2.4. Variation calling

freebayes-parallel reference/ref.txt 16 -p 1 -q 20 -m 60 \
     --min-coverage 10 -V -f \
     reference/ref.fa ERR114970.bam > ERR114970.raw.vcf

使用freebayes (Bayesian haplotype-based polymorphism discovery)鉴定 Variation, Freebayes的单倍体模型比较适合细菌基因组SNP/INDEL鉴定。

-q –min-base-quality, 最低碱基质量; -q 20; -m –min-mapping-quality,最低比对质量; -m 60; -p –ploidy 设置单倍体模型, -p 1; –min-coverage: 最小覆盖度。 –min-coverage 10;

snippy使用的是freebayes并行模式 freebayes-parallel

需要使用分割regions文件, 使用fasta_generate_regions.py 生成。

fasta_generate_regions.py reference/ref.fa.fai 80866 > reference/ref.txt

Variation 过滤:

snippy-vcf_filter --minqual 10 --mincov 10 --minfrac 0.9  ERR114970.raw.vcf > ERR114970.filt.vcf

使用下面参数控制variation 质量:

–mincov :位点最低覆盖度, 10 –minfrac :alt碱基最低比率,0.9

2.5. SNP/INDE 注释

创建配置文件和snpEff build库文件,然后使用snpEff

snpEff build -c reference/snpeff.config -dataDir . -gff3 ref
snpEff ann -no-downstream -no-upstream -no-intergenic -no-utr -c reference/snpeff.config -dataDir . -noStats ref ERR114970.filt.vcf > ERR114970.vcf

使用-no-downstream -no-upstream -no-intergenic -no-utr 过滤掉基因上游、下游、间区以及UTR区域的变异位点, 生成注释的VCF文件。

至此我们已经可以获得从SRA文件到VCF的细节过程。

2.6. reference guided assembly

snippy 提供了另外一个我们称之为”reference guided”的拼装,实现方式为:

vcf-consensus ERR114970.vcf.gz < reference/ref.fa > ERR114970.consensus.fa

三、总结

3.1. 设计程序时一定要注意区分:

标准输出流: stdout 标准错误流: stderr

准则: 如果不是结果文件,一律使用标准错误流打印。

3.2 Pipeline 设计

很多工作流,其实细节都很好实现,工作流主要起的作用就是封装,让接口更加简单, snippy 集成了 bioperl、bgzip、bwa、freebayes、samtools、tabix、vcfallelicprimitives、vcfstreamsort、vcfuniq等外部程序,使用这些程序可以使用相对路径调用,也可以设置PATH路径,但是最佳实践是环境变量只有程序能够访问,这样就不会造成程序版本污染。

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

Last update:2017-11-14 7:46 PM

使用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. 烧制

win32diskimager

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

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

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

5. 系统安装

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