数据分析：扩增子-16s rRNA分析snakemake流程

news2024/10/6 8:31:29

介绍

扩增子测序是分析环境微生物的常见手段，通常使用的是16s rRNA片段。16srRNA分析主要有质控、去冗余、聚类OTU、去嵌合体、生成OTU表和物种注释等步骤。更多知识分享请到 https://zouhua.top/。

先看看前期数据处理的可视化图。

数据

18份来自宏基因组公众号的双端16s rRNA原始下机数据。

samplefq1fq2
KO1data/KO1_1.fq.gzdata/KO1_2.fq.gz
KO2data/KO2_1.fq.gzdata/KO2_2.fq.gz
KO3data/KO3_1.fq.gzdata/KO3_2.fq.gz
KO4data/KO4_1.fq.gzdata/KO4_2.fq.gz
KO5data/KO5_1.fq.gzdata/KO5_2.fq.gz
KO6data/KO6_1.fq.gzdata/KO6_2.fq.gz
OE1data/OE1_1.fq.gzdata/OE1_2.fq.gz
OE2data/OE2_1.fq.gzdata/OE2_2.fq.gz
OE3data/OE3_1.fq.gzdata/OE3_2.fq.gz
OE4data/OE4_1.fq.gzdata/OE4_2.fq.gz
OE5data/OE5_1.fq.gzdata/OE5_2.fq.gz
OE6data/OE6_1.fq.gzdata/OE6_2.fq.gz
WT1data/WT1_1.fq.gzdata/WT1_2.fq.gz
WT2data/WT2_1.fq.gzdata/WT2_2.fq.gz
WT3data/WT3_1.fq.gzdata/WT3_2.fq.gz
WT4data/WT4_1.fq.gzdata/WT4_2.fq.gz
WT5data/WT5_1.fq.gzdata/WT5_2.fq.gz
WT6data/WT6_1.fq.gzdata/WT6_2.fq.gz

步骤

  1. 质控
  2. 去冗余
  3. 聚类OTU
  4. 去嵌合体
  5. 生成OTU表
  6. 物种注释
import os
import sys
import shutil
import pandas as pd

configfile: "config.yaml"
samples = pd.read_csv(config["samples"], sep="\t", index_col=["sample"])

rule all:
    input:
        expand("{taxonomy}/taxonomy.{{ref}}.{{type}}.{res}", 
                        taxonomy=config["results"]["taxonomy"],
                        ref=["sliva", "RDP"], 
                        type=["cluster","unoise3"],
                        res=["biom","mothur","txt"])

include: "rules/00.fastqc.snk"
include: "rules/01.trim.snk"
include: "rules/02.deredundancy.snk"
include: "rules/03.clusterOTU.snk"
include: "rules/04.rechimeras.snk"
include: "rules/05.OTUtable.snk" 
include: "rules/06.assign_taxonomy.snk"
#include: "rules/07.picrust.snk"  # 这一步还没有实现
质控

质控包括剔除质量低的reads和切除带有barcodes的接头等

检查reads质量情况
  • 设置通配符函数，匹配样本名字，wildcards是snakemake的通配符函数。
  • fastqc和multiqc软件可以生成可视化的reads质量分布情况等的网页版文件。
  • config是snakemake内调用配置文件的参数。
  • fastqc和multiqc的结果可以确定质控过程reads应该保留的长度。
去冗余
  • 移除出现频率低的reads，这些reads可能是测序错误引起的，minuniquesize参数越大，相对而言计算量越小，正常建议设置成2，只去除单次出现的序列。
  • Discard_singletons排序后再去除单次出现的序列。
rule dereplicate:
    input:
        os.path.join(config["results"]["trim"], "summary_trimmed.fa")
    output:
        temp(os.path.join(config["results"]["deredundancy"], "deredundancy.fa"))
    params:
        name = config["params"]["deredundancy"]["name"],
        mins = config["params"]["deredundancy"]["mins"]
    log:
        os.path.join(config["logs"], "02.deredundancy.log")
    shell:
        '''
        vsearch --derep_fulllength {input} --output {output} --relabel {params.name} \
            --minuniquesize {params.mins} --sizeout 2>{log}
        '''

rule Discard_singletons:
    input:
        os.path.join(config["results"]["deredundancy"], "deredundancy.fa")
    output:
        os.path.join(config["results"]["deredundancy"], "sorted.fa")
    params:
        size = config["params"]["deredundancy"]["size"]
    log:
        os.path.join(config["logs"], "02.sorted.log")
    shell:
        '''
        vsearch --sortbysize {input} --output {output} --minsize {params.size} 2>{log}
        '''
