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

介绍

扩增子测序是分析环境微生物的常见手段，通常使用的是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.mfbz.cn/a/577514.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

go 测试和文件

go 测试和文件

go 测试和文件 需求传统测试单元测试牛刀小试总结练习文件介绍打开关闭文件读文件一次性读取文件写文件文件或文件夹是否存在文件拷贝 需求 有一个函数&#xff0c;怎样确认他运行结果是正确的&#xff1f; func addUpper(n int)int {res : 0for i : 1; i < n; i {res1}r…
阅读更多...
设计模式学习笔记 - 开源实战五（下）：总结Mybatis中用到的10种设计模式

设计模式学习笔记 - 开源实战五（下）：总结Mybatis中用到的10种设计模式

概述 本章再对 Mybatis 用到的设计模式做一个总结。它用到的设计模式也不少。有些前面章节已经经过了&#xff0c;有些则比较简单。 SqlSessionFactoryBuilder&#xff1a;为什么要用建造者模式来创建 SqlSessionFactory&#xff1f; 在《Mybatis如何权衡易用性、性能和灵活性…
阅读更多...
【算法基础实验】图论-UnionFind连通性检测之quick-find

【算法基础实验】图论-UnionFind连通性检测之quick-find

Union-Find连通性检测之quick-find 理论基础 在图论和计算机科学中&#xff0c;Union-Find 或并查集是一种用于处理一组元素分成的多个不相交集合&#xff08;即连通分量&#xff09;的情况&#xff0c;并能快速回答这组元素中任意两个元素是否在同一集合中的问题。Union-Fin…
阅读更多...
编译支持播放H265的cef控件

编译支持播放H265的cef控件

接着在上次编译的基础上增加h265支持编译支持视频播放的cef控件&#xff08;h264&#xff09; 测试页面&#xff0c;直接使用cef_enhancement,里边带着的那个html即可&#xff0c;h265视频去这个网站下载elecard,我修改的这个版本参考了里边的修改方式&#xff0c;不过我的这个…
阅读更多...
用友政务财务系统FileDownload接口存在任意文件读取漏洞

用友政务财务系统FileDownload接口存在任意文件读取漏洞

声明&#xff1a; 本文仅用于技术交流&#xff0c;请勿用于非法用途 由于传播、利用此文所提供的信息而造成的任何直接或者间接的后果及损失&#xff0c;均由使用者本人负责&#xff0c;文章作者不为此承担任何责任。 简介 用友政务财务系统是由用友软件开发的一款针对政府机…
阅读更多...
maven-idea新建和导入项目

maven-idea新建和导入项目

全局配置 新建项目 需要新建的文件夹 src/testsrc/test/javasrc/main/java 注&#xff1a;1、新建Java-class&#xff0c;输入.com.hello.hellomaven 2、快捷键psvm显示 public static void main(String[] args) {.... } package com.hello;public class hellomaven {publ…
阅读更多...
初学python记录：力扣1146. 快照数组

初学python记录：力扣1146. 快照数组

题目&#xff1a; 实现支持下列接口的「快照数组」- SnapshotArray&#xff1a; SnapshotArray(int length) - 初始化一个与指定长度相等的 类数组 的数据结构。初始时&#xff0c;每个元素都等于 0。void set(index, val) - 会将指定索引 index 处的元素设置为 val。int sna…
阅读更多...
Git泄露和hg泄露原理理解和题目实操

Git泄露和hg泄露原理理解和题目实操

一.Git泄露 1.简介 Git是一个开源的分布式版本控制系统&#xff0c;它可以实现有效控制应用版本&#xff0c;但是在一旦在代码发布的时候&#xff0c;存在不规范的操作及配置&#xff0c;就很可能将源代码泄露出去。那么&#xff0c;一旦攻击者发现这个问题之后&#xff0c;就…
阅读更多...
【算法基础实验】图论-基于DFS的连通性检测

【算法基础实验】图论-基于DFS的连通性检测

基于DFS的连通性检测 理论基础 在图论中&#xff0c;连通分量是无向图的一个重要概念&#xff0c;特别是在处理图的结构和解析图的组成时。连通分组件表示图中的一个子图&#xff0c;在这个子图中任意两个顶点都是连通的&#xff0c;即存在一条路径可以从一个顶点到达另一个顶…
阅读更多...
如何消除浏览器SmartScreen对网站“不安全”提示？

如何消除浏览器SmartScreen对网站“不安全”提示？

面对互联网时代用户对网站安全性和可信度的严苛要求&#xff0c;网站运营者时常遭遇Microsoft Defender SmartScreen&#xff08;SmartScreen&#xff09;提示网站不安全的困扰。本文将剖析SmartScreen判定网站不安全的原因&#xff0c;并为运营者提供应对策略&#xff0c;以恢…
阅读更多...
codeforce#933 题解

