2

我正在尝试一个 SnakeMake 管道,但我遇到了一个我真的不明白的错误。

我有一个目录(raw_data),其中有输入文件:

ll /home/nico/labo/etudes/Optimal/data/raw_data
total 41M
drwxrwxr-x  2 nico nico 4,0K mars   6 16:09 ./
drwxrwxr-x 11 nico nico 4,0K mars   6 16:14 ../
-rw-rw-r--  1 nico nico  15M févr. 27 12:21 sampleA_R1.fastq.gz
-rw-rw-r--  1 nico nico  19M févr. 27 12:22 sampleA_R2.fastq.gz
-rw-rw-r--  1 nico nico 3,4M févr. 27 12:21 sampleB_R1.fastq.gz
-rw-rw-r--  1 nico nico 4,3M févr. 27 12:22 sampleB_R2.fastq.gz

此目录包含 2 个样本的 4 个文件。我为 SnakeMake 管道创建了一个配置 json 文件,名为config_snakemake_Optimal_mapping_BaL.json

{
    "fastqExtension": "fastq.gz",
    "fastqDir": "/home/nico/labo/etudes/Optimal/data/raw_data",
    "outputDir": "/home/nico/labo/etudes/Optimal/data/mapping_BaL",
    "logDir": "logs",
    "reference": {
        "fasta": "/home/nico/labo/references/genomes/HIV1/BaL_AY713409/BaL_AY713409.fasta",
        "index": "/home/nico/labo/references/genomes/HIV1/BaL_AY713409/BaL_AY713409.fasta.bwt"
    }
}

最后是 SnakeMake 文件snakefile_bwa_samtools.py

import subprocess
from os.path import join

### Globals ---------------------------------------------------------------------

# A Snakemake regular expression matching fastq files.

SAMPLES, = glob_wildcards(join(config["fastqDir"], "{sample}_R1."+config["fastqExtension"]))
print(SAMPLES)

### Rules -----------------------------------------------------------------------

# Pipeline output files
rule all:
    input: expand(join(config["outputDir"], "{sample}.bam.bai"), sample=SAMPLES)

# Reads alignment on reference genome and BAM file creation
rule bwa_mem_to_bam:
    input:
        index = config["reference"]["index"],
        fasta = config["reference"]["fasta"],
        fq1_ID = "{sample}_R1."+config["fastqExtension"],
        fq2_ID = "{sample}_R2."+config["fastqExtension"],
        fq1 = join(config["fastqDir"], "{sample}_R1."+config["fastqExtension"]),
        fq2 = join(config["fastqDir"], "{sample}_R2."+config["fastqExtension"])
    output:
        temp(join(config["outputDir"], "{sample}.bamUnsorted"))
    version:
        subprocess.getoutput(
        "man bwa | tail -n 1 | cut -d ' ' -f 1 | cut -d '-' -f 2"
        )
    log:
        join(config["outputDir"], config["logDir"], "{sample}.bwa_mem.log")
    message:
        "Alignment of {input.fq1_ID} and {input.fq2_ID} on {input.fasta} with BWA version {version}."
    shell:
        "bwa mem {input.fasta} {input.fq1} {input.fq2} 2> {log} | samtools view -Sbh - > {output}"

# Sorting the BAM files on genomic positions
rule bam_sort:
    input:
        join(config["outputDir"], "{sample}.bamUnsorted")
    output:
        join(config["outputDir"], "{sample}.bam")
    log:
        join(config["outputDir"], config["logDir"],  "{sample}.samtools_sort.log")
    version:
        subprocess.getoutput(
            "samtools --version | "
            "head -1 | "
            "cut -d' ' -f2"
        )
    message:
        "Genomic sorting of {input} with samtools version {version}."
    shell:
        "samtools sort -f {input} {output} 2> {log}"

# Indexing the BAM files
rule bam_index:
    input:
        join(config["outputDir"], "{sample}.bam")
    output:
        join(config["outputDir"], "{sample}.bam.bai")
    message:
        "Indexing {input}."
    shell:
        "samtools index {input}"

我运行这个管道:

snakemake --cores 3 --snakefile /home/nico/labo/scripts/pipeline_illumina/snakefile_bwa_samtools.py --configfile /home/nico/labo/etudes/Optimal/data/snakemake_config_files/config_snakemake_Optimal_mapping_BaL.json

我有以下错误输出:

['sampleB', 'sampleA']
MissingInputException in line 18 of /home/nico/labo/scripts/pipeline_illumina/snakefile_bwa_samtools.py:
Missing input files for rule bwa_mem_to_bam:
sampleB_R1.fastq.gz
sampleB_R2.fastq.gz

或取决于时刻:

['sampleB', 'sampleA']
PeriodicWildcardError in line 40 of /home/nico/labo/scripts/pipeline_illumina/snakefile_bwa_samtools.py:
The value _unsorted in wildcard sample is periodically repeated (sampleB_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted_unsorted). This would lead to an infinite recursion. To avoid this, e.g. restrict the wildcards in this rule to certain values.

样本出现在列表中时被正确检测到(第一行输出类型),我肯定在规则中弄乱了通配符bwa_mem_to_bam,但我真的不明白为什么..有什么线索吗?

4

3 回答 3

1

我迅速查看了您的代码。

为什么第一个没有成功?