聚类OTU
  • 根据序列相似性聚类（一般阈值通常为97%）能有效避免测序错误，将聚类后的序列集合称为可操作分类单元（operational taxonomic units OTUs），每一个OTU代表一类物种。常用的方法有uparse、unoise3等，但这些方法没有考虑因单碱基多态性而存在的序列多样性，于是deblur和DADA2算法被开放出来针对聚类和单碱基多态性。这里我们暂时用uparse和unoise3算法。
rule cluster_vsearch:
    input:
        os.path.join(config["results"]["deredundancy"], "sorted.fa")
    output:
        fa     = os.path.join(config["results"]["clusterOTU"], "OTU.cluster.fa"),
        biom   = os.path.join(config["results"]["clusterOTU"], "OTU.cluster.biom"),
        mothur = os.path.join(config["results"]["clusterOTU"], "OTU.cluster.mothur"),
        tab    = os.path.join(config["results"]["clusterOTU"], "OTU.cluster.txt"),
        uc     = os.path.join(config["results"]["clusterOTU"], "OTU.cluster.uc")
    params:
        identity = config["params"]["clusterOTU"]["identity"],
        threads  = config["params"]["clusterOTU"]["threads"],
        name     = config["params"]["clusterOTU"]["name"]
    log:
        os.path.join(config["logs"], "03.clusterOTU.vsearch.log")
    shell:
        '''
        vsearch --threads {params.threads} --cluster_fast {input} --id {params.identity} \
            --centroids {output.fa} --biomout {output.biom} \
            --mothur_shared_out {output.mothur} --otutabout {output.tab} \
            --relabel {params.name} --uc {output.uc}  2>{log}
        '''

rule cluster_uparse:
    input:
        os.path.join(config["results"]["deredundancy"], "sorted.fa")
    output:
        os.path.join(config["results"]["clusterOTU"], "OTU.uparse.fa")
    params:
        name = config["params"]["clusterOTU"]["name"]
    log:
        os.path.join(config["logs"], "03.clusterOTU.uparse.log")
    shell:
        '''
        usearch11 -cluster_OTUs {input} -otus {output} -relabel {params.name} 2>{log}
        '''

rule cluster_unoise3:
    input:
        os.path.join(config["results"]["clusterOTU"], "OTU.uparse.fa")
    output:
        os.path.join(config["results"]["clusterOTU"], "OTU.unoise3.fa")
    log:
        os.path.join(config["logs"], "03.clusterOTU.unoise3.log")
    shell:
        '''
        usearch11 -unoise3 {input} -zotus {output} 2>{log}
        '''
去嵌合体
  • 嵌合体的产生是因为PCR过程可能某两端DNA序列连接在一起而扩增。
  • 去嵌合体的原理是将序列比对到对应的参考数据库，16s rRNA的参考数据库主要有RDP、GreenGenes和sliva，其中常用的是RDP和sliva。
  • 去嵌合体后保留的序列即为代表序列。
rule rechimeras_sliva:
    input:
        expand("{cluster}/OTU.{{type}}.fa", 
                cluster=config["results"]["clusterOTU"], 
                type=["cluster","unoise3"])
    output:
        borderline  = os.path.join(config["results"]["rechimeras"], "chimeric.{type}.borderline"),
        nonchimeras = os.path.join(config["results"]["rechimeras"], "OTU.rechimera_silva.{type}.fa"),
        chimeras    = os.path.join(config["results"]["rechimeras"], "chimeric_silva.{type}.sequence")
    params:
        db = config["params"]["rechimeras"]["silva"] 
    log:
        os.path.join(config["logs"], "04.OTU.rechimera_silva.{type}.log")
    shell:
        '''
        vsearch --uchime_ref {input} --db {params.db} --borderline {output.borderline} \
            --chimeras {output.chimeras} --nonchimeras {output.nonchimeras} 2>{log}
        '''

rule rechimeras_RDP:
    input:
        expand("{cluster}/OTU.{{type}}.fa", 
                cluster=config["results"]["clusterOTU"], 
                type=["cluster","unoise3"])
    output:
        os.path.join(config["results"]["rechimeras"], "OTU.rechimera_RDP.{type}.fa")
    params:
        db = config["params"]["rechimeras"]["rdp"]
    log:
        os.path.join(config["logs"], "04.OTU.rechimera_RDP.{type}.log")
    shell:
        '''
        usearch11 -uchime2_ref {input} -db {params.db} -chimeras {output} -strand plus -mode balanced 2>{log}
        '''
