Archives

Categories

合并和拆分双端序列:seqtk_interleave和seqtk_deinterleave

seqtk_interleaveseqtk_deinterleaveseqtk_utils 系列工具中的一组对应工具,即合并PE双端序列与拆分Interleaved格式的fasta/q文件。

1. seqtk_interleave

很多应用程序需要 Interleaved 格式序列作为输入文件,比如 SortMeRNA 很多工具都可以解决这个问题,比如: seqtk mergepe 可以很方便的将两个fq文件的R1和R2转换成下列文件格式:

@HISEQ:573:HLWMHBCXY:1:1101:11739:2308 1:N:0:TGGACTAT
CCTATGGGGTGCAGCAGTGGGGAATATTGGACAATGGGGGGAACCCTGAT
+
IIIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIII
@HISEQ:573:HLWMHBCXY:1:1101:11739:2308 2:N:0:TGGACTAT
GCTGTGAGGACTACTTGGGTATCTAATCCTGTTTGCTACCCACGCTTTCG
+
BDDDDHHIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHIIIIIIIIIIIII

但是有时我们需要下面格式。

@HISEQ:573:HLWMHBCXY:1:1101:11739:2308/1
CCTATGGGGTGCAGCAGTGGGGAATATTGGACAATGGGGGGAACCCTGAT
+
IIIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIII
@HISEQ:573:HLWMHBCXY:1:1101:11739:2308/2
GCTGTGAGGACTACTTGGGTATCTAATCCTGTTTGCTACCCACGCTTTCG
+
BDDDDHHIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHIIIIIIIIIIIII

seqtk_interleave 可以输出这样old fashion格式。

命令行接口:

$ seqtk_interleave
Usage: seqtk_interleave [options] <in1.fq> <in2.fq>
Options:
  -o  old fashion interleave fastq format
  -v  print version number

使用 -o 切换到 /1 和 /2 模式

示例:

seqtk_interleave  OS2_1.R1.fq  OS2_1.R2.fq
seqtk_interleave  -o  OS2_1.R1.fq  OS2_1.R2.fq

2. seqtk_deinterleave

seqtk_deinterleave 是 seqtk_interleave 反操作, 即使interleave的结果,拆分成R1和R2双端模式。

命令行接口:

$seqtk_deinterleave
Usage: seqtk_deinterleave [options] <fasta/q> <prefix>
Options:
  -f  specify input file as fasta format
  -v  print version number

使用 -f 设定输入文件为fasta文件格式, 后缀名会默认使用 .fasta

示例:

seqtk_interleave  OS2_1.R1.fq  OS2_1.R2.fq > OS2.fastq
seqtk_deinterleave OS2.fastq  OS2
seqtk seq  -1  OS2.fastq  >OS2-1.fastq
seqtk seq  -2  OS2.fastq  >OS2-2.fastq
seqtk seq  -A  OS2.fastq  >OS2.fasta
seqtk_deinterleave -f OS2.fastq  OS2

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

Last Upate: 2017/10/21 22:47

OTU表子集操作:otutable_subset 和 otutable_subsamples

otutable_subsetotutable_subsamples 程序设计的目的是解决获取OTU表子集的操作方式, 比如:

1. 提取给定OTU,并按照指定顺序生成新的OTU;
2. 提取特定样本,并按照特定顺序生成新的OTU;

otutable_subset

Usearchotutab_otu_subset 可以完成类似的事情, 但是不能很好解决一些问题, 比如不支持带有注释的OTU表,支持streaming 模式。

命令行接口:

$ otutable_subset
Usage: otutab_subset [options] <otu-table>
Options:
  -r  reorder the OTUs by OTUs in list
  -v  print version number

默认按照 OTU表中出现的顺序输出, 如指定 -r 可选项可以按照输入文件顺序。

实例:

文件 otu_ids.txt 为我们需要提取的OTU列表:

OTU_2
OTU_1
OTU_3

命令

otutable_subset -r  otutable_counts_ann.txt   otu_ids.txt
otutable_subset -r  otutable_counts.txt   otu_ids.txt

使用场景:比如我们对注释的结果去除线粒体或者叶绿体来源的序列,需要重新生成OTU表。

otutab_otu_subset 命令行模式:

usearch -otutab_otu_subset otutable_counts.txt -labels otu_ids.txt -output subset.txt

测试场景中对OTU表要求比较高,比如OTU的编号顺序以及是否带注释等。

