接受略有不同的 Snakemake 规则输入(.fq 与 .fq.gz)

2024-04-17

我是 Snakemake 的新手,希望能够选择一对.fq文件或一对.fq.gz文件并运行它们trim_galore得到一对修剪过的.fq.gz输出文件。在不提供所有 Snakefile 的情况下,我得到了下面丑陋的解决方案,我只是复制了规则并更改了输入。什么是更好的解决方案?

#Trim galore paired end trimming rule for unzipped fastqs:
rule trim_galore_unzipped_PE:
    input:
        r1=join(config['fq_in_path'], '{sample}1.fq'),
        r2=join(config['fq_in_path'], '{sample}2.fq'),
    output:
        r1=join(config['trim_out_path'], '{sample}1_val_1.fq.gz'),
        r2=join(config['trim_out_path'], '{sample}2_val_2.fq.gz'),
    params:
        out_path=config['trim_out_path'],
    conda:
        'envs/biotools.yaml',
    shell:
        'trim_galore --gzip -o {params.out_path} --paired {input.r1} {input.r2}'

#Trim galore paired end trimming rule for gzipped fastqs:
rule trim_galore_zipped_PE:
    input:
        r1=join(config['fq_in_path'], '{sample}1.fq.gz'),
        r2=join(config['fq_in_path'], '{sample}2.fq.gz'),
    output:
        r1=join(config['trim_out_path'], '{sample}1_val_1.fq.gz'),
        r2=join(config['trim_out_path'], '{sample}2_val_2.fq.gz'),
    params:
        out_path=config['trim_out_path'],
    conda:
        'envs/biotools.yaml',
    shell: 
        'trim_galore --gzip -o {params.out_path} --paired {input.r1} {input.r2}'

使用输入函数可能是最好的解决方案,如下所示:

  1. 将通配符传递给输入函数
  2. 使用已知的 YAML 值,使用该示例名称构建理论文件名。
  3. 使用 python 函数检查哪个文件(技术上是文件后缀)有效
  4. 建立有效文件列表
  5. 返回并解压有效文件列表。

Notes:

  • 输入和输出应该具有相同的通配符,如果没有,则会导致问题
  • 在输入函数中,确保它不能返回空字符串,因为 Snakemake 将此解释为“缺少输入”要求,这不是您想要的。
  • 如果您采纳这些建议,请更新规则名称,我忘了。

蛇文件:

 configfile: "config.yaml"

 from os.path import join
 from os.path import exists

 rule all:
     input:
         expand("{trim_out_path}/{sample}.{readDirection}.fq.gz",
             trim_out_path=config["trim_out_path"],
             sample=config["sampleList"],
             readDirection=['1','2'])


 def trim_galore_input_determination(wildcards):
     potential_file_path_list = []
     # Cycle through both suffix possibilities:
     for fastqSuffix in [".fq", ".fq.gz"]:

         # Cycle through both read directions
         for readDirection in ['.1','.2']:

             #Build the list for ech suffix
             potential_file_path = config["fq_in_path"] + "/" + wildcards.sample + readDirection + fastqSuffix

             #Check if this file actually exists
             if exists(potential_file_path):

                 #If file is legit, add to list of acceptable files
                 potential_file_path_list.append(potential_file_path)

     # Checking for an empty list
     if len(potential_file_path_list):
         return potential_file_path_list
     else:
         return ["trim_galore_input_determination_FAILURE" + wildcards.sample]

 rule trim_galore_unzipped_PE:
     input:
         unpack(trim_galore_input_determination)
     output:
         expand("{trim_out_path}/{{sample}}.{readDirection}.fq.gz",
             trim_out_path=config["trim_out_path"],
             readDirection=['1','2'])
     params:
         out_path=config['trim_out_path'],
     conda:
         'envs/biotools.yaml',
     shell:
         'trim_galore --gzip -o {params.out_path} --paired {input}'

配置.yaml:

fq_in_path: input/fq
trim_out_path: output
sampleList: ["mySample1", "mySample2"]

$tree:

|-- [tboyarsk      1540 Sep  6 15:17]  Snakefile
|-- [tboyarsk        82 Sep  6 15:17]  config.yaml
|-- [tboyarsk       512 Sep  6  8:55]  input
|   |-- [tboyarsk       512 Sep  6  8:33]  fq
|   |   |-- [tboyarsk         0 Sep  6  7:50]  mySample1.1.fq
|   |   |-- [tboyarsk         0 Sep  6  8:24]  mySample1.2.fq
|   |   |-- [tboyarsk         0 Sep  6  7:50]  mySample2.1.fq
|   |   `-- [tboyarsk         0 Sep  6  8:24]  mySample2.2.fq
|   `-- [tboyarsk       512 Sep  6  8:55]  fqgz
|       |-- [tboyarsk         0 Sep  6  7:50]  mySample1.1.fq.gz
|       |-- [tboyarsk         0 Sep  6  8:32]  mySample1.2.fq.gz
|       |-- [tboyarsk         0 Sep  6  8:33]  mySample2.1.fq.gz
|       `-- [tboyarsk         0 Sep  6  8:32]  mySample2.2.fq.gz
`-- [tboyarsk       512 Sep  6  7:55]  output

$snakemake -dry(输入:fg)

 rule trim_galore_unzipped_PE:
     input: input/fq/mySample1.1.fq, input/fq/mySample1.2.fq
     output: output/mySample1.1.fq.gz, output/mySample1.2.fq.gz
     jobid: 1
     wildcards: sample=mySample1


 rule trim_galore_unzipped_PE:
     input: input/fq/mySample2.1.fq, input/fq/mySample2.2.fq
     output: output/mySample2.1.fq.gz, output/mySample2.2.fq.gz
     jobid: 2
     wildcards: sample=mySample2


 localrule all:
     input: output/mySample1.1.fq.gz, output/mySample2.1.fq.gz, output/mySample1.2.fq.gz, output/   mySample2.2.fq.gz
     jobid: 0

 Job counts:
         count   jobs
         1       all
         2       trim_galore_unzipped_PE
         3

$snakemake -dry(输入:fgqz)

 rule trim_galore_unzipped_PE:
     input: input/fqgz/mySample1.1.fq.gz, input/fqgz/mySample1.2.fq.gz
     output: output/mySample1.1.fq.gz, output/mySample1.2.fq.gz
     jobid: 1
     wildcards: sample=mySample1


 rule trim_galore_unzipped_PE:
     input: input/fqgz/mySample2.1.fq.gz, input/fqgz/mySample2.2.fq.gz
     output: output/mySample2.1.fq.gz, output/mySample2.2.fq.gz
     jobid: 2
     wildcards: sample=mySample2


 localrule all:
     input: output/mySample1.1.fq.gz, output/mySample1.2.fq.gz, output/mySample2.1.fq.gz, output/   mySample2.2.fq.gz
     jobid: 0

 Job counts:
         count   jobs
         1       all
         2       trim_galore_unzipped_PE
         3

有多种方法可以使其更加通用,但由于您声明并使用 YAML 配置来构建大部分文件名,因此我将避免在答案中讨论它。只是说这是可能的并且有点令人鼓舞。

“--paired {input}”将扩展以提供这两个文件。由于 for 循环,1 总是在 2 之前。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