生成OTU表
  • 聚类生成的OTU比对到去嵌合体后的OTU得到OTU表。
  • 设置比对阈值 identity: 97%。
rule make_otutab_vsearch:
    input:
        otu_ref = expand("{rechimeras}/OTU.rechimera_{{ref}}.{{type}}.fa", 
                        rechimeras=config["results"]["rechimeras"],
                        ref=["sliva", "RDP"], 
                        type=["cluster","unoise3"]),
        trim_fa = os.path.join(config["results"]["trim"], "summary_trimmed.fa")
    output:
        aln    = os.path.join(config["results"]["OTUtable"], "OTU_table.{ref}.{type}.shr.aln"),
        biom   = os.path.join(config["results"]["OTUtable"], "OTU_table.{ref}.{type}.biom"),
        mothur = os.path.join(config["results"]["OTUtable"], "OTU_table.{ref}.{type}.mothur"),
        uc     = os.path.join(config["results"]["OTUtable"], "OTU_table.{ref}.{type}.uc"),
        txt    = os.path.join(config["results"]["OTUtable"], "OTU_table.{ref}.{type}.txt")
    params:
        identity = config["params"]["OTUtable"]["identity"],
        threads  = config["params"]["OTUtable"]["threads"]
    log:
        os.path.join(config["logs"], "05.OTU_table.{ref}.{type}.log")
    shell:
        '''
        vsearch --usearch_global {input.trim_fa} --db {input.otu_ref} --id {params.identity}\
            --strand plus --threads {params.threads} --alnout {output.aln} --biomout {output.biom} \
            --mothur_shared_out {output.mothur} --uc {output.uc} --otutabout {output.txt}   2>{log}
        '''

rule OTU_tight_clusters:
    input:
        otu_ref = expand("{rechimeras}/OTU.rechimera_{{ref}}.{{type}}.fa", 
                        rechimeras=config["results"]["OTUtable"],
                        ref=["sliva", "RDP"], 
                        type=["cluster","unoise3"])
    output:
        fa = os.path.join(config["results"]["OTUtable"], "OTU.rechimera.{ref}.{type}_new.fa"),
        uc = os.path.join(config["results"]["OTUtable"], "OTU.rechimera.{ref}.{type}_hits.uc")
    params:
        identity = config["params"]["OTUtable"]["identity"],
        accepts  = config["params"]["OTUtable"]["accepts"],
        rejects  = config["params"]["OTUtable"]["rejects"]
    log:
        os.path.join(config["logs"], "05.OTU_table.{ref}.{type}.log")
    shell:
        '''
        usearch11 -cluster_fast {input.otu_ref} -id {params.identity} \
            -maxaccepts {params.accepts} -maxrejects {params.rejects} \
            -top_hit_only -uc {output.uc} -centroids {output.fa} 2>{log}
        '''
物种注释
  • 代表序列比对到16s rRNA参考数据库。
  • 控制比对identity参数。
rule assign_taxonomy:
    input:
        otu_ref = expand("{rechimeras}/OTU.rechimera_{{ref}}.{{type}}.fa", 
                        rechimeras=config["results"]["rechimeras"],
                        ref=["sliva", "RDP"], 
                        type=["cluster","unoise3"])
    output:
        biom   = os.path.join(config["results"]["taxonomy"], "taxonomy.{ref}.{type}.biom"),
        mothur = os.path.join(config["results"]["taxonomy"], "taxonomy.{ref}.{type}.mothur"),
        txt    = os.path.join(config["results"]["taxonomy"], "taxonomy.{ref}.{type}.txt")
    params:
        identity = config["params"]["taxonomy"]["identity"],
        threads  = config["params"]["taxonomy"]["threads"],
        database = config["params"]["taxonomy"]["silva"]
    log:
        os.path.join(config["logs"], "06.taxonomy.{ref}.{type}.log")
    shell:
        '''
        vsearch --usearch_global {input.otu_ref} --db {params.database} \
                --biomout {output.biom} --mothur_shared_out {output.mothur} --otutabout {output.txt} \
                --id {params.identity} --threads {params.threads} 2>{log}
        '''

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/1627491.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