otutable_subsamples

Usearchotutab_otu_subset 可以完成类似的事情, 但是不能很好解决一些问题, 比如排序问题。

命令行接口:

$ otutable_subsamples
Usage: otutable_subsamples <otu-table>  <samples| F1,F2,F3,F4>
version: 0.0.1

Usearch 同样有类似命令 otutab_sample_subset

实例:
文件 sample_ids.txt 为我们需要提取的OTU列表;

Con1d-2
Con1d-1
Con1d-3

我们只关注着三个样本的信息, 我们可以这样操作:

otutable_subsamples  otutable_counts_ann.txt  sample_ids.txt
otutable_subsamples  otutable_counts_ann.txt  sample_ids.txt | otutable_filter - 0
cat  cat  sample_ids.txt  |  otutable_subsamples otutable_counts_ann.txt  - | otutable_filter - 0

otutab_sample_subset 命令行模式:

usearch -otutab_sample_subset  otutable_counts.txt   -labels sample_ids.txt -output subset.txt

测试场景中不支持带有注释的OTU表。

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

Last Upate: 2017/10/21 21:07

OTU表过滤操作:otutable_filter

otutable_filter程序主要实现了对OTU表中的丰度信息进行过滤,比如过滤符合一定规则的OTU(比如丰度和小于设定阈值)。

命令行接口:

$ otutable_filter
Usage: otutable_abundance [options] <otu-table> <cutoff:2>
Options:
  -v  print version number

程序接受两个参数, OTU表和设置阈值,支持未注释的OTU表和注释的OTU表。

实例

otutable_filter  otutable_freq.txt  0.1

过滤掉特定OTU(所有的样本相对丰度和小于0.1)的OTU。

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

Last Upate: 2017/10/21 18:51

生成不同水平系统分类丰度表:otutable_abundance 和 otutable_level

这篇文章介绍 otutable_abundanceotutable_level, 对注释的OTU表统计不同系统分类水平丰度信息表,
支持OTU频数表和相对丰度表,注释信息按照单字母缩写模式:

比如:d:Bacteria,p:Cyanobacteria/Chloroplast,c:Chloroplast,f:Chloroplast,g:Streptophyta

1. otutable_level

命令行接口

$ otutable_level
Usage: otutable_level [options] <otu-table:annotated>
Options:
  -l  taxonomy level in [dkcofgs] Default:['g']
  -v  print version number

通过 -l 指定分类水平, 默认是 genus 水平,映射表如下:

k   Kingdom o   Order
d   Domain  f   Family
p   Phylum  g   Genus
c   Class   s   Species

程序接口简单接受一个参数: OTU注释表。

实例

otutable_level  -l g  otutable_counts_ann.txt

2. otutable_abundance

otutable_abundance 多增加一个参数,用于计算单个样本指定分类水平的丰度信息。

命令行接口

$ otutable_abundance
Usage: otutable_abundance [options] <otu-table:annotated> <sample>
Options:
  -l  taxonomy level in [dkcofgs] Default:['g']
  -v  print version number

实例

$ otutable_abundance  -l g  otutable_counts_ann.txt  WG7d-1 >/dev/null

生成丰度表可以制作各种柱状图、热图以及各种统计分析。

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

Last Upate: 2017/10/21 18:35

OTU表生成 Krona 输入文件:otutable_krona

这篇文章介绍如何根据注释的 OTU Table 进行 Krona 可视化

Krona

1. Krona介绍

Krona是一种基于HTML5的交互式数据可视化方式,兼容大部分现代浏览器,比如 Chrome,Firefox,IE等。

参考文件:Interactive metagenomic visualization in a Web browser Github仓库地址: Krona: Interactively explore metagenomes and more from a web browser.

KronaTools 可本地安装,并包含各种程序进行导入不同工具的输出结果,比如:RDP、BLAST、MEGAN等, Krona的各种使用我们通过另外一篇文章进行描述,现在主要使用Krona的基本功能。 ktImportText

主要命令行接口:

$ ktImportText

ktImportText \
   [options] \
   text_1[,name_1] \
   [text_2[,name_2]] \
   ...

   [-o <string>]  Output file name. [Default: 'text.krona.html']

   [-n <string>]  Name of the highest level. [Default: 'all']

   [-q]           Files do not have a field for quantity.

   [-c]           Combine data from each file, rather than creating separate datasets within the chart.

   [-u <string>]  URL of Krona resources to use instead of bundling them with the chart (e.g. "http://krona.sourceforge.net"). \
                  Reduces size of charts and allows updates, though charts will not work
                  without access to this URL.

