标签: snakemake

我第一次使用snakemake 是为了使用cutadapt、bwa 和GATK（修剪；映射；调用）构建一个基本的管道。我想在目录中包含的每个 fastq 文件上运行此管道，而无需在蛇文件或配置文件中指定它们的名称或任何内容。我想成功地做到这一点。

前两个步骤（cutadapt 和 bwa/修剪和映射）运行良好，但我在使用 GATK 时遇到了一些问题。

首先，我必须从 bam 文件生成 g.vcf 文件。我正在使用以下规则执行此操作：

configfile: "config.yaml"

import os
import glob

rule all:
    input:
        "merge_calling.g.vcf"

rule cutadapt:
    input:
        read="data/Raw_reads/{sample}_R1_{run}.fastq.gz", 
        read2="data/Raw_reads/{sample}_R2_{run}.fastq.gz" 
    output:
        R1=temp("trimmed_reads/{sample}_R1_{run}.fastq.gz"),
        R2=temp("trimmed_reads/{sample}_R2_{run}.fastq.gz") 
    threads:
        10
    shell:
        "cutadapt -q {config[Cutadapt][Quality_value]} -m {config[Cutadapt][min_length]} -a {config[Cutadapt][forward_adapter]} -A  {config[Cutadapt][reverse_adapter]} -o {output.R1} -p '{output.R2}' {input.read} {input.read2}"

rule bwa_map:
    input:
        genome="data/genome.fasta",
        read=expand("trimmed_reads/{{sample}}_{pair}_{{run}}.fastq.gz", pair=["R1", "R2"]) 
    output:
        temp("mapped_bam/{sample}_{run}.bam")
    threads:
        10
    params:
        rg="@RG\\tID:{sample}\\tPL:ILLUMINA\\tSM:{sample}"
    shell:
        "bwa mem -t 2 -R '{params.rg}' {input.genome} {input.read} | samtools view -Sb - > {output}" …

我有一个规则，它遍历文件会拉出 Fastq 文件路径并在 Fastq 文件上运行 trimGalore。但是，某些文件已损坏/被截断，因此 trimGalore 无法处理它们。它继续在剩余文件上运行，但整体规则失败并删除包含成功处理文件的输出文件夹。如何保留输出文件夹？

我尝试更改 shell 命令以忽略退出状态，但 snakemake 似乎set -euo pipefail在运行的 shell 元素中强制执行。