C# WinForm —— 08 Form初始化、布局、注册事件

C# WinForm —— 08 Form初始化、布局、注册事件

Form 初始化 Form初始化的时候会调用 Designer.cs 里的 InitializeComponent(); 函数，在InitializeComponent(); 函数里面有Load Form语句时会调用 FrmLogin_Load()函数 Form布局 两种方式： 拖控件到窗体，设置属性在Load事件中写代码添加…
阅读更多...
线性神经网络示例

线性神经网络示例

通过5个条件判定一件事情是否会发生，5个条件对这件事情是否发生的影响力不同，计算每个条件对这件事情发生的影响力多大，写一个线性神经网络模型pytorch程序,最后打印5个条件分别的影响力。 一 在这个场景中，一个线性神经网络&…
阅读更多...
knife4j swagger 使用笔记

knife4j swagger 使用笔记

1.接口访问的端口跟后台设置的不一致，接口请求无反应 处理办法 2.响应参数不显示问题 （1）返回的参数里面一定要有响应的参数对象，如下： （2）TableDataInfo 定义成泛型类 TableDataInfo package…
阅读更多...
Int4：Lucene 中的更多标量量化

Int4：Lucene 中的更多标量量化

作者：来自 Elastic Benjamin Trent, Thomas Veasey 在 Lucene 中引入 Int4 量化 在之前的博客中，我们全面介绍了 Lucene 中标量量化的实现。 我们还探索了两种具体的量化优化。 现在我们遇到了一个问题：int4 量化在 Lucene 中是如何工作的以…
阅读更多...
软件需求管理规程（Word原件2024）

软件需求管理规程（Word原件2024）

软件开发人员及用户往往容易忽略信息沟通，这导致软件开发出来后不能很好地满足用户的需要，从而造成返工。而返工不仅在技术上给开发人员带来巨大的麻烦，造成人力、物力的浪费，而且软件的性能也深受影响。所以在软件项目开发周期的…
阅读更多...
单片机为什么有多组VDD？

单片机为什么有多组VDD？

以前我在画尺寸小的PCB时，比较头痛，特别是芯片引脚又多的，芯片底下，又不能打太多过孔。 可能有些老铁也比较好奇，为什么一个单片机芯片，有这么多组VDD和VSS。 比如下面这个100个引脚的STM32单片机。 有5组…
阅读更多...
前端实现将当前页面内容下载成图片（图片可做到高清画质）

前端实现将当前页面内容下载成图片（图片可做到高清画质）

插件背景： html2canvas可以把你想要转变的元素变为图片，使用file-saver下载图片。 1、安装html2canvas、file-saver npm install html2canvasnpm install file-saver --save 2、在Vue组件中引入并使用html2canvas、file-saver import html2canvas fro…
阅读更多...
智慧旅游开启智慧出行新时代，科技引领旅行新风尚：以科技为引领，推动旅游业智慧化升级，为旅行者提供更加便捷、高效的旅行服务

智慧旅游开启智慧出行新时代，科技引领旅行新风尚：以科技为引领，推动旅游业智慧化升级，为旅行者提供更加便捷、高效的旅行服务

一、引言 随着信息技术的飞速发展，智慧旅游作为一种全新的旅游形态，正逐渐改变着人们的出行方式。它利用现代科技手段，实现旅游资源的智能化管理、旅游信息的智能化传播和旅游服务的智能化提供，为旅行者带来更加便捷、高效的旅行…
阅读更多...
Qt下使用OpenCV截取图像并在QtableWidget表格上显示

Qt下使用OpenCV截取图像并在QtableWidget表格上显示

文章目录 前言一、在QLabel上显示图片并绘制矩形框二、保存矩形框数据为CSV文件三、保存截取图像四、将截取图像填充到表格五、图形视图框架显示图像六、示例完整代码总结 前言 本文主要讲述了在Qt下使用OpenCV截取绘制的矩形框图像，并将矩形框数据保存为CSV文件&a…
阅读更多...
气膜仓库：现代化仓储新选择—轻空间

气膜仓库：现代化仓储新选择—轻空间

气膜仓库，作为现代化仓储的新选择，越来越受到人们的青睐。相比传统料仓，气膜仓库具有诸多优势，使其成为各行各业的首选储存解决方案。 1. 高效节能 气膜仓库的建设周期短，基础简单，安装快捷，能耗…
阅读更多...
C#命名空间常用函数