我们只需要生成下面格式文件,就可以将数据进行Krona可视化:

3   d:Bacteria  p:Firmicutes    c:Bacilli   o:Bacillales    f:Paenibacillaceae_1    g:Cohnella

制表符分隔, 第一列为数值,后面为等级分类;

2. otutable_krona

我们实现了otutable_krona, ( 程序链接 )可直接从OTU表,生成指定文件的 Krona输入文件

otutable_krona 命令行接口:

$ otutable_krona
Usage: otutable_krona <otu-table:annotated> <sample>
version: 0.0.2

需要指定两个参数, 第一个为带有注释的OTU表,第二个参数为样本;

命令行示例: otutable_krona otutable_freq_ann.txt Con1d-1 >Con1d-1.txt ktImportText Con1d-1.txt:Con1d-1 -o Con1d-1.html

otutable_krona 对生成的krona 文件做了特殊处理, 尝试对于一个分类如果包含子分类,当前分类会添加_unclassified 标识。 这些分类一般代表一些比较新的、未分类的菌,对于研究寻找新菌来说特别有帮助。

unclassified

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

Last Upate: 2017-10-20 4:49 PM

对OTU表进行注释:otutable_annotation

这篇文章介绍小工具 otutable_annotation, 这个工具设计比较通用,可以对 usearch sintax 或者 usearch utax 、RDP、 DADA2的分类结果进行OTU, 不过需要需要一些适配器才能完成工作,下面介绍如何结合其它工具完成各种OTU表注释。

1. 输入文件格式:

OTU表格式为标准QIIMME OTU 表格式, 包含以 # 开头的注释行;

#OTU ID F_5 F_6 F_7 F_8
OTU_1   1763    8408    1420    2674
OTU_2   1582    7132    1234    2317
OTU_3   0   10  1   15

注释文件需要这样格式:

OTU_1   d:Bacteria,p:Cyanobacteria/Chloroplast,c:Chloroplast,f:Chloroplast,g:Streptophyta

我们可以通过特定适配器灵活完成注释。

otutable_counts.txt  OTU表;
otu.sintax.txt  sintax分类结果;
otu.rdp.txt 分类结果;
dada2.taxa.txt 分类结果;

2. 命令行接口:

$otutable_annotation
Usage: otutable_annotation <sintax|utax> <otu-table>
version: 0.0.2

3. 各种分类器实例:

3.1 sintax|utax 分类注释

sintax 结果为四列文件, 第四列为序列分类文件;

格式:

第一列:OTU编号     OTU_1
第二列:分类信息(括号包含可信度值)
        d:Bacteria(1.0000),p:Cyanobacteria/Chloroplast(1.0000),c:Chloroplast(1.0000), \
        f:Chloroplast(1.0000),g:Streptophyta(1.0000)       
第三列:strand信息  +
第四列:分类结果,使用了sintax_cutoff进行了过滤;
        d:Bacteria,p:Cyanobacteria/Chloroplast,c:Chloroplast,f:Chloroplast,g:Streptophyta

我们只需要第一列和第四列, 使用 cut 获取目标列;

执行示例:

cut -f1,4 otu.sintax.txt | otutable_annotation - otutable_counts.txt >otutable_counts_ann.txt

3.2 rdp 分类注释

rdp 分类结果格式比较复杂,我们是需要使用 otutable-rdp 完成数据准备; rdp 输出结果:

OTU_1       Root    rootrank    1.0 Bacteria    domain  1.0 Cyanobacteria/Chloroplast   \
            phylum  1.0 Chloroplast class   1.0 Chloroplast family  1.0                 \
            Streptophyta    genus   1.0

otutable-rdp 命令接口:

$ otutable_rdp
Usage: otutab_rdp [options] <text>
Options:
  -c  confidence cutoff Default:[0.8]
  -v  print version number

可以通过:

$otutable_rdp -c 0.8 otu.rdp.txt | head -n 1
OTU_1   d:Bacteria,p:Cyanobacteria/Chloroplast,c:Chloroplast,f:Chloroplast,g:Streptophyta

完成格式化。

结合 otutable_annotation 可以对 OTU 表进行注释