rule trimGalore:
    """
    This module takes in the temporary file created by parse sampleFile rule and determines if libraries are single end or paired end.
    The appropriate step for trimGalore is then ran and a summary of the runs is produced in summary_tg.txt
    """
    input:
        rules.parse_sampleFile.output[1]+"singleFile.txt", rules.parse_sampleFile.output[1]+"pairFile.txt"
    output:
        directory(projectDir+"/trimmed_reads/")
    log:
        projectDir+"/logs/"+stamp+"_trimGalore.log"
    params:
        p = trimGaloreParams
    shell:
        """
        (awk -F "," '{{print …

我有一个奇怪的问题，它来来去去，我真的不知道什么时候以及为什么。
我正在运行这样的蛇形管道：

conda activate $myEnv    
snakemake -s $snakefile --configfile test.conf.yml --cluster "python $qsub_script" --latency-wait 60 --use-conda -p -j 10 --jobscript "$job_script"

我在 conda 环境中安装了 snakemake 5.9.1（也尝试降级到 5.5.4）。
如果我只运行这个命令，这工作正常，但是当我将此命令 qsub 到我正在使用的 PBS 集群时，我收到一个错误。我的 qsub 脚本如下所示：

#PBS stuff...

source ~/.bashrc
hostname
conda activate PGC_de_novo

cd $workDir
snakefile="..."
qsub_script="pbs_qsub_snakemake_wrapper.py"
job_script="..."
snakemake -s $snakefile --configfile test.conf.yml --cluster "python $qsub_script" --latency-wait 60 --use-conda -p -j 10 --jobscript "$job_script" >out 2>err

我得到的错误信息是：

...
Traceback (most recent call last):
  File "/path/to/pbs_qsub_snakemake_wrapper.py", line 6, in <module>
    from snakemake.utils …

我在建设有一个管道snakemake，并使用conda和singularity环境以确保可重复性。我遇到了一个错误，我的绘图上的文本被矩形替换

在 Linux 和 Mac 系统上试验管道并禁用奇异容器后，问题似乎源于缺少字体库，因为当我仅在 Mac 上仅使用conda( --use-conda)运行管道时，文本绘制正常。

奇点容器是从这个使用 Debian GNU/Linux 的miniconda docker镜像构建的。我设法创建了一个最小的示例管道，其中没有绘制文本。

# Snakefile
singularity: "docker://continuumio/miniconda3"

rule all:
    input:
        "mtcars-plot.png"

rule plot_mtcars:
    output:
        "mtcars-plot.png"
    conda:
        "minimal.yaml"
    script:
        "mtcars-test.R"

# mtcars-test.R
library(ggplot2)
png("mtcars-plot.png")
ggplot(mtcars, aes(factor(cyl), mpg)) + geom_boxplot()
dev.off()

# minimal.yaml
channels:
    - bioconda
    - conda-forge
    - defaults
dependencies:
    - r-base =3.6
    - r-ggplot2

要绘制断开的图，请运行管道

snakemake --use-conda --use-singularity

在 Debian GNU/Linux 上使用 R 正确绘制文本时，我可能缺少哪些软件包/库？

我有一个相当复杂的工作流程，有 750 个样本和大约 18.000 个作业，一开始 Snakemake 运行得很好，但在大约 4.000 个作业之后它突然冻结，重新启动后它会挂起“选择要执行的作业...”24 小时，之后我终止了它。不过，最初的 DAG 构建大约需要 2-3 分钟。

当我使用该--verbose选项运行 Snakemake（v5.32.0 和 v5.32.1）时，我得到大量与此类似的行：

Cbc0010I After 600 nodes, 304 on tree, -52534.791 best solution, best possible -52538.194 (7.08 seconds

我试图删除该.snakemake文件夹，希望那里会发生什么事情，但不幸的是，事实并非如此。对我来说，CBC MILP 求解器似乎在某种程度上没有收敛，并且它不断前进，并将最好的和最好的可能解决方案更紧密地结合在一起！？

现在我不知道如何继续并解决问题。我可能的解决方案是以某种方式改变收敛标准或求解器本身。在手册中我找到了该选项--scheduler-ilp-solver，但它显然只有一个选项，即默认选项COIN_CMD。

终止（较短的）运行后，我得到了这个详细的输出

Result - User ctrl-cuser ctrl-c

Objective value:                52534.79114334
Upper bound:                    52538.202
Gap:                            -0.00
Enumerated nodes:               186926
Total iterations:               1807277
Time (CPU seconds):             1181.97
Time (Wallclock seconds):       1188.11

接下来，我将尝试限制工作流程中的样本数量，看看这是否有任何影响（对于具有 500 个样本的其他数据集，它运行没有任何问题（使用 …

从文档中，我知道我可以定义一个配置文件并通过以下方式使用它

snakemake --profile <name>

我可以将这些配置文件之一设置为默认配置文件（在输入命令时使用snakemake）吗？

在 .bashrc 中定义别名是一种解决方法。然而，我想知道是否有“官方”snakemake 解决方案。

旁注：我对 Snakemake 默认配置文件的最初动机是定义默认的核心数量，这样我就不必总是键入

snakemake -j [cores]

默认核心数为“无”（对于我使用 conda 安装的 Snakemake）。

我有一个 Snakemake 配方，其中包含一个非常昂贵的准备步骤，对于所有调用来说都很常见。这是用于演示的伪规则：

rule sample:
    input:
        "{name}.config"
    output:
        "{name}.npz"
    run:
        import somemodule
        
        data = somemodule.Loader("some_big_data")  # expensive
        np.savez(output, data.process(input))  # also expensive

目前，每个目标都从头data加载，这不是最理想的。我怎样才能让它只加载一次？

rule sample:
    input:
        "{name}.config"
    output:
        "{name}.npz"
    setup:
        import somemodule
        
        data = somemodule.Loader("some_big_data")  # expensive
    run:
        np.savez(output, data.process(input))  # also expensive

rule sample:
    input:
        "{name}.config"
    output:
        "{name}.npz"
    run:
        import somemodule

        data = somemodule.Loader("some_big_data")  # expensive
        
        for job in jobs:
            np.savez(job.output,
                     data.process(job.input))  # also expensive