Mash: 使用MinHash快速估算基因距離

工具介紹

Mash擴展了MinHash降維技術,使其成對的突變距離和P值顯著性檢驗,從而可以有效地聚類和搜索大量序列集合。混搭將大序列和序列集還原為較小的代表性草圖,進而可以快速估計基因組之前全局突變距離。MASH能用于聚類NCBI RefSeq中所有的基因組,基于組裝或者它的reads的實時的基因組聚類,聚類大量的宏基因組數據集。 該工具最終發(fā)表在Genome Biology上。

工具下載與基本使用

Github上很貼心,可以直接下載已經編譯好的版本:

wget https://github.com/marbl/Mash/releases/download/v2.2/mash-Linux64-v2.2.tar

解壓之后就能直接使用了

tar -zxvf mash-Linux64-v2.2.tar
cd mash-Linux64-v2.2/
./mash

我們簡單看看其幫助文檔:

Mash version 2.2

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 mixture of
            sequences.

  sketch    Create sketches (reduced representations for fast operations).

  triangle  Estimate a lower-triangular distance matrix.

Mash 提供兩個基本功能用于序列之間的比對 sketchdist。 sketch 功能將序列轉化為哈希結構圖。 dist 功能比對 序列之前的sketches結果并且返回 Jaccard index的估計值,P值, 還有Mash距離, 這些值能用于估計一個簡單進化模型的序列突變率。

具體使用例子

理論說了一大通,現在到大家最關鍵的環(huán)節(jié),這個工具究竟可以怎么用?

例子一:基因組之間的遺傳距離對比

下載兩個Ecoil的測試數據:

wget https://gembox.cbcb.umd.edu/mash/genome1.fna
wget https://gembox.cbcb.umd.edu/mash/genome2.fna

由于基因組比較小,這里直接調用mash dist的功能:

mash dist genome1.fna genome2.fna

查看結果,不出意料兩個基因組遺傳距離非常接近:

###Reference-ID, Query-ID, Mash-distance, P-value, and Matching-hashes
genome1.fna   genome2.fna     0.0222766       0       456/1000

當然如果我們是要比對差別很大的,兩個或者多個很大的基因組,直接使用dist功能就會很慢。這里mash的手冊上推薦先使用sketch功能轉化基因序列:

mash sketch genome1.fna
mash sketch genome2.fna
mash dist genome1.fna.msh genome2.fna.msh

除了基因組序列外,該功能還能直接應用到序列的reads中(fastq文件)。通過計算不同fastq文件之間的遺傳距離,并將遺傳距離結果整合可視化,你可以快速的得到你的不同個體測序序列之間的遺傳關系:

關于如何具體建樹我這里就不細說,但是給大家推薦一個非常實用的Github網站,這個網站已經將Mash構樹的過程封裝好為mashtree,使用起來非常方便。Github鏈接:https://github.com/lskatz/mashtree

例子二:構建本地序列的mash數據庫進行遺傳距離分析

除了上面介紹的例子之外,mash還支持使用已知的序列構建本地的數據庫,然后利用該構建的數據庫,來對其它序列進行遺傳距離的分析。

使用剛剛下載好的genome1.fna和genome2.fna進行本地mash數據庫的構建:

mash sketch -o reference genome1.fna genome2.fna

構建好后可以使用info功能對構建好的數據庫進行檢查:

mash info reference.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:                      2

Sketches:
  [Hashes]  [Length]  [ID]         [Comment]

  1000      4639675   genome1.fna  gi|49175990|ref|NC_000913.2| Escherichia coli str. K-12 substr. MG1655, complete genome

  1000      5498450   genome2.fna  gi|47118301|dbj|BA000007.2| Escherichia coli O157:H7 str. Sakai DNA, complete genome

下載新的序列,與構建好的mash數據庫進行遺傳距離比對:

wget https://gembox.cbcb.umd.edu/mash/genome3.fna
mash dist reference.msh genome3.fna

###結果顯示genome3.fna與數據庫的比對結果,可以發(fā)現genome1和genome3遺傳具體是完全一致的(同一個樣本):
genome1.fna     genome3.fna     0       0       1000/1000
genome2.fna     genome3.fna     0.0222766       0       456/1000

例子三:與已構建好的RefSeq sketch數據庫進行遺傳距離比對

除了自己構建mash sketch數據庫外,我們當然可以下載一些已經構建好的數據庫。最常用的要數RefSeq sketch數據庫。

下載好RefSeq數據庫:

wget https://gembox.cbcb.umd.edu/mash/refseq.genomes.k21s1000.msh

由于mash沒有給出測試文件,這里從我自己服務器中大豆的數據中抽取一部分進行測試:

head -n 1600 SRR7817178_1.fq > test_1.fq
head -n 1600 SRR7817178_2.fq > test_2.fq
cat test_1.fq test_2.fq >reads.fastq