C#命名空间常用函数

在C#中，不同命名空间下有各种常用函数，下面列举一些常见的函数及其对应的命名空间： System命名空间： Console.WriteLine()：用于向控制台输出信息。Convert.ToInt32()：用于将其他数据类型转换为整数类型。 S…
阅读更多...
Kafka 3.x.x 入门到精通（05）——对标尚硅谷Kafka教程

Kafka 3.x.x 入门到精通（05）——对标尚硅谷Kafka教程

Kafka 3.x.x 入门到精通（05）——对标尚硅谷Kafka教程 2. Kafka基础2.1 集群部署2.2 集群启动2.3 创建主题2.4 生产消息2.5 存储消息2.6 消费消息2.6.1 消费消息的基本步骤2.6.2 消费消息的基本代码2.6.3 消费消息的基本原理2.6.3.1消费者组2.6.3.1.1 消费…
阅读更多...
凹凸技术揭秘·羚珑智能设计平台·逐梦设计数智化

凹凸技术揭秘·羚珑智能设计平台·逐梦设计数智化

从技术和功能形态层面，我们把设计数智化分成了两个方向，一个方向是「模板化设计」，另一个方向是「程序化设计」。 2、模板化设计— 「模板化设计」的核心目标：是实现线下设计物料的数字化，在数字化设计资产的基础之上…
阅读更多...
WildCard开通GitHub Copilot

WildCard开通GitHub Copilot

更多AI内容请关注我的专栏：《体验AI》 期待您的点赞👍收藏⭐评论✍ WildCard开通GitHub Copilot GitHub Copilot 简介主要功能工作原理 开通过程1、注册Github账号2、准备一张信用卡或虚拟卡3、进入github copilot页4、选择试用5、选择支付方式6、填写卡…
阅读更多...
C语言：插入排序

C语言：插入排序

插入排序 1.解释2.步骤3.举例分析示例结果分析 1.解释 插入排序是一种简单直观的排序算法，它的工作原理是通过构建有序序列，对于未排序数据，在已排序序列中从后向前扫描，找到相应位置并插入。插入排序在实现上，通常采…
阅读更多...
SSH新功能揭秘：远程工作提升指南【AI写作】

SSH新功能揭秘：远程工作提升指南【AI写作】

首先，这篇文章是基于笔尖AI写作进行文章创作的，喜欢的宝子，也可以去体验下，解放双手，上班直接摸鱼~ 按照惯例，先介绍下这款笔尖AI写作，宝子也可以直接下滑跳过看正文~ 笔尖Ai写作：…
阅读更多...
如何实现直播声卡反向给手机充电功能呢？

如何实现直播声卡反向给手机充电功能呢？

在数字化时代的浪潮中，声卡作为多媒体系统的核心组件，扮演着声波与数字信号相互转换的关键角色。它不仅能够将来自各类音源的原始声音信号转换为数字信号，进而输出到各类声响设备，更能够通过音乐设备数字接口(MIDI)发出合成乐器的…
阅读更多...
多家企业机密数据遭Lockbit3.0窃取，亚信安全发布《勒索家族和勒索事件监控报告》

多家企业机密数据遭Lockbit3.0窃取，亚信安全发布《勒索家族和勒索事件监控报告》

本周态势快速感知 本周全球共监测到勒索事件87起，与上周相比勒索事件大幅下降。美国依旧为受勒索攻击最严重的国家，占比45%。 本周Cactus是影响最严重的勒索家族，Lockbit3.0和Bianlian恶意家族紧随其后，从整体上看Lockbit3.0依旧…
阅读更多...
BERT-CRF 微调中文 NER 模型

BERT-CRF 微调中文 NER 模型

文章目录 数据集模型定义数据集预处理BIO 标签转换自定义Dataset拆分训练、测试集 训练验证、测试指标计算推理其它相关参数CRF 模块 数据集 CLUE-NER数据集：https://github.com/CLUEbenchmark/CLUENER2020/blob/master/pytorch_version/README.md 模型定义 imp…
阅读更多...
如何安全进行速卖通自养号测评操作？

如何安全进行速卖通自养号测评操作？

对于新加入的卖家而言，进行销量测评显得尤为关键。速卖通平台上的新店往往难以获得活动的扶持，且初始流量相当有限。因此，开店的首要任务便是积极展开测评工作，努力积累初始的评论和销售记录。测评的益处颇为显著，它不…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023