标签: snakemake

我不知道如何使用 Snakemake 规则来删除已变得无用的 Snakemake 输出文件。

具体来说，我有一条规则bwa_mem_sam创建一个名为{sample}.sam. 我有另一个规则，bwa_mem_bam它创建一个名为{sample.bam}.

如果这两个文件包含不同格式的相同信息，我想删除第一个文件，但无法成功执行此操作。

任何帮助将非常感激。本.

rule bwa_mem_map:
    input:
        sam="{sample}.sam",
        bam="{sample}.bam"
    shell:
        "rm {input.sam}"

# Convert SAM to BAM.
rule bwa_mem_map_bam:
    input:
        rules.sam_to_bam.output


# Use bwa mem to map reads on a reference genome.
rule bwa_mem_map_sam:
    input:
        reference=reference_genome(),
        index=reference_genome_index(),
        fastq=lambda wildcards: config["units"][SAMPLE_TO_UNIT[wildcards.sample]],
    output:
        "mapping/{sample}.sam"
    threads: 12
    log: 
        "mapping/{sample}.log"
    shell:
        "{BWA} mem -t {threads} {input.reference} {input.fastq} > {output} 2> {log} "\
        "|| (rc=$?; cat {log}; exit $rc;)"


rule sam_to_bam:
    input: …

据我了解蛇制作规则部分的文档log，必须“手动”将内容发送到日志文件。在我看来，使用本output节中定义的文件可以达到相同的结果。

这两种可能的方法之间的重要区别是什么？

该部分的真正用处是什么log？

我正在尝试使用 Snakemake 制作一个简单的管道，从网络上下载两个文件，然后将它们合并到一个输出中。

我认为可行的是以下代码：

dwn_lnks = {
    '1': 'https://molb7621.github.io/workshop/_downloads/sample.fa',
    '2': 'https://molb7621.github.io/workshop/_downloads/sample.fa'       
    }
import os

# association between chromosomes and their links
def chromo2link(wildcards):
    return dwn_lnks[wildcards.chromo]

rule all:
    input:
        os.path.join('genome_dir', 'human_en37_sm.fa')

rule download:
    output:
        expand(os.path.join('chr_dir', '{chromo}')),
    params:
        link=chromo2link,
    shell:
        "wget {params.link} -O {output}"


