前情提要如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章 sourmash比较数据集https://2017-cicese-metagenomics./en/latest/sourmash_compare.html sourmash教程主页：https://sourmash./en/latest/ Mash: fast genome and metagenome distance estimation using MinHash. Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, Phillippy AM. Genome Biol. 2016 Jun 20;17(1):132. doi: 10.1186/s13059-016-0997-x. sourmash是一款超快、轻量级核酸搜索和比对软件，方法叫MinHash。这这一个命令行使用的python包，可以从DNA序列中快速分析kmer，并进行样品比较和绘图。目的是快速且准确的比较大数据集。在2016年6月20日发表，目前已经被引用61次。 目的
采用k-mer Jaccard distance进行样品比较什么是k-mer k-mer就是长度为k的DNA序列 ATTG - a 4-mer
ATGGAC - a 6-mer比如，一个长度为16的序列可以提取11个长度为6的 k-mers AGGATGAGACAGATAG
AGGATG
GGATGA
GATGAG
ATGAGA
TGAGAC
GAGACA
AGACAG
GACAGA
ACAGAT
CAGATA
AGATAG需要记住的概念：
K-mers与组装图三大基因组拼接方法之一，就是将基因组打断成k-mer再拼接，通过k-mer步移实现拼接 为什么采用k-mers而不是全长序列拼接 计算机喜欢k-mer，因为匹配准确快速。 长k-mers存在物种特异性图1. 以k=40为例，kmer很容易按属水平分开细菌 采用k-mers比较样品安装sourmash # 设置工作目录，这是我的目录，学习时请修改为你工作的目录
pwd=~/test/metagenome17
cd ${pwd}
# 依赖关系，之前安装成功可跳过
sudo apt-get -y update && sudo apt-get install -y python3.5-dev python3.5-venv make libc6-dev g++ zlib1g-dev
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer
# 安装程序包
pip install -U https://github.com/dib-lab/sourmash/archive/master.zip产生Illumina reads的signature mkdir work
cd work
curl -L -o SRR1976948.abundtrim.subset.pe.fq.gz https:///k3sq7/download
# 此文件很大，如果继续以之前的分析，可以执行下行链接至此目录
# ln ../data/SRR1976948.abundtrim.subset.pe.fq.gz ./
curl -L -o SRR1976948.megahit.abundtrim.subset.pe.assembly.fa https:///dxme4/download
curl -L -o SRR1976948.spades.abundtrim.subset.pe.assembly.fa https:///zmekx/download图2. sourmash工作流程 mkdir ../sourmash
cd ../sourmash
# 计算过滤序列的k51特征
sourmash compute -k51 --scaled 10000 ../work/SRR1976948.abundtrim.subset.pe.fq.gz -o SRR1976948.reads.scaled10k.k51.sig比较序列与组装结果 sourmash compute -k51 --scaled 10000 ../work/SRR1976948.spades.abundtrim.subset.pe.assembly.fa -o SRR1976948.spades.scaled10k.k51.sig
sourmash compute -k51 --scaled 10000 ../work/SRR1976948.megahit.abundtrim.subset.pe.assembly.fa -o SRR1976948.megahit.scaled10k.k51.sig图3. 评估污染比例 sourmash search SRR1976948.reads.scaled10k.k51.sig SRR1976948.megahit.scaled10k.k51.sig --containment
sourmash search SRR1976948.reads.scaled10k.k51.sig SRR1976948.spades.scaled10k.k51.sig --containment结果如下： loaded query: SRR1976948.abundtrim.subset.pe... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
48.7% SRR1976948.megahit.abundtrim.subset.pe.assembly.fa
loaded query: SRR1976948.abundtrim.subset.pe... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
47.5% SRR1976948.spades.abundtrim.subset.pe.assembly.fa为什么只有不到一半的Kmer在拼接结果中？ 我们看看拼接结果中有多少kmer在原始序列中 sourmash search SRR1976948.megahit.scaled10k.k51.sig SRR1976948.reads.scaled10k.k51.sig --containment
sourmash search SRR1976948.spades.scaled10k.k51.sig SRR1976948.reads.scaled10k.k51.sig --containment结果如下： loaded query: SRR1976948.megahit.abundtrim.s... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
99.8% SRR1976948.abundtrim.subset.pe.fq.gz
loaded query: SRR1976948.spades.abundtrim.su... (k=51, DNA)
loaded 1 signatures and 0 databases total.
1 matches:
similarity match
---------- -----
99.9% SRR1976948.abundtrim.subset.pe.fq.gz基本都全部可以找到。这是因为原始序列中包含大量illumina测序的错误，而在拼接中被丢弃或校正。 特征比较 为了更深刻理解，我们采用osf下载更多数据集比较 pip install osfclient
osf -p ay94c fetch osfstorage/signatures/SRR1977249.megahit.scaled10k.k51.sig SRR1977249.megahit.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977249.reads.scaled10k.k51.sig SRR1977249.reads.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977249.spades.scaled10k.k51.sig SRR1977249.spades.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.megahit.scaled10k.k51.sig SRR1977296.megahit.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.reads.scaled10k.k51.sig SRR1977296.reads.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/SRR1977296.spades.scaled10k.k51.sig SRR1977296.spades.scaled10k.k51.sig
osf -p ay94c fetch osfstorage/signatures/subset_assembly.megahit.scaled10k.k51.sig subset_assembly.megahit.scaled10k.k51.sig图4. 样品间比较原理 # 比较所有signature文件
sourmash compare *sig -o Hu_metagenomes
# 比较结果绘图
sourmash plot --labels Hu_metagenomes生成 图5. 样品间kmer聚类结果 此软件还可以做什么？
想学习看该软件的官方帮助文档吧。https://sourmash./en/latest/ |
|
|