otutable_rdp -c 0.8 otu.rdp.txt | otutable_annotation – otutable_counts.txt >otutable_counts_ann.txt

3.3 dada2 分类注释

DADA2 工作环境是R平台, 在R平台, 可以使用 cbind 合成 丰度表和注释信息,下面介绍下 otutable-dada2 如何转换成,我们需要的格式:

$ otutable_dada2
Usage: otutab_dada2 [options] <text>
Options:
  -i  ignore the title line;
  -v  print version number

DADA2 结果使用下面命令产生:

taxa.fmt <- cbind('#seq'=rownames(taxa), taxa)
write.table(taxa.fmt , "dada2.taxa.txt", sep="\t", quote=F, row.names = F)

为了普适性, 我们包含了表头信息;
格式化处理:

otutable_dada2 -i  dada2.taxa.txt | head

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

使用SINTAX对16S序列进行分类:SINTAX

一、SINTAX介绍

先前介绍过RDP分类器的使用:使用RDP对16S序列进行分类:rdp_classifier, 现在介绍另外一个基于非朴素贝叶斯分类算法: SINTAX(SImple Non-Bayesian TAXonomy) , 文章见:

SINTAX: a simple non-Bayesian taxonomy classifier for 16S and ITS sequences.

SINTAX使用 Kmer (默认k=8, 放回式采样32 kmers)去计算和参考库共享Kmer, 确定最佳Hits, 迭代100次, 确定每个Level的分类(出现次数最多的分类)和可信度(出现的频率)。由此可见并不需要训练参考序列库。

后面会继续介绍 DADA2 assignTaxonomyassignSpecies
QIIME2 的 q2-feature-classifier,新文章: Optimizing taxonomic classification of marker gene sequences

二、SINTAX 实践:

1 数据库下载

SINTAX提供了 RDP training set 16 (13k seqs, with species names )SILVA 123 (1.6M seqs.)以及LTP库和Greengenes 13.5 (1.2M seqs.)

SINTAX库不需要额外训练,但是需要按照特定格式组织(具体要求),比如:

>X71857_S000021696;tax=d:Bacteria,p:Firmicutes,c:Clostridia,o:Clostridiales,f:Clostridiaceae_1,g:Clostridium_sensu_stricto,s:Clostridium_puniceum;

按照下面表进行缩写可以给下游分析提供很大便利,后续会涉及的物种分类都会按照这种格式组织,包括RDP的分类结果。

k   Kingdom o   Order
d   Domain  f   Family
p   Phylum  g   Genus
c   Class   s   Species

2 格式化数据库

使用makeudb_usearch 对下载的数据库进行构建索引UDB文件。

usearch -makeudb_sintax   rdp_16s_v16_sp.fa  -output  rdp_16s_v16_sp.udb

rdp_16s_v16_sp.fa 为推荐使用库,大部分的物种分类虽然不能准确预测,但是对于一些物种还是可以进行预测,特别是人类微生物组数据,DADA2 也提供了 assignSpecies 函数预测种水平信息。

3 序列分类

SINTAX分类器命令行接口:

usearch -sintax otus.fa -db rdp_16s_v16_sp.udb      \
        -threads 24 -strand both -sintax_cutoff 0.8 \
        -tabbedout otu.sintax 

几个参数:

-threads 指定线程数,默认使用10线程
-strand  plus或者both
-sintax_cutoff 使用0.8
-tabbedout table格式输出

输出文件为四列:

第一列:OTU编号     OTU_1
第二列:分类信息(括号包含可信度值)
        d:Bacteria(1.0000),p:Cyanobacteria/Chloroplast(1.0000),c:Chloroplast(1.0000), \
        f:Chloroplast(1.0000),g:Streptophyta(1.0000)       
第三列:strand信息  +
第四列:分类结果,使用了sintax_cutoff进行了过滤;
        d:Bacteria,p:Cyanobacteria/Chloroplast,c:Chloroplast,f:Chloroplast,g:Streptophyta

从本例可以看到分类结果为来自叶绿体序列,正常分析流程,应该包含去除线粒体和叶绿体的参数。

4 物种分类统计

使用sintax_summary可以统计不同分类水平的丰度信息:

usearch -sintax_summary otu.sintax  -otutabin otutable_counts.txt -output phylum_summary.txt -rank p

后续可以进行各种数据可视化分析

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

