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

Question

为了提供一些背景知识，我正在尝试建立一个管道来分析计算机预测 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 将引发错误。

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

谢谢， 爱德华

Answer 1

如果您需要将样本与 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}"

Answer 2

你显示的配置文件应该在连字符后有一个空格，即：

AMP:
  - ATAAAAACCATACAC...
  -GATGTGCTAGAG... <- should be "- GATGT..." 

但我想这不是问题。

AMP 和 gRNA 是列表，因此如果您想使用第一项，请使用 AMP=config['AMP'][0]。 （那么我猜你需要一些逻辑来将 fastq 文件对分配给正确的 AMP 和 gRNA。）

Answer 3

您在another answer下方评论了以下内容：

有没有办法为列表的所有元素表示位置 1:n 并在 shell 命令中循环

假设我了解您想要做什么，您可以执行以下操作：

rule crispresso_run:
    input:
        R1="data/{sample}_L001_R1_001.fastq.gz",
        ...
    output:
        ...
    params:
        AMP=config['AMP'],
        gRNA=config['gRNA']
    run:
        for i in range(0, len(params.AMP)):
            AMP = params.AMP[i]
            gRNA = params.gRNA[i]
            shell('/path/to/crispresso2 -a %s -g %s -r1 {input.R1}...' % (AMP, gRNA))

（但我怀疑如果您将通配符分配给 AMP 和 gRNA 并让 snakemake 处理循环，可能会有更好的解决方案）