时间戳

首页 > 解决方案 > 使用 Snakemake 在 Shell 命令中迭代字符串

snakemake - 使用 Snakemake 在 Shell 命令中迭代字符串

问题描述

为了提供一些背景知识，我正在尝试建立一个管道来分析计算机预测的 CRISPR 脱靶的深度测序结果。我在 50 个不同位置从基因组中扩增了一个已知序列，每个扩增子都包含一个预测的脱靶位点，我的原始 CRISPR 指南可能会在该位点结合。我将两个配对的末端 NGS 文件连同扩增子序列和脱靶指南一起输入到 CRISPResso 程序中。50 个位点的两个配对末端文件、扩增子序列和 gRNA 各不相同。我需要使用多个条件、捐助者和复制品来做到这一点，以便数字快速加起来。

我在下面创建了蛇形工作流程：

configfile: "config.yaml"

singularity:
    "docker://pinellolab/crispresso2"

SAMPLES, =glob_wildcards("data/{sample}_L001_R1_001.fastq.gz")

rule all:
    input:
        expand("outfiles/{sample}", sample=SAMPLES)

rule crispresso_run:
    input:
        R1="data/{sample}_L001_R1_001.fastq.gz",
        R2="data/{sample}_L001_R2_001.fastq.gz"

    params:
        AMP=config['AMP'],
        gRNA=config['gRNA']

    output:
        directory("outfiles/{sample}")

    shell:
        "docker run -v ${{PWD}}:/DATA -w /DATA -i pinellolab/crispresso2 \
            CRISPResso -r1 {input.R1} \
            -r2 {input.R2} \
            -a {params.AMP} \
            -g {params.gRNA} \
            -q 30 -s 30 --min_bp_quality_or_N 30 -w 3 -o {output}"

使用配置文件：

AMP:
  - ATAAAAACCATACACATTCAGTGGGAAACCTTCAGCCATAGAGAAGTATAGGCAGGGTGCAGCTGATTGCTCTGTCTTTGGGCAATTTAGCTTTTAGGCCAGAGGCCACAGATGGGTAGCCTGGTGTGTGCCTAGGGTGTTTTTGTTTGGCTGGCGCAATATTTTTTAAAACTGTAAGTTTATTGCCAGCATTTAA
  -GATGTGCTAGAGATGAGAAAGGATGTGGCAGAAGAAGTACCTATCTCTTGAGGGATGAAGTGGCCTCATTTCACCTACTGAGAGTCAGGAAGTGCCCCATCTGCAGCCTCTGGGCTGGTTGGGGTCAGTCTGCAGATTTTCCTTGCTTTCTCCCATGCCCTTGTCTTTCTCTCCCCTGTAGAGAAAGACACTGATGTTGCTGTTGTTCTAGGAACAGTGGAGACAACTG
gRNA:
  - TTAGGCCAGAGGCCACAGATGGG
  - CCAGCCCAGAGGCTGCAGATGGG

shell 命令中的 -a 和 -g 参数分别是应该放置扩增子和 gRNA 的位置。

当我在配置文件中只有一个扩增子和 gRNA 时，一切都很好。但是当我在配置文件中包含两个时，snakemake 只是连接两个扩增子或 gRNA，这会引发错误。我猜我错误地使用了 params 函数，但是找不到任何显示如何以这种方式使用它的东西（即在 shell 命令中迭代字符串）。

此外，为了澄清，配对的末端文件需要以与每个相应的 gRNA 和扩增子对相同的顺序运行，否则 CRISPResso 将引发错误。

我对这一切都很陌生，对自己解决问题的了解还不够。任何帮助将不胜感激。

谢谢，爱德华

标签： snakemake

解决方案

如果您需要将样本与 AMP 和 gRNA 序列配对，您需要了解哪些样本与哪些样本配对的信息，因此您需要知道哪些样本名称是可能的。

假设样本名为“sample1”和“sample2”，您将更改配置文件如下（未经测试）：

{
    "AMP" : {
        "sample1" : "ATAAAAACCATACACATTCAGTGGGAAACCTTCAGCCATAGAGAAGTATAGGCAGGGTGCAGCTGATTGCTCTGTCTTTGGGCAATTTAGCTTTTAGGCCAGAGGCCACAGATGGGTAGCCTGGTGTGTGCCTAGGGTGTTTTTGTTTGGCTGGCGCAATATTTTTTAAAACTGTAAGTTTATTGCCAGCATTTAA",
        "sample2" : "GATGTGCTAGAGATGAGAAAGGATGTGGCAGAAGAAGTACCTATCTCTTGAGGGATGAAGTGGCCTCATTTCACCTACTGAGAGTCAGGAAGTGCCCCATCTGCAGCCTCTGGGCTGGTTGGGGTCAGTCTGCAGATTTTCCTTGCTTTCTCCCATGCCCTTGTCTTTCTCTCCCCTGTAGAGAAAGACACTGATGTTGCTGTTGTTCTAGGAACAGTGGAGACAACTG"
    },
    "gRNA" : {
        "sample1" : "TTAGGCCAGAGGCCACAGATGGG",
        "sample2" : "CCAGCCCAGAGGCTGCAGATGGG"
    }
}

以下是如何修改您的蛇文件以使用此配对信息：

configfile: "config.yaml"

singularity:
    "docker://pinellolab/crispresso2"

# SAMPLES, =glob_wildcards("data/{sample}_L001_R1_001.fastq.gz")
SAMPLES = config["AMP"].keys()

rule all:
    input:
        expand("outfiles/{sample}", sample=SAMPLES)

def get_AMP(wildcards):
    return config["AMP"][wildcards.sample]

def get_gRNA(wildcards):
    return config["gRNA"][wildcards.sample]

rule crispresso_run:
    input:
        R1="data/{sample}_L001_R1_001.fastq.gz",
        R2="data/{sample}_L001_R2_001.fastq.gz"

    params:
        AMP=get_AMP,
        gRNA=get_gRNA

    output:
        directory("outfiles/{sample}")

    shell:
        "docker run -v ${{PWD}}:/DATA -w /DATA -i pinellolab/crispresso2 \
            CRISPResso -r1 {input.R1} \
            -r2 {input.R2} \
            -a {params.AMP} \
            -g {params.gRNA} \
            -q 30 -s 30 --min_bp_quality_or_N 30 -w 3 -o {output}"

推荐阅读