codeforce#933 题解

E. Rudolf and k Bridges 题意不讲了&#xff0c;不如去题干看图。 传统dp&#xff0c;每个点有两个选择&#xff0c;那么建桥要么不建。需要注意的是在状态转移的时候&#xff0c;桥是有长度的&#xff0c;如果不建需要前d格中建桥花费最少的位置作为状态转移的初态。 #incl…
阅读更多...
发那科FANUC机器人R-2000iB平衡缸维修攻略

发那科FANUC机器人R-2000iB平衡缸维修攻略

在发那科机器人中&#xff0c;平衡缸扮演着稳定机械臂运动的关键角色。它通过内部的压力调节来平衡负载&#xff0c;保证机器人的精准定位和平稳操作。一旦出现法兰克机械手平衡缸故障或损坏&#xff0c;机器人的性能可能会大打折扣&#xff0c;因此及时且正确的FANUC机械手平衡…
阅读更多...
uniapp获取当前位置及检测授权状态

uniapp获取当前位置及检测授权状态

uniapp获取当前位置及检测授权定位权限 文章目录 uniapp获取当前位置及检测授权定位权限效果图创建js文件permission.jslocation.js 使用 效果图 Android设备 点击 “设置”&#xff0c;跳转应用信息&#xff0c;打开“权限即可”&#xff1b; 创建js文件 permission.js 新建…
阅读更多...
HTTP基础知识

HTTP基础知识

1. HTTP常见的状态码有哪些&#xff1f; 常见状态码&#xff1a; 200&#xff1a;服务器已成功处理了请求。 通常&#xff0c;这表示服务器提供了请求的网页。 301 &#xff1a; (永久移动) 请求的网页已永久移动到新位置。 服务器返回此响应(对 GET 或 HEAD 请求的响应)时&a…
阅读更多...
Java使用SpringBoot和EasyExcel 实现动态数据导出实战

Java使用SpringBoot和EasyExcel 实现动态数据导出实战

Java使用SpringBoot和EasyExcel 实现动态数据导出实战 1、前言2、【资源地址】3、代码示例(demo)4、目前Java实现数据导出为Excel方式5、依赖6、总结 1、前言 工作中有用到将数据导出为Excel的场景&#xff0c;在此记录下。在日常开发中&#xff0c;Excel文件处理是一项常见的…
阅读更多...
LeetCode 面试题 08.02——迷路的机器人

LeetCode 面试题 08.02——迷路的机器人

阅读目录 1. 题目2. 解题思路3. 代码实现 1. 题目 2. 解题思路 此题就是一个典型的图搜索题&#xff0c;一种就是广度优先搜索&#xff0c;一种就是深度优先搜索。 3. 代码实现 class Solution { public:vector<vector<int>> pathWithObstacles(vector<vecto…
阅读更多...
软件需求管理规程（Word原件2024）

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

软件开发人员及用户往往容易忽略信息沟通&#xff0c;这导致软件开发出来后不能很好地满足用户的需要&#xff0c;从而造成返工。而返工不仅在技术上给开发人员带来巨大的麻烦&#xff0c;造成人力、物力的浪费&#xff0c;而且软件的性能也深受影响。所以在软件项目开发周期的…
阅读更多...
Bellman Ford算法：解决负权边图的最短路径问题

Bellman Ford算法：解决负权边图的最短路径问题

Bellman Ford算法的介绍 在计算机科学的世界中&#xff0c;Bellman Ford算法是一种解决单源最短路径问题的算法&#xff0c;它可以处理有负权边的图。这个算法的名字来源于两位科学家Richard Bellman和Lester Randolph Ford&#xff0c;他们是这个算法的发明者。 这个算法的主…
阅读更多...
hive启动beeline报错

hive启动beeline报错

问题一在zpark启动集群报错 出现上面的问题执行以下代码 chmod 777 /opt/apps/hadoop-3.2.1/logs 问题二启动beeline报错 执行 cd /opt/apps/hadoop-3.2.1 bin/hadoop dfsadmin -safemode leave 问题三执行查询语句报错 执行 set hive.exec.mode.local.autotrue;
阅读更多...
公考相丽君政治素养研习课

公考相丽君政治素养研习课

公考相丽君政治素养研习课&#xff0c;是广大公考学子提升政治素养、深化政治理解的宝贵课程。相丽君老师以其深厚的政治理论功底和丰富的教学经验&#xff0c;为学员们呈现了一堂生动而深刻的政治课。课程中&#xff0c;相老师深入浅出地讲解了政治理论的基本概念和核心思想&a…
阅读更多...
最新文章

Copyright @ 2022~2023

友情链接： 我的编程经验分享 一条长河网 尧图网站开发 尧图建网站