Last Update: 2017-10-20 11:19 AM

多tsv文件组合操作: tabtk_join 和 tabtk_agg

今天介绍两个合并 tsv文件 的工具: tabtk_jointabtk_agg

1. tabtk_join   组合多个文件, key设置为第一列
2. tabtk_agg    指定Key列(可多列)和目标列,将多个文件的目标列合并到同一个文件

下面以下图为示例介绍如何使用tabtk_join和tabtk_agg完成类似操作

tabtk_utiks

一、tabtk_join 表合并操作

我们经典使用场景:

组合不同16S多样性实验的结果(比如不同文章的不同16S分区的物种分类信息),不推荐使用Closed-OTU方式构建新的OTU表, 最合适的方式是单独构建OTU表,然后合并不同水平的分类结果。

命令行接口:

$tabtk_join
Usage: tabtk_join [options] 
Options:
  -p CHAR  placehold for missing value, default: ['-'];
  -v       print version number

可选项:

-p 占位符,数据缺失时使用占位符代替, 数值型可以使用'0', 文本型使用'-'

示例:

$cat A-phylum.txt

#level  F_1 F_2
Fusobacteria    0.001872    0.000248
Firmicutes  0.85635 0.827091
Proteobacteria  0.082811    0.092849
Acidobacteria   0.001045    0.002109
Bacteroidetes   0.004354    0.003442
Above_phylum    0.001915    0.000713
Chloroflexi 0.000044    0.000093
Actinobacteria  0.051158    0.071581
Candidatus_Saccharibacteria 0.000435    0.000372
Gemmatimonadetes    0   0.000031

$tabtk_join -p 0  A-phylum.txt  B-phylum.txt

约定:

文件必须使用'#'作为表注释, 这个是很好的习惯,而且后续所有程序都会默认忽略掉以'#'开头的注释行

二、tabtk_agg 表汇总操作

我们经典使用场景:

tabtk_agg 使用场景更多,比如多做样本数据分析是,单样本分析完,需要对分组样本进行如表数值汇总操作,比如:执行完非冗余基因集合丰度定量,需要合并成多样本表丰度表,这样可以使用tabtk_agg轻松完成。

命令行接口:

$ tabtk_agg
Usage: tabtk_agg [options] [label:text ...]
Options:
  -k STR  the keys fields pattern: 1:2:3, default: [1];
  -t STR  the titles for keys: key_1:key_2:key_3, default: [catalog];
  -c INT  the target column default: [2];
  -p CHAR placehold for missing value: default ['0'];
  -i      ignore the head line;
  -v      print version number

tabtk_agg 命令行接口稍微复杂了点,可选项多点:

-k 指定(组合) 列作为合并的主 key , 默认第一列;
-t 指定(组合) 列的表头,单列默认catalog,如果指定多列作为key,默认使用key_1, key_2 \          
   描述表头,可以显示指定. -k的个数要和-t的个数对应;
-c 汇总目标列,默认使用第二列;
-p 占位符, 填充缺失数据
-i 忽略标题行

大部分情况需样本信息,而结果中一般没有样本名称信息,所以在指定文件的时候需要额外添加一个样本信息。

示例(汇总多个样本kallisto丰度定量结果):

$cat  A-abundance.tsv | head -n 4
target_id       length  eff_length      est_coun
nr_1    309     16.3657 0       0
nr_2    318     17.8465 0       0
nr_3    351     40.6759 0       0

$time  tabtk_agg -i -p 0 -k 1 -t tpm -c 4  A:A-abundance.tsv  B:B-abundance.tsv |  head -n 2
#tpm    A       B
nr_1    0       0

real    0m1.539s
user    0m1.415s
sys     0m0.126s

2秒可以完成100万行列汇总操作;

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

Last Update: 2017-09-25 04:34 PM

计算OTU表相对丰度:otutable_freq

本文介绍otutable-utils系列的一个小工具:otutable_freq (Github 链接 ), 就是计算OTU表的相对丰度(频数转频率).

1. 命令行接口

$ otutable_freq
Usage:  otutable_freq <matrix>
version 0.0.2

2. 命令行接口

otutable_freq otutable_counts.txt >otutable_freqs.txt

3. 其它实现

3.1 otutab_counts2freqs

otutab_counts2freqs 是 Usearch 9.2 版本才加入的命令行工具,接口也比较简单