当您声明fq1_IDand时fq1,请查看示例 2。您没有分配相同的字符串。对于fq1您为 fq1_ID 不存在的文件添加一个目录,因此snakemake 正在工作目录(如果-d未设置选项,则为当前目录)中搜索带有您的字符串的文件名。因为这些变量在输入部分。

因此,通过删除两个 fq1/2_ID,它将清除所有文件搜索问题。

雨果

于 2017-03-07T04:49:01.123 回答
0

感谢 Hugo 检查我的代码和您的解释,这很有意义!我终于在今天早上醒来时突然想到(最好的),并意识到我忽略了params规则的一部分,fq1_ID 和 fq2_ID 不是输入而是参数。我将代码更改为:

rule bwa_mem_to_bam:
input:
    index = config["reference"]["index"],
    fasta = config["reference"]["fasta"],
    fq1 = join(config["fastqDir"], "{sample}_R1.fastq.gz"),
    fq2 = join(config["fastqDir"], "{sample}_R2.fastq.gz")
output:
    temp(join(config["outputDir"],"{sample}_unsorted.bam"))
params:
    fq1_ID = "{sample}_R1.fastq.gz",
    fq2_ID = "{sample}_R2.fastq.gz",
    ref_ID = os.path.basename(config["reference"]["fasta"])
version:
    subprocess.getoutput(
    "man bwa | tail -n 1 | cut -d ' ' -f 1 | cut -d '-' -f 2"
    )
log:
    join(config["outputDir"], config["logDir"], "{sample}.bwa_mem.log")
message:
    "Alignment of {params.fq1_ID} and {params.fq2_ID} on {params.ref_ID} with BWA version {version}."
shell:
    "bwa mem {input.fasta} {input.fq1} {input.fq2} 2> {log} | samtools view -Sbh - > {output}"

它工作得很好!

snakemake --cores 3 --snakefile /home/nico/labo/scripts/pipeline_illumina/snakefile_bwa_samtools.py --configfile /home/nico/labo/etudes/Optimal/data/snakemake_config_files/config_snakemake_Optimal_mapping_BaL.json 
Provided cores: 3
Rules claiming more threads will be scaled down.
Job counts:
    count   jobs
    1   all
    2   bam_index
    2   bam_sort
    2   bwa_mem_to_bam
    7
Alignment of sampleB_R1.fastq.gz and sampleB_R2.fastq.gz on BaL_AY713409.fasta with BWA version 0.7.12.

Alignment of sampleA_R1.fastq.gz and sampleA_R2.fastq.gz on BaL_AY713409.fasta with BWA version 0.7.12.

1 of 7 steps (14%) done
Genomic sorting of sampleB_unsorted.bam with samtools version 1.2.

Removing temporary output file /home/nico/labo/etudes/Optimal/data/mapping_BaL/sampleB_unsorted.bam.
2 of 7 steps (29%) done
Indexing sampleB.bam.

3 of 7 steps (43%) done
4 of 7 steps (57%) done
Genomic sorting of sampleA_unsorted.bam with samtools version 1.2.

Removing temporary output file /home/nico/labo/etudes/Optimal/data/mapping_BaL/sampleA_unsorted.bam.
5 of 7 steps (71%) done
Indexing sampleA.bam.

6 of 7 steps (86%) done
localrule all:
    input: /home/nico/labo/etudes/Optimal/data/mapping_BaL/sampleB.bam.bai, /home/nico/labo/etudes/Optimal/data/mapping_BaL/sampleA.bam.bai

7 of 7 steps (100%) done

最后得到我正确的信息:

  • 将 BaL_AY713409.fasta 上的 sampleB_R1.fastq.gz 和 sampleB_R2.fastq.gz 与 BWA 版本 0.7.12 对齐。
  • 将 BaL_AY713409.fasta 上的 sampleA_R1.fastq.gz 和 sampleA_R2.fastq.gz 与 BWA 版本 0.7.12 对齐。
于 2017-03-07T08:37:04.743 回答
0

最后,我成功地通过管道删除了规则中的fq1_IDandfq2_ID变量rule bwa_mem_to_bam并替换message了规则input.fq1_IDinput.fq2_IDby中的input.fq1and input.fq2

消息不太优雅,但管道运行正常。还是不明白到底哪里错了,如果有人能解释一下,我还在听!

的正确代码rule bwa_mem_to_bam

rule bwa_mem_to_bam:
input:
    index = config["reference"]["index"],
    fasta = config["reference"]["fasta"],
    fq1 = join(config["fastqDir"], "{sample}_R1."+config["fastqExtension"]),
    fq2 = join(config["fastqDir"], "{sample}_R2."+config["fastqExtension"])
output:
    temp(join(config["outputDir"], "{sample}.bamUnsorted"))
version:
    subprocess.getoutput(
    "man bwa | tail -n 1 | cut -d ' ' -f 1 | cut -d '-' -f 2"
    )
log:
    join(config["outputDir"], config["logDir"], "{sample}.bwa_mem.log")
message:
    "Alignment of {input.fq1} and {input.fq2} on {input.fasta} with BWA version {version}."
shell:
    "bwa mem {input.fasta} {input.fq1} {input.fq2} 2> {log} | samtools view -Sbh - > {output}"
于 2017-03-06T17:51:09.837 回答