宏基因组实战5. sourmash基于Kmer比较数据集

前情提要如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章宏基因组分析理论教程微生物组入门圣经+宏基因组分析实操课程1背景知识-Shell入门与本地blast实战2...

前情提要

如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章

测试数据

刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：

下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。

最后送上本教程使用到的所有文件同步共享文件夹链接：http://pan.baidu.com/s/1hsIjosk 密码：y0tb 。

sourmash比较数据集

https://2017-cicese-metagenomics.readthedocs.io/en/latest/sourmash_compare.html

sourmash教程

主页：https://sourmash.readthedocs.io/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次。

目的

比对reads至拼装结果
比较数据集并构建聚类图

采用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-mer是很不同的
长k-mer具有很强的物种特异性

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://osf.io/k3sq7/download
# 此文件很大，如果继续以之前的分析，可以执行下行链接至此目录
# ln ../data/SRR1976948.abundtrim.subset.pe.fq.gz ./
curl -L -o SRR1976948.megahit.abundtrim.subset.pe.assembly.fa https://osf.io/dxme4/download
curl -L -o SRR1976948.spades.abundtrim.subset.pe.assembly.fa https://osf.io/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

生成Hu_metagenomes.dendro.png和Hu_metagenomes.matrix.png两个文件，分别是树形图，还有热图+树形图，可以用Filezilla下载查看。

图5. 样品间kmer聚类结果

此软件还可以做什么？

更多的kmer研究
物种分类
功能分析

想学习看该软件的官方帮助文档吧。https://sourmash.readthedocs.io/en/latest/

宏基因组实战5. sourmash基于Kmer比较数据集

你可能感兴趣的文章

相关问题

0 条评论

作家榜 »