usearch -otutab_counts2freqs otutable_counts.txt -output otutable_freqs.txt

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

Last Upate: 2017-10-18 6:35 PM

通过 minhash 指纹快速进行菌种鉴定:MASH 篇

前面一个文章提到了如何下载数据,如何并行下载NCBI拼装数据库的Genbank和Refseq数据, 现在开始要讲讲如何使用这些数据。

一、 mash介绍

发表在 Genome Biology 杂志 Mash: fast genome and metagenome distance estimation using MinHash。 可以通过Github Mash 仓库 下载可执行文件;

wget https://github.com/marbl/Mash/releases/download/v2.0/mash-Linux64-v2.0.tar
tar xvf mash-Linux64-v2.0.tar
mv  mash-Linux64-v2.0  mash-2.0

MinHash介绍可以参考:Wikipedia MinHash

Mash原理图

原理图

通过构建序列(基因组,原始reads)Kmer集合,并计算出kmer 的Hash集合, 然后通过杰卡德指数Jaccard Index 估计相似度, 即 A和B共享的Hash集合(Hash值可以通过Hash函数计算)与A和B所有的不同Hash集合的比值。

另外: 传统的DDH杂交算法也可以很好的进行菌株鉴定,但是速度没有使用MinHash快,下次介绍几种DDH相关工具。

二、命令行接口

主程序为单一可执行文件程序,分发比较方便,注册主程序命令接口:

$ mash

Mash version 2.0
Type 'mash --license' for license and copyright information.

Usage:
  mash <command> [options] [arguments ...]
Commands:
  bounds  Print a table of Mash error bounds.
  dist    Estimate the distance of query sequences to references.
  info    Display information about sketch files.
  paste   Create a single sketch file from multiple sketch files.
  screen  Determine whether query sequences are within a larger pool of
          sequences.
  sketch  Create sketches (reduced representations for fast operations).

现在主要介绍的功能模块有: sketch, paste,dist

sketch 创建sketches指纹;
paste  合并已创建指纹;
dist   计算指纹之间的相似度;

dist子命令产生的结果可以用来构建系统进化树,该内容这篇文章不会涉及; sketch/paste/dist/screen 都支持文件列表输入;

三、用例介绍

这里使用我们测的Mycobacterium abscessus几个菌株为示例:

使用其中6株菌,2/3/7/8/9 作为库,使用129进行查询;

$ ls examples
129.fasta  2.fasta  3.fasta  7.fasta  8.fasta  9.fasta

3.1 构建sketch库

$ cd examples
$ mash sketch -p 20 -o  examples.msh   2.fasta  3.fasta  7.fasta  8.fasta  9.fasta
$ mash info  examples.msh

Header:
  Hash function (seed):          MurmurHash3_x64_128 (42)
  K-mer size:                    21 (64-bit hashes)
  Alphabet:                      ACGT (canonical)
  Target min-hashes per sketch:  1000
  Sketches:                      5
Sketches:
  [Hashes]  [Length]  [ID]     [Comment]
  1000      5030882   2.fasta  -
  1000      4918032   3.fasta  -
  1000      5519232   7.fasta  -
  1000      4777861   8.fasta  -
  1000      5134600   9.fasta  -

使用mash info 可以查看sketch信息。

3.2 计算菌株相似度

mash dist  -v 0.05 -p 20  examples.msh   129.fasta | sort -k3n |   mash_dist -  | head
9.fasta 129.fasta       0.00904645      0       705/1000
2.fasta 129.fasta       0.0103562       0       673/1000
7.fasta 129.fasta       0.011621        0       644/1000
3.fasta 129.fasta       0.0136496       0       601/1000
8.fasta 129.fasta       0.0255655       0       413/1000

mash dist 结果, 第三列为距离值,越小越好, 第四列为P-value值,越小可信度越高, 第五列为共享hashes数,使用距离信息可以判断菌种相似度。

四、构建NCBI Genbank/Refseq sketch库

基因组比较多,需要使用脚本进行并行构建sketch库

4.1 构建sketch库

mash-sketch gca.txt  /project/org-inspect/db/GCA  GCA
mash-sketch gcf.txt  /project/org-inspect/db/GCF  GCF

4.2 合并sketch文件

find  GCA  -name '*.msh' | perl -ane '$t=`readlink -f $F[0]`; print $t'  >gca-mash.list
mash-paste-batch  gca-mash.list  gca.msh

