使用 snakemake 对端 bwa 比对

Question

我刚开始使用snakemake，在snakemake 中做映射时遇到了一个简单的问题。我有一对 _1.fastq.gz 和 _2.fastq.gz，我想为大约 10 对 fastq.gz 做对端映射。所以我为此写了一个snakemake文件。

import os
import snakemake.io
import glob

(SAMPLES,READS,) = glob_wildcards("raw/{sample}_{read}.fastq.gz")
READS=["1","2"]
REF="/data/data/reference/refs/ucsc.hg19.fasta.gz"

rule all:
    input: expand("raw/{sample}.bam",sample=SAMPLES)

rule bwa_map:
    input:
        ref=REF,
        r1=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS),
        r2=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS)

    output: "raw/{sample}.bam"

    shell: "bwa mem -M -t 8 {input.ref} {input.r1} {input.r2} | samtools view -Sbh - > {output}"

错误：

RuleException:
CalledProcessError in line 17 of /data/data/Samples/snakemake-example/WGS-test/step2.smk:
Command ' set -euo pipefail;  bwa mem -M -t 8 /data/data/reference/refs/ucsc.hg19.fasta.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz raw/sub1_1.fastq.gz raw/sub1_2.fastq.gz raw/sub2_1.fastq.gz raw/sub2_2.fastq.gz | samtools view -Sbh - > raw/sub2.bam ' returned non-zero exit status 1.
  File "/data/data/Samples/snakemake-example/WGS-test/step2.smk", line 17, in __rule_bwa_map
  File "/root/miniconda3/envs/bioinfo/lib/python3.6/concurrent/futures/thread.py", line 56, in run

我想要的输出就像生成所有 10 个 bam 文件一样

sub1.bam sub2.bam sub3.bam ...

好像把所有的fastq文件放到一个命令里了。我如何在不使用硬编码方法的情况下自动将它们和运行成对分开。请指教

Answer 1

第一条规则（此处 rule all）指定您希望在 snakemake 工作流程中创建的文件。

对于给定的文件f，在rule all::input中，snakemake会遍历所有的规则并尝试找到一个可以创建f的（基于模式匹配的output 每条规则的片段）。

假设f = raw/my_sample.bam

一旦 snakemake 找到可以创建 f 的规则，它将确定创建该文件所需的所有输入文件。

所以在这里，snakemake 发现 f = raw/my_sample.bam 可以通过 rule bwa_map 创建（因为 f 匹配模式 raw/<anything>.bam）然后确定需要制作哪些文件f 基于 input 段。

Snakemake 认为：如果我有，我可以做出 raw/my_sample.bam 文件 ref="/data/data/reference/refs/ucsc.hg19.fasta.gz" 文件 r1=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS) 和文件 r2=expand("raw/{sample}_{read}.fastq.gz",sample=SAMPLES,read=READS)

在 expand 中，r1 将 sample 扩展到 SAMPLES 中的每个值，将 read 扩展到 READS 中的每个值。但是您在 SAMPLES 中有 10 个值，在 READS 中有 2 个值，因此 r1 为它尝试制作的每个输出文件扩展到 20 个不同的文件路径。它忽略 output 子句中存在的 sample 通配符（因为您已在 expand 调用中覆盖它）。

您必须让输出子句中定义的通配符冒泡到输入子句:

import os
import snakemake.io
import glob

(SAMPLES,READS,) = glob_wildcards("raw/{sample}_{read}.fastq.gz")
READS=["1","2"]
REF="/data/data/reference/refs/ucsc.hg19.fasta.gz"

rule all:
    input: expand("raw/{sample}.bam",sample=SAMPLES)

rule bwa_map:
    input:
        ref=REF,
        # determine `r1` based on the {sample} wildcard defined in `output`
        # and the fixed value `1` to indicate the read direction
        r1="raw/{sample}_1.fastq.gz",
        # determine `r2` based on the {sample} wildcard similarly
        r2="raw/{sample}_2.fastq.gz"

    output: "raw/{sample}.bam"

    # better to pass in the threads than to hardcode them in the shell command
    threads: 8

    shell: "bwa mem -M -t {threads} {input.ref} {input.r1} {input.r2} | samtools view -Sbh - > {output}"

我建议你看看 bwa 对齐规则是如何写在 snakemake wrappers 资源中的（你也可以考虑使用 wrapper）：https://snakemake-wrappers.readthedocs.io/en/stable/wrappers/bwa/mem.html

Off-topic：从代码审查的角度来看，我质疑为什么要将对齐的数据输出到 raw 目录？将对齐的数据输出到 align 或 aligned 会更有意义吗？您似乎也在从不使用的包中导入。

使用 snakemake 对端 bwa 比对

use snakemake pair-end bwa alignment

python

bioinformatics

snakeyaml

snakemake