這里使用-m 2參數來進一步提高結果的質量,這個參數會忽略掉低質量的單拷貝的k-mer:

mash sketch -m 2 reads.fastq

使用已下載好的RefSeq數據庫,對序列的reads進行遺傳相關系查詢:

mash dist refseq.genomes.k21s1000.msh reads.fastq.msh > distances.tab

查看結果:

GCF_000004515.4_Glycine_max_v2.0_genomic.fna.gz reads.fastq     0.295981        2.36893e-05     1/1000
GCF_000297375.1_ASM29737v1_genomic.fna.gz       reads.fastq     0.295981        2.25234e-05     1/1000
GCF_000400935.1_ASM40093v1_genomic.fna.gz       reads.fastq     0.295981        2.16808e-05     1/1000
GCF_000400955.1_ASM40095v1_genomic.fna.gz       reads.fastq     0.295981        2.16542e-05     1/1000
GCF_000837185.1_ViralProj14097_genomic.fna.gz   reads.fastq     0.295981        1.63576e-05     1/1000
GCF_001411555.1_wgs.5d_genomic.fna.gz   reads.fastq     0.295981        2.36883e-05     1/1000
GCF_001578535.1_ASM157853v1_genomic.fna.gz      reads.fastq     0.295981        1.98188e-05     1/1000
GCF_001662445.1_ASM166244v1_genomic.fna.gz      reads.fastq     0.295981        2.31897e-05     1/1000

首先可以看到這些reads是和大豆(G.max)有較近的遺傳距離,另外reads里面還含有一些污染,與其它細菌病毒有關系。

這里我們可以發(fā)現,mash是可以將reads的污染檢測出來。事實上,mash提供了screen功能專門用于污染檢測:

mash screen -w -p 4 refseq.genomes.k21s1000.msh|sort -gr |head

檢測結果如下:

0.799067        9/1000  1       6.84538e-42     GCF_001343725.1_ViralMultiSegProj188731_genomic.fna.gz  [2 seqs] NC_020234.1 Rosellinia necatrix partitivirus 2 RdRp gene for RNA-dependent RNA polymerase, co                                                             mplete cds [...]
0.798763        6/672   1       2.34994e-28     GCF_000899295.1_ViralProj176433_genomic.fna.gz  NC_018671.1 Sauropus leaf curl disease associated DNA beta, complete genome
0.783786        6/1000  3       2.57053e-27     GCF_001401365.1_ViralProj299179_genomic.fna.gz  NC_028095.1 Torulaspora delbrueckii dsRNA Mbarr-1 killer virus strain EX1180 Kbarr-1 killer toxin gene, comple                                                             te cds
0.758374        2/666   1       2.73245e-09     GCF_000922435.1_ViralProj259986_genomic.fna.gz  NC_024777.1 Small begomovirus-associated satellite isolate Sa19-S1, complete sequence
0.758338        3/1000  1       2.27756e-13     GCF_000911175.1_ViralMultiSegProj225924_genomic.fna.gz  [2 seqs] NC_022616.1 Red clover cryptic virus 1 isolate IPP_Nemaro segment 1 RNA-dependent RNA polymer                                                             ase gene, complete cds [...]
0.743837        2/1000  27      6.16326e-09     GCF_001590135.1_ViralProj314638_genomic.fna.gz  NC_029628.1 Lake Sarah-associated circular virus-21 isolate LSaCV-21-LSMU-2013, complete sequence
0.743837        2/1000  2       6.16326e-09     GCF_000928835.1_ViralProj268860_genomic.fna.gz  NC_025790.1 Black sea bass polyomavirus 1 isolate 2835, complete genome
0.743837        2/1000  2       6.16326e-09     GCF_000888115.1_ViralProj45931_genomic.fna.gz   NC_013801.1 Croton yellow vein mosaic alphasatellite, complete genome
0.743837        2/1000  2       6.16326e-09     GCF_000848385.1_ViralProj14678_genomic.fna.gz   NC_001782.1 Saccharomyces cerevisiae killer virus M1, complete genome
0.743837        2/1000  2       6.16326e-09     GCF_000004515.4_Glycine_max_v2.0_genomic.fna.gz [1191 seqs] NC_016088.2 Glycine max cultivar Williams 82 chromosome 1, Glycine_max_v2.0, whole genome shotgun   

可以看到相對上面的結果,screen功能能夠輸出更加清晰的污染檢測結果,包括了identity, shared-hashes, median-multiplicity, p-value, query-ID, query-comment等一系列的信息,更便于污染檢測。

小結

mash的使用到這里就結束了,我盡量介紹了mash最常用的功能,更多的具體用法還需要大家回到其工作手冊進行查詢(原文鏈接)。總的來說,mash最大的優(yōu)點就是速度快,使用方便,準確率高,是一款非常實用的小軟件。

?著作權歸作者所有,轉載或內容合作請聯系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

友情鏈接更多精彩內容