使用蛇形对端 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 模式）然后根据input 段确定制作f 需要哪些文件.

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 包装器资源中的（您也可以考虑使用包装器）：https://snakemake-wrappers.readthedocs.io/en/stable/wrappers/bwa/mem.html

题外话：从代码审查的角度来看，我质疑为什么要将对齐的数据输出到 raw 目录？将对齐的数据输出到align 或aligned 是否更有意义？您似乎还从不使用的包中导入。