接受略有不同的 Snakemake 规则输入(.fq 与 .fq.gz) 的相关文章

  • NCBI Genbank核苷酸序列数据库检索基因序列解读

    核酸数据库 Genbank数据库 Nucleotide数据库 一 基因序列注释内容解析 以dut基因编码的大肠杆菌酶dutpase为例 在Nucleotide数据库search X01714或者dutpase 检索链接https www n
  • BDS - Chapter - 3 - Remedial Unix Shell

    This book assumes you re familiar with basic topics such as what a terminal is what the shell is the Unix filesystem hie
  • Python 中的递归生成器

    我编写了一个函数来返回一个生成器 其中包含给定长度的子字符串的每个唯一组合 这些子字符串包含主字符串中的 n 个以上元素 举例来说 如果我有 abcdefghi 和长度为 2 的探针 并且每个列表有 4 个元素的阈值 我想得到 ab cd
  • fasta.gz 上的 SeqIO.parse

    编码新手 Pytho biopython 新手 这是我在网上的第一个问题 如何打开压缩的 fasta gz 文件以提取信息并在我的函数中执行计算 这是我正在尝试执行的操作 我尝试了不同的方法 以及错误是什么的简化示例 我正在使用的 gzip
  • 当集群(slurm)取消作业时 Snakemake 挂起

    也许答案对很多人来说都是显而易见的 但我很惊讶我找不到关于这个主题的问题 这对我来说是一个主要问题 我将不胜感激的提示 当在 slurm 管理的集群上提交作业时 如果队列管理器取消该作业 例如 由于资源或时间不足 snakemake 似乎不
  • 查询 DNS 服务记录以查找主机名和 TCP/IP

    在一篇关于生命科学标识符 see LSID Tester 用于测试生命科学标识符解析服务的工具 罗德里克 DM 佩奇博士写道 给定 LSID urn lsid ubio org namebank 11815 在 DNS 中查询 SRV 记录
  • 在 Snakemake HTML 报告中包含参数和源代码

    我想在我的html报告中包含shell命令以及snakemake规则的外部脚本的源代码 我看到人们在RULE序列的表中包含这些 下面的示例是文档中基本示例的一部分 https snakemake readthedocs io en stab
  • 读取 .fasta 序列以提取核苷酸数据,然后写入 TabDelimited 文件

    在继续之前 我想请读者参考我之前使用 Perl 时遇到的问题 因为我是这一切的初学者 以下是我这几天发的帖子 按时间顺序排列 如何平均制表符分隔数据中的列值 Solved 为什么我在输出文件中看不到计算结果 Solved 使用 fasta
  • R-Shiny 中的自动多文件下载

    我正在尝试弄清楚如何获得data frame对其自身进行子集化 然后为每个子集编写一个 csv 文件 我正在写一个shiny应用程序将为不同的仪器生成模板文件 我需要能够为每个批次 板 任何内容获取一个文件 显然 我们可以进行手动排序 但这
  • “通配符”对象没有属性“输出”

    我收到一条相当简单的规则的错误 我必须为另一个程序编写一个任务文件 需要一个 tsv 文件 我从配置文件中读取一定数量的参数 并使用 shell 命令将它们写入文件中 Code rule create tasks output temp t
  • 在 Snakemake 脚本中使用 argparse

    是否可以将自定义命令行参数传递给snakemake脚本 我已经尝试过了 但是用以下命令执行 Snakefileargparse结果出错snakemake error unrecognized arguments zz 下面是一个示例脚本 i
  • 根据同一列表中的下一个项目从列表中删除项目

    我刚刚开始学习 python 这里有一个蛋白质序列的排序列表 总共 59 000 个序列 其中一些是重叠的 我在这里列出了一个玩具清单 例如 ABCDE ABCDEFG ABCDEFGH ABCDEFGHIJKLMNO CEST DBTSF
  • awk 命令在 Snakemake --use-singularity 中失败

    我正在尝试将 Snakemake 与 Singularity 结合起来 我注意到一个简单的awk使用奇点时命令不再起作用 这 1最后一行被 bash 替换 而不是被用作第一个字段awk 这是一个最小的工作示例 蛇形锉刀 singularit
  • 用于生物信息学/生物统计学/医学研究的 Clojure 或 Scala [关闭]

    Closed 这个问题是基于意见的 help closed questions 目前不接受答案 我不是一个专业的程序员 我的领域是医学研究 但我对C C 和各种脚本语言相当有能力 不久前我对 Lisp 很感兴趣 但一直没有时间认真学习它 短
  • Snakemake 通配符:使用目录输出中的通配符文件

    我是 Snakemake 的新手 并尝试在规则中使用特定文件 来自directory 克隆 git 存储库的另一个规则的输出 目前 这给了我一个错误Wildcards in input files cannot be determined
  • 如何快速识别 Snakemake 中的规则是否需要输入函数

    我正在关注其文档页面上的 Snakemake 教程 并且确实陷入了输入函数的概念https snakemake readthedocs io en stable tutorial advanced html step 3 input fun
  • Snakemake - 如何使用文件的每一行作为输入?

    我需要使用文件的每一行tissuesused txt作为snakemake中并行规则的输入 我想总共大约有 48 个工作机会 for line in cat tissuesused txt do echo Sorting line phen
  • 当我从网络运行 CGI 脚本时,为什么 python 找不到某些模块?

    我不知道这里可能有什么问题 我有一些来自 Biopython 的模块 当使用交互式提示或通过命令行执行 python 脚本时 我可以轻松导入这些模块 问题是 当我尝试在 Web 可执行 cgi 脚本中导入相同的 biopython 模块时
  • Snakemake 声明规则以非零退出代码退出,即使使用“|| true”?

    我的 Snakemake 管道断言 每当我运行任何规则时 我的代码都会引发非零退出代码 即使我的代码在我手动运行相同的代码时返回错误代码 0 并且在 Snakemake 中运行时它可以正常工作 根据建议这个问题 https stackove
  • 将文本文件转换为 plink PED 和 MAP 格式

    我有以下数据 其中的一小部分 名为 short2 pre snp tumor txt rs987435 C G 1 1 1 0 2 rs345783 C G 0 0 1 0 0 rs955894 G T 1 1 2 2 1 rs608879

随机推荐