find  GCF  -name '*.msh' | perl -ane '$t=`readlink -f $F[0]`; print $t'  >gcf-mash.list
mash-paste-batch  gcf-mash.list  gcf.msh

sketch 文件比较多,需要采用分而治之策略,每1000个合并成个1个sketch文件, 然后对合并后的文件再做最后一次合并生成最终的mash sketch文件用于菌株鉴定。

五、菌株鉴定

现在使用我们已经测的Mycobacterium abscessus菌株为例,介绍如何通过MASH进行菌株鉴定,可以使用拼装的基因组也可以使用原始reads;

5.1 拼装基因组模式:

$mash dist  -v 0.05 -p 20   ../gca.msh  129.fasta | sort -k3n |   mash_dist -  | head -n 5
GCA_900136885.1_10625_5_56_genomic.fna.gz   129.fasta   0.00103256  0   958/1000
GCA_001217125.1_7396_7_42_genomic.fna.gz    129.fasta   0.00108342  0   956/1000
GCA_900139505.1_10625_5_44_genomic.fna.gz   129.fasta   0.00144387  0   942/1000
GCA_900139565.1_10625_5_51_genomic.fna.gz   129.fasta   0.00144387  0   942/1000
GCA_900140145.1_10665_2_61_genomic.fna.gz   129.fasta   0.00181236  0   928/1000

从结果可以看到排在第一位的为: GCF_900136885.1, Mycobacterium abscessus subsp. abscessus (high GC Gram+)

对12万个sketches指纹库进行搜索速度也是很快。

$time  mash dist  -v 0.05 -p 20   ../gca.msh  129.fasta | sort -k3n |   mash_dist -  >/dev/null

real    0m5.044s
user    0m8.101s
sys     0m2.547s

5.2 reads 模式:

合并双端序列:

seqtk mergepe 129_R1.fastq.gz  129_R2.fastq.gz  >129.fastq

计算相似度:

可以分两步执行,第一步构建sketch, 第二步计算相似度。

$time  mash sketch -r   -o 129.msh  129.fastq
Sketching 129.fastq...
Estimated genome size: 9.0162e+07
Estimated coverage:    10.721
Writing to 129.msh...

real    1m37.118s
user    1m36.602s
sys     0m0.521s

因为原始reads肯定含有一些噪音kmer, 拼装过程基本可以过滤这些kmer(一般低丰度,会造成估计基因组过大),推荐先拼装再鉴定菌株;

$ time mash dist -r -v 0.05 -p 20   ../gca.msh  129.msh  | sort -k3n |   mash_dist -  | head -n 5
GCA_001217125.1_7396_7_42_genomic.fna.gz        129.fastq       0.111073        5.90289e-219    51/1000
GCA_900139505.1_10625_5_44_genomic.fna.gz       129.fastq       0.111073        1.82741e-218    51/1000
GCA_900139565.1_10625_5_51_genomic.fna.gz       129.fastq       0.111073        1.83521e-218    51/1000
GCA_900136885.1_10625_5_56_genomic.fna.gz       129.fastq       0.11197 2.07645e-214    50/1000
GCA_001606335.1_ASM160633v1_genomic.fna.gz      129.fastq       0.112887        7.21657e-209    49/1000

real    0m2.773s
user    0m7.404s
sys     0m1.335s

从结果可以看出效果没有拼装后再鉴定效果好,我们也可以使用 screen 进行查看reads的主要分类:

$ time mash dist -r -v 0.05 -p 20   ../gca.msh  129.msh  | sort -k3n  |   mash_dist - | head -n 5
GCA_001217125.1_7396_7_42_genomic.fna.gz        129.fastq       0.111073        5.90289e-219    51/1000
GCA_900139505.1_10625_5_44_genomic.fna.gz       129.fastq       0.111073        1.82741e-218    51/1000
GCA_900139565.1_10625_5_51_genomic.fna.gz       129.fastq       0.111073        1.83521e-218    51/1000
GCA_900136885.1_10625_5_56_genomic.fna.gz       129.fastq       0.11197 2.07645e-214    50/1000
GCA_001606335.1_ASM160633v1_genomic.fna.gz      129.fastq       0.112887        7.21657e-209    49/1000

real    0m3.269s
user    0m6.165s
sys     0m1.956s

结果和reads 模式 mash dist 结果差不多。

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