rule merger:
    input:
        expand(os.path.join('chr_dir', "{chromo}"), chromo=dwn_lnks.keys())
    output:
        os.path.join('genome_dir', 'human_en37_sm.fa')
    run:
        txt = open({output}, 'a+')
        with open (os.path.join('chr_dir', "{chromo}") as file:
                    line = file.readline()
                    while line:
                        txt.write(line)
                        line = file.readline()
        txt.close()

此代码返回错误： No values given for wildcard 'chromo'. in …

任何人都可以帮助我了解当样本名称未写入蛇形工作流程时是否可以从 config.yml 文件访问样本详细信息？这样我就可以为不同的项目重复使用工作流程，并且只调整配置文件。让我给你举个例子：

我有四个属于一起的样本，应该一起分析。它们被称为样本 1-4。每个样本都带有更多信息，但为了简单起见，我们可以说它只是一个名称标签，例如 S1、S2 等。

我的 config.yml 文件可能如下所示：

samples: ["sample1","sample2","sample3","sample4"]

sample1:
  tag: "S1"
sample2:
  tag: "S2"
sample3:
  tag: "S3"
sample4:
  tag: "S4"

这是我们使用的蛇文件的一个例子：

configfile: "config.yaml"

rule final: 
  input: expand("{sample}.txt", sample=config["samples"])

rule rule1:
  output:  "{sample}.txt"
  params:  tag=config["{sample}"]["tag"]
  shell:   """
           touch {output}
           echo {params.tag} > {output}

rule1 试图做的是创建一个以每个样本命名的文件samples，该文件保存在配置文件的变量中。到目前为止没有问题。然后，我想将示例标签打印到该文件中。正如上面编写的代码一样，运行snakemake将失败，因为config["{sample}"]实际上会{sample}在配置文件中查找不存在的变量，因为我需要将其替换为运行规则的当前示例，例如sample1.

有谁知道这是否有可能做到，如果是，我该怎么做？

理想情况下，我想进一步压缩信息（见下文），但这是更远的道路。

samples:
    sample1:
        tag: "S1"
    sample2:
        tag: "S2"
    sample3:
        tag: "S3"
    sample4:
        tag: "S4"

我使用 Snakemake 执行一些规则，但我遇到了一个问题：

rule filt_SJ_out:
input: 
    "pass1/{sample}SJ.out.tab"

output:
    "pass1/SJ.db"

shell:''' 
gawk '$6==1 || ($6==0 && $7>2)' {input} >> {output};


'''

在这里，我只想将一些文件合并到一个通用文件中，但是通过在谷歌上搜索，我发现输入中使用的通配符也必须在输出中使用。

但我找不到解决这个问题的解决方案..

提前致谢

这是我正在尝试做的事情的一个例子：

mydictionary={
'apple': 'crunchy fruit',
'banana': 'mushy and yellow'
}

rule all:
    input:
        expand('{key}.txt', key=mydictionary.keys())

rule test:
    output: temp('{f}.txt')
    shell:
        """
        echo {mydictionary[wildcards.f]} > {output}
        cat {output}
        """

由于某种原因，我无法访问字典内容。我尝试使用双大括号，但文本文件的内容变成了文字{mydictionary[wildcards.f]}（而我想要字典中相应条目的内容）。

我有大量这样组织的输入文件：

data/
??? set1/
?   ??? file1_R1.fq.gz
?   ??? file1_R2.fq.gz
?   ??? file2_R1.fq.gz
?   ??? file2_R2.fq.gz
|   :
?   ??? fileX_R2.fq.gz
??? another_set/
?   ??? asdf1_R1.fq.gz
?   ??? asdf1_R2.fq.gz
?   ??? asdf2_R1.fq.gz
?   ??? asdf2_R2.fq.gz
|   :
?   ??? asdfX_R2.fq.gz
:   
??? many_more_sets/
    ??? zxcv1_R1.fq.gz
    ??? zxcv1_R2.fq.gz
    :
    ??? zxcvX_R2.fq.gz

如果您熟悉生物信息学 - 这些当然是来自配对末端测序运行的 fastq 文件。我正在尝试生成一个可以读取所有这些的蛇形工作流程，但我已经在第一条规则上失败了。这是我的尝试：

configfile: "config.yaml"

rule all:
    input:
        read1=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R1.fq.gz", output=config["output"]),
        read2=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R2.fq.gz", output=config["output"])

rule clip_and_trim_reads:
    input:
        read1=expand("{data}/{set}/{{sample}}_R1.fq.gz", data=config["raw_data"], set=config["sets"]),
        read2=expand("{data}/{set}/{{sample}}_R2.fq.gz", data=config["raw_data"], set=config["sets"])
    output:
        read1=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R1.fq.gz", output=config["output"]),
        read2=expand("{output}/clipped_and_trimmed_reads/{{sample}}_R2.fq.gz", output=config["output"])
    threads: 16 …

我正在尝试采用一种简单的方法来在一个规则中创建工作流程所需的所有子目录。ChildIOException但是，每当我尝试执行在工作流程顶部创建所有必需目录的规则时，我都会得到一个对我来说毫无意义的信息：

Building DAG of jobs...
ChildIOException:
File/directory is a child to another output:
/scratch/groups/xxx/xxx/neand_sQTL/filtered_vcf
/scratch/groups/xxx/xxx/neand_sQTL/filtered_vcf/merged_filtered_chr1.vcf.gz

以下是有问题的规则：

rule mkdir_vcf:
    output:
        directory("gtex_vcf/"),
        directory("kg_vcf/"),
        directory("merged/"),
        directory("filtered_vcf/"),
        touch(".mkdir.chkpnt")
    shell:
        "mkdir -p {output}"

rule vcf_split1_23:
    input:
        vcf=config["vcf"],
        chk=".mkdir.chkpnt"
    output:
        "gtex_vcf/gtex_chr{i}.vcf"
    threads:
        23
    shell:
        "tabix -h {input.vcf} chr{wildcards.i} > {output}"

我尝试使用该directory()函数来查看是否有助于解决该错误，但没有。不知道在这里还能做什么。我不能包含mkdir其中，vcf_split1_23因为这是一项并行工作，制定一条成功创建目录一次但错误创建目录 22 次的规则是很糟糕的形式。我当然想mkdir_vcf在其他规则之前运行。

我处于一种情况，我想将我的工作流程分散到可变数量的块中，我事先不知道。也许通过具体来解释问题最容易：

有人递给我使用bcl2fastq该no-lane-splitting选项解复用的 FASTQ 文件。我想根据车道拆分这些文件，单独映射每个车道，然后最后再次收集所有内容。但是，我事先不知道车道数。

理想情况下，我想要这样的解决方案，

rule split_fastq_file: (...)  # results in N FASTQ files
rule map_fastq_file: (...)  # do this N times
rule merge_bam_files: (...)  # merge the N BAM files

但我不确定这是可能的。该expand功能要求我知道车道数，并且无法看到如何为此使用通配符。

我应该说我对 Snakemake 比较陌生，而且我可能完全误解了 Snakemake 的工作原理。我花了一些时间来习惯于通过关注输出文件然后倒退来思考“颠倒”的事情。

关于snakemake和使用多个配置文件的快速问题。

我正在为基因组预处理创建一个相当大的管道，其想法是它是用户友好的，并且在用户能力方面占据了最低的共同点。

因此，不允许用户在主配置文件中定义某些工具使用的线程。我想实现snakemake的“workflow.cores * Percentage”功能，即；我可以为每个工具指定总核心的百分比（如在命令行中使用 --cores # 定义的）。

这使事情变得简单，并且不会让用户对主配置文件中的所有选项感到困惑。但我仍然希望允许用户根据需要微调线程数量。我的想法是我有一个没有线程号的主配置文件，以及一个带有线程号的第二个配置文件。

如果用户决定使用自己决定的线程数，他们将简单地（在主配置中）在键中提供肯定的信息，如下所示：manualThreadChoice：“yes”

在 Snakefile 中有一个简单的 if 语句，如果 manualThreadChoice 是肯定的，则相应地向相应的规则提供线程数，否则仅使用自动定义的线程百分比。

根据biostars上的这个问题：https ://github.com/yanalab/celseq2/issues/33

现在可以定义多个配置文件。

见评论：

“感谢您报告此问题。这是因为snakemake更新了他们的API，现在支持多个配置文件。根据他们的日志，配置文件已更改为configfiles。请参阅：snakemake/snakemake@23624ee#diff-88e96378bf2405c8a8f8ac971519039e。”

因此，不要使用以下命令调用我们的配置文件

configfile: "path/to/config.yaml"

我们可以用：

configfiles: 

我的问题是，我们是否提供两个单独的路径作为文件列表的一部分：即

configfiles: ["path/to/config1.yaml", "path/to/config2.yaml"]

然后我们如何从不同的配置访问密钥。由于使用单个配置文件，我们将使用：

config['key']

我尝试通过索引访问不同的配置文件：

config[0]['key']

但这行不通。

我正在使用 Snakemake 5.7.0，我相信它是一个具有多个配置文件功能的版本。